Fix the bad atomic radii assigned to ghost atom, which causes a bad Becke partition and a negative BSSE

gpu4pyscf/dft/gen_grid.py

@@ -35,6 +35,8 @@
 from pyscf.dft import gen_grid as gen_grid_cpu
 from gpu4pyscf.lib import utils
 from pyscf.gto.eval_gto import BLKSIZE, NBINS, CUTOFF, make_screen_index
+from pyscf.data.elements import _std_symbol_without_ghost
+from pyscf.data.elements import charge as elements_proton
 from pyscf import __config__
 from gpu4pyscf.lib import logger
 from gpu4pyscf.dft import radi
@@ -217,7 +219,7 @@ def gen_atomic_grids(mol, atom_grid={}, radi_method=radi.gauss_chebyshev,
         symb = mol.atom_symbol(ia)
 
         if symb not in atom_grids_tab:
-            chg = gto.charge(symb)
+            chg = elements_proton(_std_symbol_without_ghost(symb))
             atom_config = atom_grid.get(symb, default)
             if atom_config is not None:
                 n_rad, n_ang = atom_config

gpu4pyscf/dft/radi.py

@@ -22,6 +22,7 @@
 import pyscf
 from pyscf.data import radii
 from pyscf.data.elements import charge as elements_proton
+from pyscf.data.elements import _std_symbol_without_ghost
 
 BRAGG_RADII = radii.BRAGG
 COVALENT_RADII = radii.COVALENT
@@ -35,7 +36,7 @@ def treutler_atomic_radii_adjust(mol, atomic_radii):
 # i > j
 # fac(i,j) = \frac{1}{4} ( \frac{ra(j)}{ra(i)} - \frac{ra(i)}{ra(j)}
 # fac(j,i) = -fac(i,j)
-    charges = [elements_proton(x) for x in mol.elements]
+    charges = [elements_proton(_std_symbol_without_ghost(x)) for x in mol.elements]
     rad = cupy.sqrt(atomic_radii[charges]) + 1e-200
     rr = rad.reshape(-1,1) * (1./rad)
     a = .25 * (rr.T - rr)
@@ -55,7 +56,7 @@ def get_treutler_fac(mol, atomic_radii):
     # fac(i,j) = \frac{1}{4} ( \frac{ra(j)}{ra(i)} - \frac{ra(i)}{ra(j)}
     # fac(j,i) = -fac(i,j)
     '''
-    charges = [elements_proton(x) for x in mol.elements]
+    charges = [elements_proton(_std_symbol_without_ghost(x)) for x in mol.elements]
     #atomic_radii = cupy.asarray(atomic_radii[charges])
     rad = numpy.sqrt(atomic_radii[charges]) + 1e-200
     rr = rad.reshape(-1,1) * (1./rad)
@@ -65,7 +66,7 @@ def get_treutler_fac(mol, atomic_radii):
     return cupy.asarray(a)
 
 def get_becke_fac(mol, atomic_radii):
-    charges = [elements_proton(x) for x in mol.elements]
+    charges = [elements_proton(_std_symbol_without_ghost(x)) for x in mol.elements]
     atomic_radii = numpy.asarray(atomic_radii[charges])
     rad = atomic_radii[charges] + 1e-200
     rr = rad.reshape(-1,1) * (1./rad)

gpu4pyscf/dft/tests/test_rks.py

@@ -194,6 +194,63 @@ def test_roks(self):
         ref = mf.to_cpu().run()
         self.assertAlmostEqual(mf.e_tot, ref.e_tot, 8)
 
+    def test_ghost_atom_close_to_real_atom(self):
+        mol = pyscf.M(
+            atom = """
+            O      0.000000     0.000000     0.000000
+            H      0.960000     0.000000     0.000000
+            H     -0.240000     0.930000     0.000000
+            ghost:H     -0.240000    -0.310000     0.880000
+            """,
+            basis = "def2-TZVPD", # Large and diffuse basis required to reproduce the bug
+            verbose = 0,
+        )
+
+        mf = rks.RKS(mol, xc = "pbe").density_fit(auxbasis = "def2-universal-jkfit")
+        mf.grids.atom_grid = (99,590)
+        mf.conv_tol = 1e-12
+
+        test_energy = mf.kernel()
+
+        gobj = mf.Gradients()
+        gobj.grid_response = True
+        test_gradient = gobj.kernel()
+
+        ### Q-Chem reference input
+        # $molecule
+        # 0 1
+        #     O      0.000000     0.000000     0.000000
+        #     H      0.960000     0.000000     0.000000
+        #     H     -0.240000     0.930000     0.000000
+        #     @H     -0.240000    -0.310000     0.880000
+        # $end
+
+        # $rem
+        # JOBTYPE force
+        # METHOD PBE
+        # XC_GRID       000099000590
+        # BASIS def2-TZVPD
+        # SYMMETRY      FALSE
+        # SYM_IGNORE    TRUE
+        # MAX_SCF_CYCLES 100
+        # PURECART 1111
+        # SCF_CONVERGENCE 10
+        # THRESH        14
+        # ri_j        True
+        # ri_k        True
+        # aux_basis RIJK-def2-TZVP
+        # $end
+        ref_energy = -76.3811987938
+        ref_gradient = np.array([
+            [ 0.0084456, -0.0114582,  0.0030382, -0.0000256],
+            [ 0.0102732,  0.0003065, -0.0105467, -0.0000329],
+            [ 0.0000791, -0.0000347, -0.0000343, -0.0000100],
+        ]).T
+
+        ### The check threshold reflects the difference between Q-Chem and GPU4PySCF without ghost atom
+        assert np.abs(test_energy - ref_energy) < 1e-7
+        assert np.max(np.abs(test_gradient - ref_gradient)) < 3e-6
+
 if __name__ == "__main__":
     print("Full Tests for dft")
     unittest.main()

gpu4pyscf/grad/tests/test_rks_grad.py

@@ -252,8 +252,8 @@ def test_ghost_atom_grad_rks(self):
                 H     -1.430605   -0.826175   -0.900000
                 H      1.430605   -0.826175    0.900000
                 H      1.430605   -0.826175   -0.900000
-                ghost:Kr 0 0 1.5
-                ghost:Kr 0 0 -1.5
+                ghost:Ar 0 0 1.5
+                ghost:Ar 0 0 -1.5
             """, # Cyclopropane
             basis = "sto-3g",
             verbose = 0,
@@ -284,8 +284,8 @@ def test_ghost_atom_grad_rks(self):
         # H     -1.430605   -0.826175   -0.900000
         # H      1.430605   -0.826175    0.900000
         # H      1.430605   -0.826175   -0.900000
-        # @Kr 0 0 1.5
-        # @Kr 0 0 -1.5
+        # @Ar 0 0 1.5
+        # @Ar 0 0 -1.5
         # $end
 
         # $rem
@@ -304,11 +304,11 @@ def test_ghost_atom_grad_rks(self):
         # ri_k        True
         # aux_basis RIJK-def2-TZVP
         # $end
-        ref_energy = -116.3235770816
+        ref_energy = -116.3194912032
         ref_gradient = numpy.array([
-            [-0.0000000,  0.0620765, -0.0620765,  0.0000000,  0.0000000, -0.0442042, -0.0442042,  0.0442042,  0.0442042, -0.0000000,  0.0000000],
-            [ 0.0030145,  0.0424556,  0.0424556,  0.0068925,  0.0068925, -0.0257694, -0.0257694, -0.0257694, -0.0257694,  0.0006834,  0.0006834],
-            [-0.0000000,  0.0000000, -0.0000000, -0.0011744,  0.0011744,  0.0261449, -0.0261449,  0.0261449, -0.0261449,  0.0021925, -0.0021925],
+            [ 0.0000000,  0.0660107, -0.0660107,  0.0000000, -0.0000000, -0.0452116, -0.0452116,  0.0452116,  0.0452116, -0.0000000, -0.0000000],
+            [-0.0004307,  0.0448647,  0.0448647,  0.0074608,  0.0074608, -0.0263691, -0.0263691, -0.0263691, -0.0263691,  0.0006281,  0.0006281],
+            [ 0.0000000, -0.0000000,  0.0000000, -0.0008212,  0.0008212,  0.0267111, -0.0267111,  0.0267111, -0.0267111,  0.0023017, -0.0023017],
         ]).T
 
         assert numpy.abs(test_energy - ref_energy) < 2e-7

gpu4pyscf/grad/tests/test_vv10_grid.py

@@ -522,16 +522,16 @@ def test_grid_response_gradient_ghost_atom_far_away(self):
         test_gradient = gobj.kernel()
 
         # Reference come from removing the ghost atom far away
-        ref_energy = -76.02662600138008
+        ref_energy = -76.02665256617135
         ref_gradient = np.array([
-            [ 0.056476163601484,  0.074442899041095,  0.000000000000001],
-            [-0.063848118077074, -0.010594483907204, -0.000000000000001],
-            [ 0.00184695592196 , -0.063226209272444, -0.000000000000001],
+            [ 0.0557732315405932,  0.0737014727900909,  0.0000000000000011],
+            [-0.0637777885905613, -0.0102100391760461, -0.0000000000000007],
+            [ 0.0028399377889354, -0.0637134405199884, -0.0000000000000008],
             [ 0, 0, 0 ],
             [ 0, 0, 0 ],
-            [ 0.002209756185228,  0.000672985764181,  0.000000000000001],
-            [ 0.000001361947461,  0.000010775254639, -0.               ],
-            [ 0.003313880421019, -0.001305966880173, -0.               ],
+            [ 0.0014909325052633, -0.0000120529615038,  0.0000000000000009],
+            [ 0.0000273273841539,  0.0000089786373412, -0.0000000000000002],
+            [ 0.0036463593716917,  0.0002250812302098, -0.0000000000000003],
             [ 0, 0, 0 ],
             [ 0, 0, 0 ],
         ])