From 31ea02b26193a919e052887205bbca93523d3a7c Mon Sep 17 00:00:00 2001 From: aochuba Date: Wed, 4 Feb 2026 22:21:07 +0530 Subject: [PATCH 1/5] Fix cubic file unit detection for negative dimensions --- src/grid/cubic.py | 10 ++++++---- src/grid/tests/test_cubic.py | 38 +++++++++++++++++++++++++++++++++++- 2 files changed, 43 insertions(+), 5 deletions(-) diff --git a/src/grid/cubic.py b/src/grid/cubic.py index 6d1a389b..d79459d3 100644 --- a/src/grid/cubic.py +++ b/src/grid/cubic.py @@ -729,10 +729,12 @@ def read_grid_line(line): 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 - # https://gaussian.com/cubegen/ - unit = "bohr" if shape0 < 0 else "angstrom" + # if shape0, shape1, shape2 are negative, the units are in angstrom + # otherwise in bohr. + # See https://gaussian.com/cubegen/ + # "Negative values of npts <=-5 specify spacing of npts*10^-3 Angstroms" + # However, typically negative simply acts as a flag for Angstroms in many parsers/tools. + unit = "angstrom" if shape0 < 0 else "bohr" # print out the units detected for clarity print(f"Cube file units detected: {unit}") diff --git a/src/grid/tests/test_cubic.py b/src/grid/tests/test_cubic.py index 1a0666f1..4ab08461 100644 --- a/src/grid/tests/test_cubic.py +++ b/src/grid/tests/test_cubic.py @@ -1256,7 +1256,7 @@ def test_uniformgrid_from_cube_negative(self): assert_allclose(cube_data["atcoords"], atcoords) assert_allclose(cube_data["atcorenums"], pseudo_numbers) assert_allclose(cube_data["data"], data_vals) - assert_equal(cube_data["unit"], "bohr") + assert_equal(cube_data["unit"], "angstrom") assert_allclose(grid.points, ref_grid.points) assert_allclose(grid.weights, ref_grid.weights) @@ -1310,3 +1310,39 @@ def test_uniformgrid_points_without_rotate(self): ] ) assert_allclose(grid.points, expected, rtol=1.0e-7, atol=1.0e-7) + + def test_cube_unit_detection(self): + r"""Test that negative vs positive shape in Cube files determines units (Angstrom vs Bohr).""" + import os + filename = "test_unit_detection_temp.cube" + + def create_dummy_cube(fname, shape_sign=1): + with open(fname, 'w') as f: + f.write("Dummy Cube File\n") + f.write("OUTER LOOP\n") + f.write(f" 1 0.000000 0.000000 0.000000\n") # Natoms, Origin + val = 10 * shape_sign + f.write(f" {val} 0.100000 0.000000 0.000000\n") # N1, Axis1 + f.write(f" {val} 0.000000 0.100000 0.000000\n") # N2, Axis2 + f.write(f" {val} 0.000000 0.000000 0.100000\n") # N3, Axis3 + f.write(" 1 0.000000 0.000000 0.000000 0.000000\n") # Atom 1 + for _ in range(10*10*10): # Data + f.write(" 0.0") + f.write("\n") + + try: + # Test 1: Positive Shape -> Bohr + create_dummy_cube(filename, shape_sign=1) + grid, data = UniformGrid.from_cube(filename, return_data=True) + assert_equal(data['unit'], 'bohr') + assert_equal(grid.shape, [10, 10, 10]) + + # Test 2: Negative Shape -> Angstrom + create_dummy_cube(filename, shape_sign=-1) + grid, data = UniformGrid.from_cube(filename, return_data=True) + assert_equal(data['unit'], 'angstrom') + assert_equal(grid.shape, [10, 10, 10]) + + finally: + if os.path.exists(filename): + os.remove(filename) From f8b183a8c43335989953dd11535a0c1e4102bb41 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Thu, 5 Feb 2026 10:31:51 -0500 Subject: [PATCH 2/5] Improve UniformGrid parsing of cube files - Apply Angstrom to Bohr conversion only to grid origin and axes - Ensure atomic coordinates are converted consistently when flagged - Improve robustness and clarity of cube atom line parsing - Remove unused unit metadata from cube data --- src/grid/cubic.py | 41 ++++++++++++++++++----------------------- src/grid/utils.py | 5 +++++ 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/src/grid/cubic.py b/src/grid/cubic.py index 6d1a389b..628035a3 100644 --- a/src/grid/cubic.py +++ b/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,12 +715,8 @@ 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()) @@ -728,16 +724,18 @@ def read_grid_line(line): 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: @@ -745,14 +743,9 @@ def read_grid_line(line): # 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 diff --git a/src/grid/utils.py b/src/grid/utils.py index dd2ffc28..2f77b092 100644 --- a/src/grid/utils.py +++ b/src/grid/utils.py @@ -21,6 +21,11 @@ import numpy as np from scipy.special import sph_harm +from scipy.constants import physical_constants, angstrom + +_BOHR_RADIUS_M = physical_constants["Bohr radius"][0] +BOHR_TO_ANGSTROM = angstrom / _BOHR_RADIUS_M +ANGSTROM_TO_BOHR = 1.0 / BOHR_TO_ANGSTROM _bragg = np.array( [ From 2b95080ca8af519c733bd0295e654409273fa97a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Thu, 5 Feb 2026 10:57:23 -0500 Subject: [PATCH 3/5] Remove cube_unit_detection test that is not longer needed --- src/grid/tests/test_cubic.py | 38 +----------------------------------- 1 file changed, 1 insertion(+), 37 deletions(-) diff --git a/src/grid/tests/test_cubic.py b/src/grid/tests/test_cubic.py index 4ab08461..4a28df9f 100644 --- a/src/grid/tests/test_cubic.py +++ b/src/grid/tests/test_cubic.py @@ -1309,40 +1309,4 @@ def test_uniformgrid_points_without_rotate(self): [1, 0, 0], ] ) - assert_allclose(grid.points, expected, rtol=1.0e-7, atol=1.0e-7) - - def test_cube_unit_detection(self): - r"""Test that negative vs positive shape in Cube files determines units (Angstrom vs Bohr).""" - import os - filename = "test_unit_detection_temp.cube" - - def create_dummy_cube(fname, shape_sign=1): - with open(fname, 'w') as f: - f.write("Dummy Cube File\n") - f.write("OUTER LOOP\n") - f.write(f" 1 0.000000 0.000000 0.000000\n") # Natoms, Origin - val = 10 * shape_sign - f.write(f" {val} 0.100000 0.000000 0.000000\n") # N1, Axis1 - f.write(f" {val} 0.000000 0.100000 0.000000\n") # N2, Axis2 - f.write(f" {val} 0.000000 0.000000 0.100000\n") # N3, Axis3 - f.write(" 1 0.000000 0.000000 0.000000 0.000000\n") # Atom 1 - for _ in range(10*10*10): # Data - f.write(" 0.0") - f.write("\n") - - try: - # Test 1: Positive Shape -> Bohr - create_dummy_cube(filename, shape_sign=1) - grid, data = UniformGrid.from_cube(filename, return_data=True) - assert_equal(data['unit'], 'bohr') - assert_equal(grid.shape, [10, 10, 10]) - - # Test 2: Negative Shape -> Angstrom - create_dummy_cube(filename, shape_sign=-1) - grid, data = UniformGrid.from_cube(filename, return_data=True) - assert_equal(data['unit'], 'angstrom') - assert_equal(grid.shape, [10, 10, 10]) - - finally: - if os.path.exists(filename): - os.remove(filename) + assert_allclose(grid.points, expected, rtol=1.0e-7, atol=1.0e-7) \ No newline at end of file From 199d365fb1906c5cf62baa73a0adcf436ebb619e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Thu, 5 Feb 2026 11:39:17 -0500 Subject: [PATCH 4/5] Fix UniformGrid test for cube loading with negative eigenvalues - Convert cube coordinates, axes, and origin to Angstroms - Adjust test tolerances to account for lower cube file precision --- .../data/tests/cubegen_ch4_6_gen_negative.cube | 18 +++++++++--------- src/grid/tests/test_cubic.py | 15 ++++++++------- 2 files changed, 17 insertions(+), 16 deletions(-) diff --git a/src/grid/data/tests/cubegen_ch4_6_gen_negative.cube b/src/grid/data/tests/cubegen_ch4_6_gen_negative.cube index 0a6c502b..6b611c44 100644 --- a/src/grid/data/tests/cubegen_ch4_6_gen_negative.cube +++ b/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 diff --git a/src/grid/tests/test_cubic.py b/src/grid/tests/test_cubic.py index 4a28df9f..026a981a 100644 --- a/src/grid/tests/test_cubic.py +++ b/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"], "angstrom") - 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.""" @@ -1309,4 +1310,4 @@ def test_uniformgrid_points_without_rotate(self): [1, 0, 0], ] ) - assert_allclose(grid.points, expected, rtol=1.0e-7, atol=1.0e-7) \ No newline at end of file + assert_allclose(grid.points, expected, rtol=1.0e-7, atol=1.0e-7) From b58237630c1ea07d5576a98bbecfe3876635a102 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Thu, 5 Feb 2026 11:50:37 -0500 Subject: [PATCH 5/5] Fix wrong BOHR_TO_ANGSTROM unit conversion --- src/grid/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/grid/utils.py b/src/grid/utils.py index 2f77b092..b7647e08 100644 --- a/src/grid/utils.py +++ b/src/grid/utils.py @@ -24,7 +24,7 @@ from scipy.constants import physical_constants, angstrom _BOHR_RADIUS_M = physical_constants["Bohr radius"][0] -BOHR_TO_ANGSTROM = angstrom / _BOHR_RADIUS_M +BOHR_TO_ANGSTROM = _BOHR_RADIUS_M / angstrom ANGSTROM_TO_BOHR = 1.0 / BOHR_TO_ANGSTROM _bragg = np.array(