From 91701729f9e3cac657824f4b058248b7e9199274 Mon Sep 17 00:00:00 2001
From: Henry Wang <henry.wang1@bytedance.com>
Date: Tue, 7 Jan 2025 08:43:36 +0800
Subject: [PATCH] Remove the numerical implementation of dVia/dx

---
 gpu4pyscf/solvent/hessian/pcm.py | 41 --------------------------------
 gpu4pyscf/solvent/hessian/smd.py |  6 ++---
 2 files changed, 3 insertions(+), 44 deletions(-)

diff --git a/gpu4pyscf/solvent/hessian/pcm.py b/gpu4pyscf/solvent/hessian/pcm.py
index 75e8ecc1..74a51932 100644
--- a/gpu4pyscf/solvent/hessian/pcm.py
+++ b/gpu4pyscf/solvent/hessian/pcm.py
@@ -153,47 +153,6 @@ def pcm_grad_scanner(mol):
     pcmobj.reset(pmol)
     return de
 
-def fd_grad_vmat(pcmobj, dm, mo_coeff, mo_occ, atmlst=None, verbose=None):
-    '''
-    dv_solv / da
-    slow version with finite difference
-    '''
-    log = logger.new_logger(pcmobj, verbose)
-    t1 = log.init_timer()
-    pmol = pcmobj.mol.copy()
-    mol = pmol.copy()
-    if atmlst is None:
-        atmlst = range(mol.natm)
-    nao, nmo = mo_coeff.shape
-    mocc = mo_coeff[:,mo_occ>0]
-    nocc = mocc.shape[1]
-    coords = mol.atom_coords(unit='Bohr')
-    def pcm_vmat_scanner(mol):
-        pcmobj.reset(mol)
-        e, v = pcmobj._get_vind(dm)
-        return v
-
-    mol.verbose = 0
-    vmat = cupy.empty([len(atmlst), 3, nao, nocc])
-    eps = 1e-3
-    for i0, ia in enumerate(atmlst):
-        for ix in range(3):
-            dv = numpy.zeros_like(coords)
-            dv[ia,ix] = eps
-            mol.set_geom_(coords + dv, unit='Bohr')
-            vmat0 = pcm_vmat_scanner(mol)
-
-            mol.set_geom_(coords - dv, unit='Bohr')
-            vmat1 = pcm_vmat_scanner(mol)
-
-            grad_vmat = (vmat0 - vmat1)/2.0/eps
-            grad_vmat = contract("ij,jq->iq", grad_vmat, mocc)
-            grad_vmat = contract("iq,ip->pq", grad_vmat, mo_coeff)
-            vmat[i0,ix] = grad_vmat
-    t1 = log.timer_debug1('computing solvent grad veff', *t1)
-    pcmobj.reset(pmol)
-    return vmat
-
 def analytic_grad_vmat(pcmobj, dm, mo_coeff, mo_occ, atmlst=None, verbose=None):
     '''
     dv_solv / da
diff --git a/gpu4pyscf/solvent/hessian/smd.py b/gpu4pyscf/solvent/hessian/smd.py
index 644c0e3f..49897d74 100644
--- a/gpu4pyscf/solvent/hessian/smd.py
+++ b/gpu4pyscf/solvent/hessian/smd.py
@@ -154,7 +154,7 @@ def make_h1(self, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None):
         h1ao = super().make_h1(mo_coeff, mo_occ, atmlst=atmlst, verbose=verbose)
         if isinstance(self.base, scf.hf.RHF):
             dm = self.base.make_rdm1(ao_repr=True)
-            dv = pcm_hess.fd_grad_vmat(self.base.with_solvent, dm, mo_coeff, mo_occ, atmlst=atmlst, verbose=verbose)
+            dv = pcm_hess.analytic_grad_vmat(self.base.with_solvent, dm, mo_coeff, mo_occ, atmlst=atmlst, verbose=verbose)
             for i0, ia in enumerate(atmlst):
                 h1ao[i0] += dv[i0]
             return h1ao
@@ -163,8 +163,8 @@ def make_h1(self, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None):
             solvent = self.base.with_solvent
             dm = self.base.make_rdm1(ao_repr=True)
             dm = dm[0] + dm[1]
-            dva = pcm_hess.fd_grad_vmat(solvent, dm, mo_coeff[0], mo_occ[0], atmlst=atmlst, verbose=verbose)
-            dvb = pcm_hess.fd_grad_vmat(solvent, dm, mo_coeff[1], mo_occ[1], atmlst=atmlst, verbose=verbose)
+            dva = pcm_hess.analytic_grad_vmat(solvent, dm, mo_coeff[0], mo_occ[0], atmlst=atmlst, verbose=verbose)
+            dvb = pcm_hess.analytic_grad_vmat(solvent, dm, mo_coeff[1], mo_occ[1], atmlst=atmlst, verbose=verbose)
             for i0, ia in enumerate(atmlst):
                 h1aoa[i0] += dva[i0]
                 h1aob[i0] += dvb[i0]