native support dispersion correction (pyscf#1916)

* 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

* native support dispersion correction

* fixed a bug in pcm

* Update hf.py

* move dispersion to addons

* remove disp in RKS

* call get_dispersion in kernels

* added unit test for d4

* fixed dispersion correction in testing h2o

* updated ci

* skip unittest if dftd3 or dftd4 is missing

.github/workflows/ci_linux/python_deps.sh

@@ -7,6 +7,8 @@ version=$(python -c 'import sys; version=sys.version_info[:2]; print("{0}.{1}".f
 if [ $version != '3.12' ]; then
     pip install geometric
     pip install spglib
+    pip install dftd3
+    pip install dftd4
 fi
 
 #cppe

pyscf/dft/rks.py

@@ -243,6 +243,7 @@ def energy_elec(ks, dm=None, h1e=None, vhf=None):
     ecoul = vhf.ecoul.real
     exc = vhf.exc.real
     e2 = ecoul + exc
+
     ks.scf_summary['e1'] = e1
     ks.scf_summary['coul'] = ecoul
     ks.scf_summary['exc'] = exc

pyscf/dft/test/test_h2o.py

@@ -19,6 +19,18 @@
 from pyscf import lib
 from pyscf import dft
 
+
+import sys
+try:
+    import dftd3
+except ImportError:
+    pass
+
+try:
+    import dftd4
+except ImportError:
+    pass
+
 def setUpModule():
     global h2o, h2osym, h2o_cation, h2osym_cation
     h2o = gto.Mole()
@@ -505,6 +517,14 @@ def test_camb3lyp_rsh_omega(self):
         mf2.kernel()
         self.assertAlmostEqual(mf1.e_tot, -76.36649222362115, 9)
 
+    @unittest.skipIf('dftd3' not in sys.modules, "requires the dftd3 library")
+    def test_dispersion(self):
+        mf = dft.RKS(h2o)
+        mf.xc = 'B3LYP'
+        mf.disp = 'd3bj'
+        mf.run(xc='B3LYP')
+        self.assertAlmostEqual(mf.e_tot, -76.38945547396322, 9)
+
     def test_reset(self):
         mf = dft.RKS(h2o).newton()
         mf.reset(h2osym)

pyscf/grad/__init__.py

@@ -32,6 +32,7 @@
 from . import dhf
 from . import uhf
 from . import rohf
+from . import dispersion
 from .rhf import Gradients as RHF
 from .dhf import Gradients as DHF
 from .uhf import Gradients as UHF

pyscf/grad/dispersion.py

@@ -0,0 +1,59 @@
+#!/usr/bin/env python
+# Copyright 2014-2023 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# Author: Xiaojie Wu <[email protected]>
+#
+
+'''
+gradient of dispersion correction for HF and DFT
+'''
+
+import numpy
+from pyscf.scf.hf import KohnShamDFT
+
+def get_dispersion(mf_grad, disp_version=None):
+    '''gradient of dispersion correction for RHF/RKS'''
+    if disp_version is None:
+        disp_version = mf_grad.base.disp
+    mol = mf_grad.base.mol
+    disp_version = mf_grad.base.disp
+    if disp_version is None:
+        return numpy.zeros([mol.natm,3])
+
+    if isinstance(mf_grad.base, KohnShamDFT):
+        method = mf_grad.base.xc
+    else:
+        method = 'hf'
+
+    if disp_version[:2].upper() == 'D3':
+        # raised error in SCF module, assuming dftd3 installed
+        import dftd3.pyscf as disp
+        d3 = disp.DFTD3Dispersion(mol, xc=method, version=disp_version)
+        _, g_d3 = d3.kernel()
+        return g_d3
+    elif disp_version[:2].upper() == 'D4':
+        from pyscf.data.elements import charge
+        atoms = numpy.array([ charge(a[0]) for a in mol._atom])
+        coords = mol.atom_coords()
+        from dftd4.interface import DampingParam, DispersionModel
+        model = DispersionModel(atoms, coords)
+        res = model.get_dispersion(DampingParam(method=method), grad=True)
+        return res.get("gradient")
+    else:
+        raise RuntimeError(f'dispersion correction: {disp_version} is not supported.')
+
+# Inject to Gradient
+from pyscf import grad
+grad.rhf.GradientsBase.get_dispersion = get_dispersion

pyscf/grad/rhf.py

@@ -416,6 +416,8 @@ def kernel(self, mo_energy=None, mo_coeff=None, mo_occ=None, atmlst=None):
         self.de = de + self.grad_nuc(atmlst=atmlst)
         if self.mol.symmetry:
             self.de = self.symmetrize(self.de, atmlst)
+        if self.base.disp is not None:
+            self.de += self.get_dispersion()
         logger.timer(self, 'SCF gradients', *cput0)
         self._finalize()
         return self.de

pyscf/grad/test/test_rhf.py

@@ -18,6 +18,17 @@
 from pyscf import gto, scf, lib
 from pyscf import grad
 
+import sys
+try:
+    import dftd3
+except ImportError:
+    pass
+
+try:
+    import dftd4
+except ImportError:
+    pass
+
 def setUpModule():
     global mol
     mol = gto.Mole()
@@ -71,6 +82,30 @@ def test_df_rhf_grad(self):
         e2 = mfs('O  0.  0.  0.001; H  0.  -0.757  0.587; H  0.  0.757   0.587')
         self.assertAlmostEqual(g[0,2], (e2-e1)/0.002*lib.param.BOHR, 5)
 
+    @unittest.skipIf('dftd3' not in sys.modules, "requires the dftd3 library")
+    def test_rhf_d3_grad(self):
+        mf = scf.RHF(mol)
+        mf.disp = 'd3bj'
+        g_scan = mf.nuc_grad_method().as_scanner()
+        g = g_scan(mol)[1]
+
+        mf_scan = mf.as_scanner()
+        e1 = mf_scan('O  0.  0. -0.001; H  0.  -0.757  0.587; H  0.  0.757   0.587')
+        e2 = mf_scan('O  0.  0.  0.001; H  0.  -0.757  0.587; H  0.  0.757   0.587')
+        self.assertAlmostEqual((e2-e1)/0.002*lib.param.BOHR, g[0,2], 5)
+
+    @unittest.skipIf('dftd4' not in sys.modules, "requires the dftd4 library")
+    def test_rhf_d4_grad(self):
+        mf = scf.RHF(mol)
+        mf.disp = 'd4'
+        g_scan = mf.nuc_grad_method().as_scanner()
+        g = g_scan(mol)[1]
+
+        mf_scan = mf.as_scanner()
+        e1 = mf_scan('O  0.  0. -0.001; H  0.  -0.757  0.587; H  0.  0.757   0.587')
+        e2 = mf_scan('O  0.  0.  0.001; H  0.  -0.757  0.587; H  0.  0.757   0.587')
+        self.assertAlmostEqual((e2-e1)/0.002*lib.param.BOHR, g[0,2], 5)
+
     def test_x2c_rhf_grad(self):
         h2o = gto.Mole()
         h2o.verbose = 0

pyscf/hessian/__init__.py

@@ -18,6 +18,7 @@
 
 from pyscf.hessian import rhf
 from pyscf.hessian import uhf
+from pyscf.hessian import dispersion
 from pyscf.hessian.rhf import Hessian as RHF
 from pyscf.hessian.uhf import Hessian as UHF
 from pyscf.hessian.rhf import hess_nuc

pyscf/hessian/dispersion.py

@@ -0,0 +1,94 @@
+#!/usr/bin/env python
+# Copyright 2014-2023 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# Author: Xiaojie Wu <[email protected]>
+#
+
+'''
+Hessian of dispersion correction for HF and DFT
+'''
+
+
+import numpy
+from pyscf.scf.hf import KohnShamDFT
+
+def get_dispersion(hessobj, disp_version=None):
+    if disp_version is None:
+        disp_version = hessobj.base.disp
+    mol = hessobj.base.mol
+    natm = mol.natm
+    mf = hessobj.base
+    h_disp = numpy.zeros([natm,natm,3,3])
+    if disp_version is None:
+        return h_disp
+    if isinstance(hessobj.base, KohnShamDFT):
+        method = hessobj.base.xc
+    else:
+        method = 'hf'
+
+    if mf.disp[:2].upper() == 'D3':
+        import dftd3.pyscf as disp
+        coords = hessobj.mol.atom_coords()
+        mol = mol.copy()
+        eps = 1e-5
+        for i in range(natm):
+            for j in range(3):
+                coords[i,j] += eps
+                mol.set_geom_(coords, unit='Bohr')
+                d3 = disp.DFTD3Dispersion(mol, xc=method, version=mf.disp)
+                _, g1 = d3.kernel()
+
+                coords[i,j] -= 2.0*eps
+                mol.set_geom_(coords, unit='Bohr')
+                d3 = disp.DFTD3Dispersion(mol, xc=method, version=mf.disp)
+                _, g2 = d3.kernel()
+
+                coords[i,j] += eps
+                h_disp[i,:,j,:] = (g1 - g2)/(2.0*eps)
+            return h_disp
+
+    elif mf.disp[:2].upper() == 'D4':
+        from pyscf.data.elements import charge
+        atoms = numpy.array([ charge(a[0]) for a in mol._atom])
+        coords = mol.atom_coords()
+        natm = mol.natm
+        from dftd4.interface import DampingParam, DispersionModel
+        params = DampingParam(method=method)
+        mol = mol.copy()
+        eps = 1e-5
+        for i in range(natm):
+            for j in range(3):
+                coords[i,j] += eps
+                mol.set_geom_(coords, unit='Bohr')
+                model = DispersionModel(atoms, coords)
+                res = model.get_dispersion(params, grad=True)
+                g1 = res.get("gradient")
+
+                coords[i,j] -= 2.0*eps
+                mol.set_geom_(coords, unit='Bohr')
+                model = DispersionModel(atoms, coords)
+                res = model.get_dispersion(params, grad=True)
+                g2 = res.get("gradient")
+
+                coords[i,j] += eps
+                h_disp[i,:,j,:] = (g1 - g2)/(2.0*eps)
+
+        return h_disp
+    else:
+        raise RuntimeError(f'dispersion correction: {disp_version} is not supported.')
+
+# Inject to SCF class
+from pyscf import hessian
+hessian.rhf.HessianBase.get_dispersion = get_dispersion

pyscf/hessian/rhf.py

@@ -27,7 +27,7 @@
 from pyscf import lib
 from pyscf import gto
 from pyscf.lib import logger
-from pyscf.scf import _vhf
+from pyscf.scf import _vhf, hf
 from pyscf.scf import cphf
 
 # import _response_functions to load gen_response methods in SCF class
@@ -583,6 +583,8 @@ def kernel(self, mo_energy=None, mo_coeff=None, mo_occ=None, atmlst=None):
 
         de = self.hess_elec(mo_energy, mo_coeff, mo_occ, atmlst=atmlst)
         self.de = de + self.hess_nuc(self.mol, atmlst=atmlst)
+        if self.base.disp is not None:
+            self.de += self.get_dispersion()
         return self.de
     hess = kernel
 

pyscf/hessian/test/test_rks.py

@@ -18,6 +18,17 @@
 from pyscf import gto, dft, lib
 from pyscf import grad, hessian
 
+import sys
+try:
+    import dftd3
+except ImportError:
+    pass
+
+try:
+    import dftd4
+except ImportError:
+    pass
+
 def setUpModule():
     global mol, h4
     mol = gto.Mole()
@@ -116,6 +127,40 @@ def test_finite_diff_b3lyp_hess(self):
         #FIXME: errors seems too big
         self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 3)
 
+    @unittest.skipIf('dftd3' not in sys.modules, "requires the dftd3 library")
+    def test_finite_diff_b3lyp_d3_hess(self):
+        mf = dft.RKS(mol)
+        mf.conv_tol = 1e-14
+        mf.xc = 'b3lyp'
+        mf.disp = 'd3bj'
+        e0 = mf.kernel()
+        hess = mf.Hessian().kernel()
+        self.assertAlmostEqual(lib.fp(hess), -0.7586078053657133, 6)
+
+        g_scanner = mf.nuc_grad_method().as_scanner()
+        pmol = mol.copy()
+        e1 = g_scanner(pmol.set_geom_('O  0. 0. 0.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        e2 = g_scanner(pmol.set_geom_('O  0. 0. -.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        #FIXME: errors seems too big
+        self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 3)
+
+    @unittest.skipIf('dftd4' not in sys.modules, "requires the dftd4 library")
+    def test_finite_diff_b3lyp_d4_hess(self):
+        mf = dft.RKS(mol)
+        mf.conv_tol = 1e-14
+        mf.xc = 'b3lyp'
+        mf.disp = 'd4'
+        e0 = mf.kernel()
+        hess = mf.Hessian().kernel()
+        self.assertAlmostEqual(lib.fp(hess), -0.7588415571313422, 6)
+
+        g_scanner = mf.nuc_grad_method().as_scanner()
+        pmol = mol.copy()
+        e1 = g_scanner(pmol.set_geom_('O  0. 0. 0.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        e2 = g_scanner(pmol.set_geom_('O  0. 0. -.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        #FIXME: errors seems too big
+        self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 3)
+
     def test_finite_diff_wb97x_hess(self):
         mf = dft.RKS(mol)
         mf.conv_tol = 1e-14

pyscf/scf/__init__.py

@@ -92,7 +92,7 @@
         SCF converged or not
     e_tot : float
         Total HF energy (electronic energy plus nuclear repulsion)
-    mo_energy : 
+    mo_energy :
         Orbital energies
     mo_occ
         Orbital occupancy
@@ -115,6 +115,7 @@
 from pyscf.scf import chkfile
 from pyscf.scf import addons
 from pyscf.scf import diis
+from pyscf.scf import dispersion
 from pyscf.scf.diis import DIIS, CDIIS, EDIIS, ADIIS
 from pyscf.scf.uhf import spin_square
 from pyscf.scf.hf import get_init_guess