SMD solvent model (#62)

* save

* implement SMD solvent model

* updated readme

README.md

@@ -51,6 +51,7 @@ Features
 - Nonlocal functional correction (vv10) for SCF and gradient;
 - ECP is supported and calculated on CPU;
 - PCM solvent models, analytical gradients, and semi-analytical Hessian matrix;
+- SMD solvent models and solvation free energy
 
 Limitations
 --------

examples/16-smd.py

@@ -0,0 +1,38 @@
+# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+import numpy
+import pyscf
+from pyscf import gto, df
+from gpu4pyscf import scf, dft
+from gpu4pyscf.solvent import smd
+
+atom = '''
+O       0.0000000000    -0.0000000000     0.1174000000
+H      -0.7570000000    -0.0000000000    -0.4696000000
+H       0.7570000000     0.0000000000    -0.4696000000
+'''
+mol = pyscf.M(atom=atom, basis='def2-tzvpp', verbose=4)
+
+mf = dft.rks.RKS(mol, xc='HYB_GGA_XC_B3LYP')#.density_fit()
+mf = mf.SMD()
+mf.verbose = 4
+mf.grids.atom_grid = (99,590)
+mf.small_rho_cutoff = 1e-10
+mf.with_solvent.lebedev_order = 29 # 302 Lebedev grids
+mf.with_solvent.method = 'SMD'
+mf.with_solvent.solvent = 'water'
+e_tot = mf.kernel()
+print(e_tot)

gpu4pyscf/solvent/__init__.py

@@ -13,7 +13,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-from gpu4pyscf.solvent import pcm
+from gpu4pyscf.solvent import pcm, smd
 
 def PCM(method_or_mol, solvent_obj=None, dm=None):
     '''Initialize PCM model.
@@ -35,3 +35,24 @@ def PCM(method_or_mol, solvent_obj=None, dm=None):
         return pcm.pcm_for_scf(method_or_mol, solvent_obj, dm)
     else:
         raise NotImplementedError('PCM model only support SCF')
+
+def SMD(method_or_mol, solvent_obj=None, dm=None):
+    '''Initialize SMD model.
+
+    Examples:
+
+    >>> mf = PCM(scf.RHF(mol))
+    >>> mf.kernel()
+    >>> sol = PCM(mol)
+    >>> mc = PCM(CASCI(mf, 6, 6), sol)
+    >>> mc.kernel()
+    '''
+    from pyscf import gto
+    from pyscf import scf
+
+    if isinstance(method_or_mol, gto.mole.Mole):
+        return smd.SMD(method_or_mol)
+    elif isinstance(method_or_mol, scf.hf.SCF):
+        return pcm.pcm_for_scf(method_or_mol, solvent_obj, dm)
+    else:
+        raise NotImplementedError('SMD model only support SCF')

gpu4pyscf/solvent/_attach_solvent.py

@@ -87,10 +87,21 @@ def energy_elec(self, dm=None, h1e=None, vhf=None):
                 dm = self.make_rdm1()
             if getattr(vhf, 'e_solvent', None) is None:
                 vhf = self.get_veff(self.mol, dm)
+
             e_tot, e_coul = oldMF.energy_elec(self, dm, h1e, vhf)
             e_tot += vhf.e_solvent
             self.scf_summary['e_solvent'] = vhf.e_solvent.real
-            logger.debug(self, 'Solvent Energy = %.15g', vhf.e_solvent)
+
+            if self.with_solvent.method.upper() == 'SMD':
+                if self.with_solvent.e_cds is None:
+                    e_cds = self.with_solvent.get_cds()
+                    self.with_solvent.e_cds = e_cds
+                else:
+                    e_cds = self.with_solvent.e_cds
+                e_tot += e_cds
+                self.scf_summary['e_cds'] = e_cds
+                logger.info(self, f'CDS correction = {e_cds:.15f}')
+            logger.info(self, 'Solvent Energy = %.15g', vhf.e_solvent)
             return e_tot, e_coul
 
         def nuc_grad_method(self):

gpu4pyscf/solvent/grad/pcm.py

@@ -71,9 +71,9 @@ def get_dF_dA(surface):
         riJ = cupy.linalg.norm(ri_rJ, axis=-1)
         diJ = (riJ - R_in_J) / R_sw_J
         diJ[:,ia] = 1.0
-        diJ[diJ < 1e-12] = 0.0
+        diJ[diJ < 1e-8] = 0.0
         ri_rJ[:,ia,:] = 0.0
-        ri_rJ[diJ < 1e-12] = 0.0
+        ri_rJ[diJ < 1e-8] = 0.0
 
         fiJ = switch_h(diJ)
         dfiJ = grad_switch_h(diJ) / (fiJ * riJ * R_sw_J)

gpu4pyscf/solvent/pcm.py

@@ -40,7 +40,7 @@ def pcm_for_scf(mf, solvent_obj=None, dm=None):
 
 # Inject PCM to SCF, TODO: add it to other methods later
 from gpu4pyscf import scf
-scf.hf.RHF.PCM    = scf.hf.RHF.PCM    = pcm_for_scf
+scf.hf.RHF.PCM = pcm_for_scf
 
 # TABLE II,  J. Chem. Phys. 122, 194110 (2005)
 XI = {
@@ -75,8 +75,8 @@ def pcm_for_scf(mf, solvent_obj=None, dm=None):
     5810: 4.90792902522,
 }
 
-Bondi = radii.VDW
-Bondi[1] = 1.1/radii.BOHR      # modified version
+modified_Bondi = radii.VDW.copy()
+modified_Bondi[1] = 1.1/radii.BOHR      # modified version
 #radii_table = bondi * 1.2
 PI = numpy.pi
 
@@ -91,7 +91,7 @@ def switch_h(x):
     y[x>1] = 1.0
     return y
 
-def gen_surface(mol, ng=302, vdw_scale=1.2):
+def gen_surface(mol, ng=302, rad=modified_Bondi, vdw_scale=1.2, r_probe=0.0):
     '''J. Phys. Chem. A 1999, 103, 11060-11079'''
     unit_sphere = numpy.empty((ng,4))
     libdft.MakeAngularGrid(unit_sphere.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(ng))
@@ -100,7 +100,7 @@ def gen_surface(mol, ng=302, vdw_scale=1.2):
     atom_coords = cupy.asarray(mol.atom_coords(unit='B'))
     charges = mol.atom_charges()
     N_J = ng * cupy.ones(mol.natm)
-    R_J = cupy.asarray([vdw_scale*Bondi[chg] for chg in charges])
+    R_J = cupy.asarray([rad[chg] for chg in charges])
     R_sw_J = R_J * (14.0 / N_J)**0.5
     alpha_J = 1.0/2.0 + R_J/R_sw_J - ((R_J/R_sw_J)**2 - 1.0/28)**0.5
     R_in_J = R_J - alpha_J * R_sw_J
@@ -117,20 +117,20 @@ def gen_surface(mol, ng=302, vdw_scale=1.2):
     for ia in range(mol.natm):
         symb = mol.atom_symbol(ia)
         chg = gto.charge(symb)
-        r_vdw = vdw_scale*Bondi[chg]
+        r_vdw = rad[chg]
 
         atom_grid = r_vdw * unit_sphere[:,:3] + atom_coords[ia,:]
         #riJ = scipy.spatial.distance.cdist(atom_grid[:,:3], atom_coords)
         riJ = cupy.sum((atom_grid[:,None,:] - atom_coords[None,:,:])**2, axis=2)**0.5
         diJ = (riJ - R_in_J) / R_sw_J
         diJ[:,ia] = 1.0
-        diJ[diJ<1e-12] = 0.0
+        diJ[diJ<1e-8] = 0.0
 
         fiJ = switch_h(diJ)
 
         w = unit_sphere[:,3] * 4.0 * PI
         swf = cupy.prod(fiJ, axis=1)
-        idx = w*swf > 1e-16
+        idx = w*swf > 1e-12
 
         p0, p1 = p1, p1+sum(idx).get()
         gslice_by_atom.append([p0,p1])
@@ -208,10 +208,15 @@ def get_D_S(surface, with_S=True, with_D=False):
     return D, S
 
 class PCM(ddcosmo.DDCOSMO):
+    _keys = {
+        'method', 'vdw_scale', 'surface'
+    }
     def __init__(self, mol):
         ddcosmo.DDCOSMO.__init__(self, mol)
         self.method = 'C-PCM'
         self.vdw_scale = 1.2 # default value in qchem
+        self.r_probe = 0.0
+        self.radii_table = None
         self.surface = {}
         self._intermediates = {}
 
@@ -231,13 +236,14 @@ def dump_flags(self, verbose=None):
         return self
 
     def build(self, ng=None):
-        vdw_scale = self.vdw_scale
-        self.radii_table = vdw_scale * Bondi
+        if self.radii_table is None:
+            vdw_scale = self.vdw_scale
+            self.radii_table = vdw_scale * modified_Bondi + self.r_probe
         mol = self.mol
         if ng is None:
             ng = gen_grid.LEBEDEV_ORDER[self.lebedev_order]
 
-        self.surface = gen_surface(mol, ng=ng, vdw_scale=vdw_scale)
+        self.surface = gen_surface(mol, rad=self.radii_table, ng=ng)
         self._intermediates = {}
         F, A = get_F_A(self.surface)
         D, S = get_D_S(self.surface, with_S=True, with_D=True)
@@ -251,7 +257,7 @@ def build(self, ng=None):
             f_epsilon = (epsilon - 1.0)/(epsilon + 1.0/2.0)
             K = S
             R = -f_epsilon * cupy.eye(K.shape[0])
-        elif self.method.upper() in ['IEF-PCM', 'IEFPCM']:
+        elif self.method.upper() in ['IEF-PCM', 'IEFPCM', 'SMD']:
             f_epsilon = (epsilon - 1.0)/(epsilon + 1.0)
             DA = D*A
             DAS = cupy.dot(DA, S)