Fix PBC DFT Becke grid rounding error (#1960)

pyscf · Nov 16, 2023 · e64ed31 · e64ed31
1 parent fd8ee3e
commit e64ed31
Show file tree

Hide file tree

Showing 2 changed files with 36 additions and 9 deletions.
diff --git a/pyscf/pbc/dft/gen_grid.py b/pyscf/pbc/dft/gen_grid.py
@@ -165,6 +165,7 @@ def get_becke_grids(cell, atom_grid={}, radi_method=dft.radi.gauss_chebyshev,
     b = cell.reciprocal_vectors(norm_to=1)
     supatm_idx = []
     k = 0
+    tol = 1e-15
     for iL, L in enumerate(Ls):
         for ia in range(cell.natm):
             coords, vol = atom_grids_tab[cell.atom_symbol(ia)]
@@ -174,24 +175,24 @@ def get_becke_grids(cell, atom_grid={}, radi_method=dft.radi.gauss_chebyshev,
 
             mask = np.ones(c.shape[1], dtype=bool)
             if dimension >= 1:
-                mask &= (c[0]>=-.5) & (c[0]<=.5)
+                mask &= (c[0]>-.5-tol) & (c[0]<.5+tol)
             if dimension >= 2:
-                mask &= (c[1]>=-.5) & (c[1]<=.5)
+                mask &= (c[1]>-.5-tol) & (c[1]<.5+tol)
             if dimension == 3:
-                mask &= (c[2]>=-.5) & (c[2]<=.5)
+                mask &= (c[2]>-.5-tol) & (c[2]<.5+tol)
 
             vol = vol[mask]
             if vol.size > 8:
                 c = c[:,mask]
                 if dimension >= 1:
-                    vol[c[0]==-.5] *= .5
-                    vol[c[0]== .5] *= .5
+                    vol[abs(c[0]+.5) < tol] *= .5
+                    vol[abs(c[0]-.5) < tol] *= .5
                 if dimension >= 2:
-                    vol[c[1]==-.5] *= .5
-                    vol[c[1]== .5] *= .5
+                    vol[abs(c[1]+.5) < tol] *= .5
+                    vol[abs(c[1]-.5) < tol] *= .5
                 if dimension == 3:
-                    vol[c[2]==-.5] *= .5
-                    vol[c[2]== .5] *= .5
+                    vol[abs(c[2]+.5) < tol] *= .5
+                    vol[abs(c[2]-.5) < tol] *= .5
                 coords = coords[mask]
                 coords_all.append(coords)
                 weights_all.append(vol)

diff --git a/pyscf/pbc/dft/test/test_gen_grid.py b/pyscf/pbc/dft/test/test_gen_grid.py
@@ -125,6 +125,32 @@ def test_becke_grids_2d_low_dim_ft_type(self):
         s1 = get_ovlp(cell, grids)
         self.assertAlmostEqual(numpy.linalg.norm(s1-s2), 0, 5)
 
+    def test_becke_grids_round_error(self):
+        cell = pgto.Cell()
+        cell.atom = '''
+            Cu1    1.927800000000   1.927800000000   1.590250000000
+            Cu1    5.783400000000   5.783400000000   1.590250000000
+            Cu2    1.927800000000   5.783400000000   1.590250000000
+            Cu2    5.783400000000   1.927800000000   1.590250000000
+            O1     1.927800000000   3.855600000000   1.590250000000
+            O1     3.855600000000   5.783400000000   1.590250000000
+            O1     5.783400000000   3.855600000000   1.590250000000
+            O1     3.855600000000   1.927800000000   1.590250000000
+            O1     0.000000000000   1.927800000000   1.590250000000
+            O1     1.927800000000   7.711200000000   1.590250000000
+            O1     7.711200000000   5.783400000000   1.590250000000
+            O1     5.783400000000   0.000000000000   1.590250000000
+        '''
+        cell.a = numpy.array(
+            [[7.7112,    0.,    0.],
+             [0.,    7.7112,    0.],
+             [0.,    0.,    3.1805]])
+        cell.build()
+
+        grids = gen_grid.BeckeGrids(cell)
+        grids.level = 1
+        grids.build()
+        assert abs(cell.vol - grids.weights.sum()) < .1
 
 if __name__ == '__main__':
     print("Full Tests for Becke grids")