Merge pull request #290 from Ao-chuba/fix/cubic-unit-detection

Fix cubic file unit detection for negative dimensions

Author: marco-2023

src/grid/cubic.py

@@ -25,6 +25,7 @@
 from sympy.functions.combinatorial.numbers import bell
 
 from grid.basegrid import Grid, OneDGrid
+from .utils import ANGSTROM_TO_BOHR
 
 
 class _HyperRectangleGrid(Grid):
@@ -701,7 +702,6 @@ def from_cube(cls, fname, weight="Trapezoid", return_data=False):
             - ``atcorenums``\: Pseudo-number of :math:`M` atoms in the molecule.
             - ``atcoords``\: Cartesian coordinates of :math:`M` atoms in the molecule.
             - ``data``\: the grid data stored in a flattened one-dimensional array.
-            - ``unit``\: the unit of the grid (either "angstrom" or "bohr").
 
         """
         fname = str(fname)
@@ -715,44 +715,37 @@ def from_cube(cls, fname, weight="Trapezoid", return_data=False):
 
             def read_grid_line(line):
                 """Read a number and (x, y, z) coordinate from the cube file line."""
-                words = line.split()
-                return (
-                    int(words[0]),
-                    np.array([float(words[1]), float(words[2]), float(words[3])], float),
-                    # all coordinates in a cube file are in atomic units
-                )
+                npts, *vec = line.split()
+                return int(npts), np.asarray(vec, dtype=float)
 
             # number of atoms and origin of the grid
             natom, origin = read_grid_line(f.readline())
             # number of grid points in A direction and step vector A, and so on
             shape0, axis0 = read_grid_line(f.readline())
             shape1, axis1 = read_grid_line(f.readline())
             shape2, axis2 = read_grid_line(f.readline())
-            axes = np.array([axis0, axis1, axis2], dtype=float)
-            # if shape0, shape1, shape2 are negative, the units are in bohr
-            # otherwise in angstrom
+            axes = np.asarray([axis0, axis1, axis2], dtype=float)
+            # if shape0 is negative the units are in angstroms otherwise atomic units.
+            # this was verified with cube files generated from Gaussian on february 2026.
             # https://gaussian.com/cubegen/
-            unit = "bohr" if shape0 < 0 else "angstrom"
-            # print out the units detected for clarity
-            print(f"Cube file units detected: {unit}")
+            coordinates_in_angstrom = shape0 < 0
 
             # Convert negative shape values to positive
-            shape = np.array([abs(shape0), abs(shape1), abs(shape2)], dtype=int)
+            shape = np.abs(np.asarray([shape0, shape1, shape2], dtype=int))
+            if coordinates_in_angstrom:
+                print(f"Cube file units detected: angstrom, converting to atomic units grid.")
+                axes *= ANGSTROM_TO_BOHR
+                origin *= ANGSTROM_TO_BOHR
 
             # if return_data=False, only grid is returned
             if not return_data:
                 return cls(origin, axes, shape, weight)
 
             # otherwise, return the atomic numbers, coordinates, and the grid data as well
             def read_coordinate_line(line):
-                """Read atomic number and (x, y, z) coordinate from the cube file line."""
-                words = line.split()
-                return (
-                    int(words[0]),
-                    float(words[1]),
-                    np.array([float(words[2]), float(words[3]), float(words[4])], float),
-                    # all coordinates in a cube file are in atomic units
-                )
+                """Read atomic number, charge, and (x, y, z) coordinates from cube file line."""
+                atnum, charge, *coords = line.split()
+                return int(atnum), float(charge), np.asarray(coords, dtype=float)
 
             numbers = np.zeros(natom, int)
             # get the core charges
@@ -765,6 +758,9 @@ def read_coordinate_line(line):
                 if pseudo_numbers[i] == 0.0:
                     pseudo_numbers[i] = numbers[i]
 
+            if coordinates_in_angstrom:
+                coordinates *= ANGSTROM_TO_BOHR
+
             # load data stored in the cube file
             data = np.zeros(tuple(shape), float).ravel()
             counter = 0
@@ -781,7 +777,6 @@ def read_coordinate_line(line):
                 "atcorenums": pseudo_numbers,
                 "atcoords": coordinates,
                 "data": data,
-                "unit": unit,
             }
         return cls(origin, axes, shape, weight), cube_data
 

src/grid/data/tests/cubegen_ch4_6_gen_negative.cube

@@ -1,14 +1,14 @@
 Cubefile created with THEOCHEM Grid
 OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z
-    5   -6.512793   -6.512718   -6.512750
-   -6    2.605101    0.000000    0.000000
-   -6    0.000000    2.605101    0.000000
-   -6    0.000000    0.000000    2.605101
-    6    6.000000   -0.000041    0.000035    0.000003
-    1    1.000000    1.051607    1.517467    0.942259
-    1    1.000000    1.296786   -1.543671   -0.481272
-    1    1.000000   -1.476881   -0.708836    1.269987
-    1    1.000000   -0.871678    0.735176   -1.730958
+    5   -3.446421   -3.446421   -3.446421
+   -6    1.378560    0.000000    0.000000
+   -6    0.000000    1.378560    0.000000
+   -6    0.000000    0.000000    1.378560
+    6    6.000000   -0.000021    0.000018    0.000001
+    1    1.000000    0.556486    0.803009    0.498622
+    1    1.000000    0.686229   -0.816875   -0.254678
+    1    1.000000   -0.781531   -0.375099    0.672048
+    1    1.000000   -0.461272    0.389038   -0.915983
   2.66011E-09  2.43339E-09  2.14921E-08  7.16256E-08  7.41998E-08  3.40400E-08
   1.69019E-08  2.07614E-08  2.34886E-07  9.17704E-07  4.01950E-07  1.05928E-07
   7.35455E-08  3.21547E-07  2.09291E-06  6.03301E-06  1.45834E-06  1.68969E-07

src/grid/tests/test_cubic.py

@@ -27,6 +27,7 @@
 
 from grid.cubic import Tensor1DGrids, UniformGrid, _HyperRectangleGrid
 from grid.onedgrid import GaussLaguerre, MidPoint
+from grid.utils import ANGSTROM_TO_BOHR
 
 
 class TestHyperRectangleGrid(TestCase):
@@ -1246,20 +1247,20 @@ def test_uniformgrid_from_cube_negative(self):
 
         ref_grid = UniformGrid(origin, axes, shape)
 
+        # the precision of the numbers in this file is lower than reference, large atol needed
         cubefile = files("grid") / "data" / "tests" / "cubegen_ch4_6_gen_negative.cube"
         grid, cube_data = UniformGrid.from_cube(cubefile, return_data=True)
 
-        assert_allclose(grid._axes, axes)
-        assert_allclose(grid._origin, origin)
+        assert_allclose(grid._axes, axes, atol=1e-4)
+        assert_allclose(grid._origin, origin, atol=1e-4)
         assert_allclose(grid._shape, shape)
         assert_allclose(cube_data["atnums"], atnums)
-        assert_allclose(cube_data["atcoords"], atcoords)
+        assert_allclose(cube_data["atcoords"], atcoords, atol=1e-4)
         assert_allclose(cube_data["atcorenums"], pseudo_numbers)
         assert_allclose(cube_data["data"], data_vals)
-        assert_equal(cube_data["unit"], "bohr")
 
-        assert_allclose(grid.points, ref_grid.points)
-        assert_allclose(grid.weights, ref_grid.weights)
+        assert_allclose(grid.points, ref_grid.points, atol=1e-4)
+        assert_allclose(grid.weights, ref_grid.weights, atol=1e-6)
 
     def test_uniformgrid_generate_cube(self):
         r"""Test creating uniform cubic grid from cube example."""

src/grid/utils.py

@@ -21,6 +21,11 @@
 
 import numpy as np
 from scipy.special import sph_harm_y
+from scipy.constants import physical_constants, angstrom
+
+_BOHR_RADIUS_M = physical_constants["Bohr radius"][0]
+BOHR_TO_ANGSTROM = _BOHR_RADIUS_M / angstrom
+ANGSTROM_TO_BOHR = 1.0 / BOHR_TO_ANGSTROM
 
 _bragg = np.array(
     [