add pcm solvent models (pyscf#1884)

* added solvent models

* add example for RHF

* cleanup variables

* support casci casscf and ccsd

* uncomment unittests

* change example name

* update reset

* for flake8

* fixed a bug in soscf/newton_ah.py

* updated for recent master changes

* remove whitespace

* remove whitespace

examples/solvent/05-pcm.py

@@ -0,0 +1,49 @@
+#!/usr/bin/env python
+'''
+An example of using PCM solvent models in the mean-field calculations.
+'''
+
+from pyscf import gto, scf, dft, cc, solvent, mcscf
+from pyscf.solvent import pcm
+
+mol = gto.M(atom='''
+C        0.000000    0.000000             -0.542500
+O        0.000000    0.000000              0.677500
+H        0.000000    0.9353074360871938   -1.082500
+H        0.000000   -0.9353074360871938   -1.082500
+            ''',
+            verbose = 4)
+
+cm = pcm.PCM(mol)
+cm.eps = 32.613  # methanol dielectric constant
+cm.method = 'C-PCM' # or COSMO, IEF-PCM, SS(V)PE, see https://manual.q-chem.com/5.4/topic_pcm-em.html
+cm.lebedev_order = 29 # lebedev grids on the cavity surface, lebedev_order=29  <--> # of grids = 302
+
+# Hartree-Fock
+mf = scf.RHF(mol).PCM(cm)
+mf.kernel()
+
+g = mf.nuc_grad_method()
+grad = g.kernel()
+
+# DFT
+mf = dft.RKS(mol, xc='b3lyp').PCM(cm)
+mf.kernel()
+
+# calculate gradient of PCM models
+g = mf.nuc_grad_method()
+grad = g.kernel()
+
+# CCSD, solvent for MP2 can be calcualted in the similar way
+mf = scf.RHF(mol) # default method is C-PCM
+mf.kernel()
+mycc = cc.CCSD(mf).PCM()
+mycc.kernel()
+de = mycc.nuc_grad_method().as_scanner()(mol)
+
+# CASCI
+mf = scf.RHF(mol).PCM()
+mf.kernel()
+mc = solvent.PCM(mcscf.CASCI(mf,2,2))
+de = mc.nuc_grad_method().as_scanner()(mol)
+

pyscf/solvent/__init__.py

@@ -111,3 +111,35 @@ def PE(method_or_mol, solvent_obj, dm=None):
         return pol_embed.pe_for_tdscf(method_or_mol, solvent_obj, dm)
     else:
         return pol_embed.pe_for_post_scf(method_or_mol, solvent_obj, dm)
+
+def PCM(method_or_mol, solvent_obj=None, dm=None):
+    '''Initialize PCM 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, mcscf
+    from pyscf import tdscf
+    from pyscf.solvent import pcm
+
+    if isinstance(method_or_mol, gto.mole.Mole):
+        return pcm.PCM(method_or_mol)
+
+    elif isinstance(method_or_mol, scf.hf.SCF):
+        return pcm.pcm_for_scf(method_or_mol, solvent_obj, dm)
+    elif isinstance(method_or_mol, mcscf.mc1step.CASSCF):
+        return pcm.pcm_for_casscf(method_or_mol, solvent_obj, dm)
+    elif isinstance(method_or_mol, mcscf.casci.CASCI):
+        return pcm.pcm_for_casci(method_or_mol, solvent_obj, dm)
+    elif isinstance(method_or_mol, (tdscf.rhf.TDA, tdscf.rhf.TDHF)):
+        return pcm.pcm_for_tdscf(method_or_mol, solvent_obj, dm)
+    else:
+        return pcm.pcm_for_post_scf(method_or_mol, solvent_obj, dm)
+
+PCM = PCM

pyscf/solvent/ddcosmo.py

@@ -297,7 +297,6 @@ def ddcosmo_for_tdscf(method, solvent_obj=None, dm=None):
 mcscf.casci.CASCI.ddCOSMO = mcscf.casci.CASCI.DDCOSMO = ddcosmo_for_casci
 mcscf.mc1step.CASSCF.ddCOSMO = mcscf.mc1step.CASSCF.DDCOSMO = ddcosmo_for_casscf
 
-
 # Keep gen_ddcosmo_solver for backward compatibility
 def gen_ddcosmo_solver(pcmobj, verbose=None):
     '''Generate ddcosmo function to compute energy and potential matrix
@@ -811,6 +810,7 @@ def _get_vind(self, dm):
             f_epsilon = 1
         epcm = .5 * f_epsilon * numpy.einsum('jx,jx', psi, Xvec)
         vpcm = .5 * f_epsilon * vmat
+
         return epcm, vpcm
 
     def _B_dot_x(self, dm):