Improve the RSH treatment in DFT (pyscf#2525)

* Compute LR or SR K matrices directly for certain XC functionals

* Fix libxc.rsh_coeff parser for custom RSH such as "hse06+HF"

* Better SR/LR treatments for RSH in TDDFT and scf/_response_fuctions.py

* Remove TDDFT code that were never used

* Fix various issues for previous commits and update tests

* Fix tests

* Drop invalid test

* Fix pbc.tdscf tests and update docstrings

* Improve TDDFT eigen solver

* Update _lr_eig

* error introduced in previous commit

* Bugfix and optimization for lr_eig

* Restore compatibility between lr_eig and real_eig

* Fix tdscf symmetric_orth function for complex vectors

* Reduce numerical errors in TDHF diagonalization

* Dimension check

examples/pbc/11-gamma_point_all_electron_scf.py

@@ -3,10 +3,6 @@
 '''
 Gamma point Hartree-Fock/DFT for all-electron calculation
 
-The default FFT-based 2-electron integrals may not be accurate enough for
-all-electron calculation.  It's recommended to use MDF (mixed density fitting)
-technique to improve the accuracy.
-
 See also
 examples/df/00-with_df.py
 examples/df/01-auxbasis.py
@@ -33,11 +29,6 @@
 mf = scf.RHF(cell).density_fit()
 mf.kernel()
 
-# Mixed density fitting is another option for all-electron calculations
-mf = scf.RHF(cell).mix_density_fit()
-mf.with_df.mesh = [10]*3  # Tune #PWs in MDF for performance/accuracy balance
-mf.kernel()
-
 # Or use even-tempered Gaussian basis as auxiliary fitting functions.
 # The following auxbasis is generated based on the expression
 #    alpha = a * 1.7^i   i = 0..N

examples/pbc/12-gamma_point_post_hf.py

@@ -23,9 +23,7 @@
     verbose = 4,
 )
 
-mf = scf.RHF(cell).density_fit()
-mf.with_df.mesh = [10]*3
-mf.kernel()
+mf = scf.RHF(cell).density_fit().run()
 
 #
 # Import CC, TDDFT module from the molecular implementations

pyscf/dft/gks.py

@@ -89,36 +89,40 @@ def get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1):
             logger.debug(ks, 'nelec with nlc grids = %s', n)
         t0 = logger.timer(ks, 'vxc', *t0)
 
+    incremental_jk = (ks._eri is None and ks.direct_scf and
+                      getattr(vhf_last, 'vj', None) is not None)
+    if incremental_jk:
+        _dm = numpy.asarray(dm) - numpy.asarray(dm_last)
+    else:
+        _dm = dm
     if not ni.libxc.is_hybrid_xc(ks.xc):
         vk = None
-        if (ks._eri is None and ks.direct_scf and
-            getattr(vhf_last, 'vj', None) is not None):
-            ddm = numpy.asarray(dm) - numpy.asarray(dm_last)
-            vj = ks.get_j(mol, ddm, hermi)
+        vj = ks.get_j(mol, _dm, hermi)
+        if incremental_jk:
             vj += vhf_last.vj
-        else:
-            vj = ks.get_j(mol, dm, hermi)
         vxc += vj
     else:
         omega, alpha, hyb = ni.rsh_and_hybrid_coeff(ks.xc, spin=mol.spin)
-        if (ks._eri is None and ks.direct_scf and
-            getattr(vhf_last, 'vk', None) is not None):
-            ddm = numpy.asarray(dm) - numpy.asarray(dm_last)
-            vj, vk = ks.get_jk(mol, ddm, hermi)
+        if omega == 0:
+            vj, vk = ks.get_jk(mol, _dm, hermi)
+            vk *= hyb
+        elif alpha == 0: # LR=0, only SR exchange
+            vj = ks.get_j(mol, _dm, hermi)
+            vk = ks.get_k(mol, _dm, hermi, omega=-omega)
             vk *= hyb
-            if omega != 0:
-                vklr = ks.get_k(mol, ddm, hermi, omega=omega)
-                vklr *= (alpha - hyb)
-                vk += vklr
+        elif hyb == 0: # SR=0, only LR exchange
+            vj = ks.get_j(mol, _dm, hermi)
+            vk = ks.get_k(mol, _dm, hermi, omega=omega)
+            vk *= alpha
+        else: # SR and LR exchange with different ratios
+            vj, vk = ks.get_jk(mol, _dm, hermi)
+            vk *= hyb
+            vklr = ks.get_k(mol, _dm, hermi, omega=omega)
+            vklr *= (alpha - hyb)
+            vk += vklr
+        if incremental_jk:
             vj += vhf_last.vj
             vk += vhf_last.vk
-        else:
-            vj, vk = ks.get_jk(mol, dm, hermi)
-            vk *= hyb
-            if omega != 0:
-                vklr = ks.get_k(mol, dm, hermi, omega=omega)
-                vklr *= (alpha - hyb)
-                vk += vklr
         vxc += vj - vk
 
         if ground_state:

pyscf/dft/libxc.py

@@ -410,9 +410,14 @@ def rsh_coeff(xc_code):
         if 'SR_HF' in xc_code or 'LR_HF' in xc_code or 'RSH(' in xc_code:
             check_omega = False
 
-    hyb, fn_facs = parse_xc(xc_code)
-    hyb, alpha, omega = hyb
-    beta = hyb - alpha
+    (hyb, alpha, omega), fn_facs = parse_xc(xc_code)
+    if omega == 0:
+        # SR and LR Coulomb share the same coefficients
+        # Note: this change breaks compatibility with pyscf-2.7
+        assert hyb == alpha
+        beta = 0.
+    else:
+        beta = hyb - alpha
     rsh_pars = [omega, alpha, beta]
     rsh_tmp = (ctypes.c_double*3)()
     for xid, fac in fn_facs:
@@ -639,8 +644,6 @@ def possible_c_for(key):
             # dftd3 cannot be used in a custom xc description
             assert '-d3' not in token
             parse_token(token, 'compound XC', search_xc_alias=True)
-    if hyb[2] == 0: # No omega is assigned. LR_HF is 0 for normal Coulomb operator
-        hyb[1] = 0
     return tuple(hyb), tuple(remove_dup(fn_facs))
 
 _NAME_WITH_DASH = {'SR-HF'    : 'SR_HF',

pyscf/dft/numint.py

@@ -2717,9 +2717,12 @@ def rsh_and_hybrid_coeff(self, xc_code, spin=0):
         '''
         omega, alpha, beta = self.rsh_coeff(xc_code)
         if self.omega is not None:
+            if omega == 0 and self.omega != 0:
+                raise RuntimeError(f'Not support assigning omega={self.omega}. '
+                                   f'{xc_code} is not a RSH functional')
             omega = self.omega
 
-        if abs(omega) > 1e-10:
+        if omega != 0:
             hyb = alpha + beta
         else:
             hyb = self.hybrid_coeff(xc_code, spin)

pyscf/dft/rks.py

@@ -92,36 +92,40 @@ def get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1):
             logger.debug(ks, 'nelec with nlc grids = %s', n)
         t0 = logger.timer(ks, 'vxc', *t0)
 
+    incremental_jk = (ks._eri is None and ks.direct_scf and
+                      getattr(vhf_last, 'vj', None) is not None)
+    if incremental_jk:
+        _dm = numpy.asarray(dm) - numpy.asarray(dm_last)
+    else:
+        _dm = dm
     if not ni.libxc.is_hybrid_xc(ks.xc):
         vk = None
-        if (ks._eri is None and ks.direct_scf and
-            getattr(vhf_last, 'vj', None) is not None):
-            ddm = numpy.asarray(dm) - numpy.asarray(dm_last)
-            vj = ks.get_j(mol, ddm, hermi)
+        vj = ks.get_j(mol, _dm, hermi)
+        if incremental_jk:
             vj += vhf_last.vj
-        else:
-            vj = ks.get_j(mol, dm, hermi)
         vxc += vj
     else:
         omega, alpha, hyb = ni.rsh_and_hybrid_coeff(ks.xc, spin=mol.spin)
-        if (ks._eri is None and ks.direct_scf and
-            getattr(vhf_last, 'vk', None) is not None):
-            ddm = numpy.asarray(dm) - numpy.asarray(dm_last)
-            vj, vk = ks.get_jk(mol, ddm, hermi)
+        if omega == 0:
+            vj, vk = ks.get_jk(mol, _dm, hermi)
+            vk *= hyb
+        elif alpha == 0: # LR=0, only SR exchange
+            vj = ks.get_j(mol, _dm, hermi)
+            vk = ks.get_k(mol, _dm, hermi, omega=-omega)
             vk *= hyb
-            if omega != 0:  # For range separated Coulomb
-                vklr = ks.get_k(mol, ddm, hermi, omega=omega)
-                vklr *= (alpha - hyb)
-                vk += vklr
+        elif hyb == 0: # SR=0, only LR exchange
+            vj = ks.get_j(mol, _dm, hermi)
+            vk = ks.get_k(mol, _dm, hermi, omega=omega)
+            vk *= alpha
+        else: # SR and LR exchange with different ratios
+            vj, vk = ks.get_jk(mol, _dm, hermi)
+            vk *= hyb
+            vklr = ks.get_k(mol, _dm, hermi, omega=omega)
+            vklr *= (alpha - hyb)
+            vk += vklr
+        if incremental_jk:
             vj += vhf_last.vj
             vk += vhf_last.vk
-        else:
-            vj, vk = ks.get_jk(mol, dm, hermi)
-            vk *= hyb
-            if omega != 0:
-                vklr = ks.get_k(mol, dm, hermi, omega=omega)
-                vklr *= (alpha - hyb)
-                vk += vklr
         vxc += vj - vk * .5
 
         if ground_state:

pyscf/dft/test/test_gks.py

@@ -60,7 +60,7 @@ def test_ncol_gks_lda_omega(self):
         self.assertAlmostEqual(eks4, -75.883375491657, 6)
 
         mf = mol.GKS()
-        mf.xc = 'lda + .2*HF'
+        mf.xc = 'lda + .2*SR_HF(0.3)'
         mf.collinear = 'ncol'
         mf.omega = .5
         eks4 = mf.kernel()

pyscf/dft/test/test_grids.py

@@ -44,7 +44,12 @@ def tearDownModule():
 class KnownValues(unittest.TestCase):
     @classmethod
     def setUpClass(cls):
-        radi.ATOM_SPECIFIC_TREUTLER_GRIDS = False
+        cls.original_grids = dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS
+        dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = False
+
+    @classmethod
+    def tearDownClass(cls):
+        dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = cls.original_grids
 
     def test_gen_grid(self):
         grid = gen_grid.Grids(h2o)

pyscf/dft/test/test_he.py

@@ -85,6 +85,27 @@ def test_nr_b3lypg(self):
         m.xc = 'b3lyp5'
         self.assertAlmostEqual(m.scf(), -2.89992555753, 9)
 
+    def test_camb3lyp(self):
+        self.assertAlmostEqual(mol.RKS(xc='camb3lyp').kernel(), -2.89299475730048, 9)
+        self.assertAlmostEqual(mol.GKS(xc='camb3lyp').kernel(), -2.89299475730048, 9)
+        self.assertAlmostEqual(mol.UKS(xc='camb3lyp').kernel(), -2.89299475730048, 9)
+
+    def test_wb97(self):
+        self.assertAlmostEqual(mol.RKS(xc='wb97').kernel(), -2.89430888240579, 9)
+        self.assertAlmostEqual(mol.GKS(xc='wb97').kernel(), -2.89430888240579, 9)
+        self.assertAlmostEqual(mol.UKS(xc='wb97').kernel(), -2.89430888240579, 9)
+        # The old way to compute RSH, short-range = full-range - long-range
+        xc = 'wb97 + 1e-9*HF'
+        self.assertAlmostEqual(mol.RKS(xc=xc).kernel(), -2.89430888240579, 8)
+
+    def test_hse(self):
+        self.assertAlmostEqual(mol.RKS(xc='hse06').kernel(), -2.88908568982727, 9)
+        self.assertAlmostEqual(mol.GKS(xc='hse06').kernel(), -2.88908568982727, 9)
+        self.assertAlmostEqual(mol.UKS(xc='hse06').kernel(), -2.88908568982727, 9)
+        # The old way to compute RSH, short-range = full-range - long-range
+        xc = 'hse06 + 1e-9*HF'
+        self.assertAlmostEqual(mol.RKS(xc=xc).kernel(), -2.88908568982727, 8)
+
     def test_nr_lda_1e(self):
         mf = dft.RKS(mol1).run()
         self.assertAlmostEqual(mf.e_tot, -1.936332393935281, 9)

pyscf/dft/test/test_libxc.py

@@ -73,18 +73,22 @@ def test_parse_xc(self):
         self.assertAlmostEqual(hyb, .15, 12)
 
         hyb, fn_facs = dft.libxc.parse_xc('0.6*CAM_B3LYP+0.4*B3P86V5')
-        self.assertTrue(numpy.allclose(hyb, (.08, 0, 0)))
+        self.assertTrue(numpy.allclose(hyb, (.08, .08, 0)))
         self.assertTrue(numpy.allclose(fn_facs,
                                        ((433, 0.6), (1, 0.032), (106, 0.288), (132, 0.324), (7, 0.076))))
         rsh = dft.libxc.rsh_coeff('0.6*CAM_B3LYP+0.4*B3P86V5')
-        self.assertTrue(numpy.allclose(rsh, (0.33, 0.39, -0.196)))
+        rsh1 = dft.libxc.rsh_coeff('CAM_B3LYP')
+        hyb = dft.libxc.hybrid_coeff('B3P86V5')
+        self.assertTrue(numpy.allclose(rsh, (0.33, 0.6*rsh1[1]+.4*hyb, 0.6*rsh1[2])))
 
         hyb, fn_facs = dft.libxc.parse_xc('0.4*B3P86V5+0.6*CAM_B3LYP')
-        self.assertTrue(numpy.allclose(hyb, (.08, 0, 0)))
+        self.assertTrue(numpy.allclose(hyb, (.08, .08, 0)))
         self.assertTrue(numpy.allclose(fn_facs,
                                        ((1, 0.032), (106, 0.288), (132, 0.324), (7, 0.076), (433, 0.6))))
         rsh = dft.libxc.rsh_coeff('0.4*B3P86V5+0.6*CAM_B3LYP')
-        self.assertTrue(numpy.allclose(rsh, (0.33, 0.39, -0.196)))
+        rsh1 = dft.libxc.rsh_coeff('CAM_B3LYP')
+        hyb = dft.libxc.hybrid_coeff('B3P86V5')
+        self.assertTrue(numpy.allclose(rsh, (0.33, 0.6*rsh1[1]+.4*hyb, 0.6*rsh1[2])))
 
         hyb, fn_facs = dft.libxc.parse_xc('0.5*SR-HF(0.3) + .8*HF + .22*LR_HF')
         self.assertEqual(hyb, (1.3, 1.02, 0.3))
@@ -103,9 +107,18 @@ def test_parse_xc(self):
         self.assertEqual(hyb, (1.3, 1.02, 0.3))
         self.assertEqual(fn_facs, ((106, 0.5), (132, 0.5)))
 
+        rsh = dft.libxc.rsh_coeff('0.5*HSE06+.001*HF')
+        self.assertTrue(numpy.allclose(rsh, (0.11, .001, 0.125)))
+
+        rsh = dft.libxc.rsh_coeff('0.5*wb97+0.001*HF')
+        self.assertTrue(numpy.allclose(rsh, (0.4, 0.501, -.5)))
+
         self.assertRaises(ValueError, dft.libxc.parse_xc, 'SR_HF(0.3) + LR_HF(.5)')
         self.assertRaises(ValueError, dft.libxc.parse_xc, 'LR-HF(0.3) + SR-HF(.5)')
 
+        hyb = dft.libxc.hybrid_coeff('0.5*B3LYP+0.2*HF')
+        self.assertAlmostEqual(hyb, .3, 12)
+
         hyb = dft.libxc.hybrid_coeff('M05')
         self.assertAlmostEqual(hyb, 0.28, 9)
 
@@ -119,7 +132,7 @@ def test_parse_xc(self):
         #self.assertEqual(fn_facs, [(50, 1)])
 
         hyb, fn_facs = dft.libxc.parse_xc("9.999e-5*HF,")
-        self.assertEqual(hyb, (9.999e-5, 0, 0))
+        self.assertEqual(hyb, (9.999e-5, 9.999e-5, 0))
 
         ref = ((1, 1), (7, 1))
         self.assertEqual(dft.libxc.parse_xc_name('LDA,VWN'), (1,7))
@@ -218,6 +231,26 @@ def test_lyp(self):
     #    rho_b = numpy.array([[4.53272893e-06, 4.18968775e-06,-2.83034672e-06, 2.61832978e-06, 5.63360737e-06, 8.97541777e-07]]).T
     #    e, v = dft.libxc.eval_xc('tpss,', (rho_a, rho_b), spin=1, deriv=1)[:2]
 
+    #TDOO: enable this test when https://gitlab.com/libxc/libxc/-/issues/561 is solved
+    @unittest.skip('hse03 and hse06 fxc have large numerical errors in Libxc')
+    def test_hse06(self):
+        ni = dft.numint.NumInt()
+        rho = numpy.array([.235, 1.5e-9, 2e-9, 1e-9])[:,None]
+        xc = 'hse06'
+        fxc1 = ni.eval_xc_eff(xc, rho, deriv=2, xctype='GGA')[2]
+        rho = numpy.array([rho*.5]*2)
+        fxc2 = ni.eval_xc_eff(xc, rho, deriv=2, xctype='GGA')[2]
+        fxc2 = (fxc2[0,:,0] + fxc2[0,:,1])/2
+        self.assertAlmostEqual(abs(fxc2 - fxc1).max(), 0, 12)
+
+        rho = numpy.array([.235, 0, 0, 0])[:,None]
+        xc = 'hse06'
+        fxc1 = ni.eval_xc_eff(xc, rho, deriv=2, xctype='GGA')[2]
+        rho = numpy.array([rho*.5]*2)
+        fxc2 = ni.eval_xc_eff(xc, rho, deriv=2, xctype='GGA')[2]
+        fxc2 = (fxc2[0,:,0] + fxc2[0,:,1])/2
+        self.assertAlmostEqual(abs(fxc2 - fxc1).max(), 0, 12)
+
     def test_define_xc_gga(self):
         def eval_xc(xc_code, rho, spin=0, relativity=0, deriv=1, omega=None, verbose=None):
             # A fictitious XC functional to demonstrate the usage