Fix tests

gpu4pyscf/solvent/hessian/pcm.py

@@ -365,7 +365,7 @@ def get_v_dot_d2D_dot_q(d2D, v_left, q_right, natom, gridslice):
 def get_v_dot_d2DT_dot_q(d2D, v_left, q_right, natom, gridslice):
     return get_v_dot_d2D_dot_q(d2D.transpose(0,1,3,2), v_left, q_right, natom, gridslice)
 
-def analytical_hess_solver(pcmobj, dm):
+def analytical_hess_solver(pcmobj, dm, verbose=None):
     if not pcmobj._intermediates:
         pcmobj.build()
     dm_cache = pcmobj._intermediates.get('dm', None)
@@ -374,7 +374,7 @@ def analytical_hess_solver(pcmobj, dm):
     else:
         pcmobj._get_vind(dm)
     mol = pcmobj.mol
-    log = logger.new_logger(mol, mol.verbose)
+    log = logger.new_logger(mol, verbose)
     t1 = log.init_timer()
 
     natom = mol.natm

gpu4pyscf/solvent/hessian/smd.py

@@ -22,8 +22,6 @@
 from gpu4pyscf import scf
 from gpu4pyscf.lib import logger
 from gpu4pyscf.solvent import smd
-from gpu4pyscf.solvent.grad import smd as smd_grad
-from gpu4pyscf.solvent.grad import pcm as pcm_grad
 from gpu4pyscf.solvent.hessian import pcm as pcm_hess
 from gpu4pyscf.hessian.jk import _ao2mo
 
@@ -60,45 +58,6 @@ def smd_grad_scanner(mol):
     t1 = log.timer_debug1('solvent energy', *t1)
     return hess_cds # hartree
 
-
-def hess_elec(smdobj, dm, verbose=None):
-    '''
-    slow version with finite difference
-    TODO: use analytical hess_nuc
-    '''
-    log = logger.new_logger(smdobj, verbose)
-    t1 = log.init_timer()
-    pmol = smdobj.mol.copy()
-    mol = pmol.copy()
-    coords = mol.atom_coords(unit='Bohr')
-
-    def pcm_grad_scanner(mol):
-        # TODO: use more analytical forms
-        smdobj.reset(mol)
-        e, v = smdobj._get_vind(dm)
-        #return grad_elec(smdobj, dm)
-        grad = pcm_grad.grad_nuc(smdobj, dm)
-        grad+= smd_grad.grad_solver(smdobj, dm)
-        grad+= pcm_grad.grad_qv(smdobj, dm)
-        return grad
-
-    mol.verbose = 0
-    de = np.zeros([mol.natm, mol.natm, 3, 3])
-    eps = 1e-3
-    for ia in range(mol.natm):
-        for ix in range(3):
-            dv = np.zeros_like(coords)
-            dv[ia,ix] = eps
-            mol.set_geom_(coords + dv, unit='Bohr')
-            g0 = pcm_grad_scanner(mol)
-
-            mol.set_geom_(coords - dv, unit='Bohr')
-            g1 = pcm_grad_scanner(mol)
-            de[ia,:,ix] = (g0 - g1)/2.0/eps
-    t1 = log.timer_debug1('solvent energy', *t1)
-    smdobj.reset(pmol)
-    return de
-
 def make_hess_object(hess_method):
     '''For hess_method in vacuum, add nuclear Hessian of solvent smdobj'''
     if hess_method.base.with_solvent.frozen:
@@ -140,8 +99,9 @@ def kernel(self, *args, dm=None, atmlst=None, **kwargs):
             dm = dm[0] + dm[1]
         is_equilibrium = self.base.with_solvent.equilibrium_solvation
         self.base.with_solvent.equilibrium_solvation = True
-        self.de_solvent = pcm_hess.hess_elec(self.base.with_solvent, dm, verbose=self.verbose)
-        #self.de_solvent+= hess_nuc(self.base.with_solvent)
+        self.de_solvent  =    pcm_hess.analytical_hess_nuc(self.base.with_solvent, dm, verbose=self.verbose)
+        self.de_solvent +=     pcm_hess.analytical_hess_qv(self.base.with_solvent, dm, verbose=self.verbose)
+        self.de_solvent += pcm_hess.analytical_hess_solver(self.base.with_solvent, dm, verbose=self.verbose)
         self.de_solute = super().kernel(*args, **kwargs)
         self.de_cds = get_cds(self.base.with_solvent)
         self.de = self.de_solute + self.de_solvent + self.de_cds

gpu4pyscf/solvent/tests/test_pcm_grad.py

@@ -36,6 +36,7 @@ def setUpModule():
     mol.basis = 'sto3g'
     mol.output = '/dev/null'
     mol.build(verbose=0)
+    # Warning: This system has all orbitals filled, which is FAR from physical
     mol.nelectron = mol.nao * 2
     epsilon = 35.9
     lebedev_order = 3
@@ -169,11 +170,14 @@ def test_grad_IEFPCM(self):
 
     def test_grad_SSVPE(self):
         grad = _grad_with_solvent('SS(V)PE')
-        g0 = numpy.asarray(
-            [[ 3.42479745e-15, -1.00280742e-16, -1.61117735e+00],
-            [ 1.07135985e+00, -6.97375148e-16,  8.05588676e-01],
-            [-1.07135985e+00,  7.91425487e-16,  8.05588676e-01]]
-        )
+        # Note: This reference value is obtained via finite difference with dx = 1e-5
+        #       QChem 6.1 has a bug in SSVPE gradient, they use the IEFPCM gradient algorithm
+        #       to compute SSVPE gradient, which is wrong.
+        g0 = numpy.asarray([
+            [ 2.13162821e-09,  0.00000000e+00, -5.89116532e-02],
+            [ 2.28518083e-02, -2.13162821e-09,  2.94558255e-02],
+            [-2.28518090e-02,  0.00000000e+00,  2.94558234e-02],
+        ])
         print(f"Gradient error in RHF with SS(V)PE: {numpy.linalg.norm(g0 - grad)}")
         assert numpy.linalg.norm(g0 - grad) < 1e-6