Fix pbc.tdscf tests and update docstrings

pyscf/pbc/tdscf/test/test_rks.py

@@ -98,16 +98,23 @@ def test_tddft_triplet(self):
         ref = [ 9.09496456, 11.53650896]
         self.kernel('TDDFT', ref, singlet=False)
 
-    def test_rsh_tda(self):
+    def check_rsh_tda(self, xc, place=6):
         cell = self.cell
-        for xc in ('camb3lyp', 'wb97', 'hse03'):
-            print(xc)
-            mf = cell.RKS(xc=xc).run()
-            td = mf.TDA().run(nstates=5)
-            a, b = td.get_ab()
-            no, nv = a.shape[:2]
-            eref = np.linalg.eigvalsh(a.reshape(no*nv,-1))
-            self.assertAlmostEqual(abs(td.e[:2] - eref[:2]).max(), 0, 9)
+        mf = cell.RKS(xc=xc).run()
+        td = mf.TDA().run(nstates=5, conv_tol=1e-7)
+        a, b = td.get_ab()
+        no, nv = a.shape[:2]
+        eref = np.linalg.eigvalsh(a.reshape(no*nv,-1))
+        self.assertAlmostEqual(abs(td.e[:2] - eref[:2]).max(), 0, place)
+
+    def test_camb3lyp_tda(self):
+        self.check_rsh_tda('camb3lyp')
+
+    def test_wb97_tda(self):
+        self.check_rsh_tda('wb97')
+
+    def test_hse03_tda(self):
+        self.check_rsh_tda('hse03')
 
 
 class DiamondPBEShifted(unittest.TestCase):

pyscf/pbc/tdscf/test/test_uks.py

@@ -24,26 +24,30 @@
 
 
 def diagonalize(a, b, nroots=4):
-    a_aa, a_ab, a_bb = a
-    b_aa, b_ab, b_bb = b
-    nocc_a, nvir_a, nocc_b, nvir_b = a_ab.shape
-    a_aa = a_aa.reshape((nocc_a*nvir_a,nocc_a*nvir_a))
-    a_ab = a_ab.reshape((nocc_a*nvir_a,nocc_b*nvir_b))
-    a_bb = a_bb.reshape((nocc_b*nvir_b,nocc_b*nvir_b))
-    b_aa = b_aa.reshape((nocc_a*nvir_a,nocc_a*nvir_a))
-    b_ab = b_ab.reshape((nocc_a*nvir_a,nocc_b*nvir_b))
-    b_bb = b_bb.reshape((nocc_b*nvir_b,nocc_b*nvir_b))
-    a = np.block([[ a_aa  , a_ab],
-                  [ a_ab.T, a_bb]])
-    b = np.block([[ b_aa  , b_ab],
-                  [ b_ab.T, b_bb]])
+    a = spin_orbital_block(a)
+    b = spin_orbital_block(b, True)
     abba = np.block([[a        , b       ],
                      [-b.conj(),-a.conj()]])
     e = np.linalg.eig(abba)[0]
     lowest_e = np.sort(e[e.real > 0].real)
     lowest_e = lowest_e[lowest_e > 1e-3][:nroots]
     return lowest_e
 
+def spin_orbital_block(a, symmetric=False):
+    a_aa, a_ab, a_bb = a
+    nocc_a, nvir_a, nocc_b, nvir_b = a_ab.shape
+    a_aa = a_aa.reshape((nocc_a*nvir_a,nocc_a*nvir_a))
+    a_ab = a_ab.reshape((nocc_a*nvir_a,nocc_b*nvir_b))
+    if symmetric:
+        a_ba = a_ab.T
+    else:
+        a_ba = a_ab.conj().T
+    a_bb = a_bb.reshape((nocc_b*nvir_b,nocc_b*nvir_b))
+    a = np.block([[a_aa, a_ab],
+                  [a_ba, a_bb]])
+    return a
+
+
 class DiamondM06(unittest.TestCase):
     ''' Reproduce RKS-TDSCF results
     '''
@@ -90,13 +94,7 @@ def test_tda(self):
         ref = [11.09731427, 11.57079413]
         td = self.kernel('TDA', ref)
         a, b = td.get_ab()
-        a_aa, a_ab, a_bb = a
-        nocc_a, nvir_a, nocc_b, nvir_b = a_ab.shape
-        a_aa = a_aa.reshape((nocc_a*nvir_a,nocc_a*nvir_a))
-        a_ab = a_ab.reshape((nocc_a*nvir_a,nocc_b*nvir_b))
-        a_bb = a_bb.reshape((nocc_b*nvir_b,nocc_b*nvir_b))
-        a = np.block([[ a_aa  , a_ab],
-                      [ a_ab.T, a_bb]])
+        a = spin_orbital_block(a)
         eref = np.linalg.eigvalsh(a)
         self.assertAlmostEqual(abs(td.e[:4] - eref[:4]).max(), 0, 8)
 
@@ -107,22 +105,23 @@ def test_tdhf(self):
         eref = diagonalize(a, b)
         self.assertAlmostEqual(abs(td.e[:4] - eref[:4]).max(), 0, 8)
 
-    def test_rsh_tda(self):
+    def check_rsh_tda(self, xc, place=6):
         cell = self.cell
-        for xc in ('camb3lyp', 'wb97', 'hse03'):
-            print(xc)
-            mf = cell.UKS(xc=xc).run()
-            td = mf.TDA().run(nstates=3)
-            a, b = td.get_ab()
-            a_aa, a_ab, a_bb = a
-            nocc_a, nvir_a, nocc_b, nvir_b = a_ab.shape
-            a_aa = a_aa.reshape((nocc_a*nvir_a,nocc_a*nvir_a))
-            a_ab = a_ab.reshape((nocc_a*nvir_a,nocc_b*nvir_b))
-            a_bb = a_bb.reshape((nocc_b*nvir_b,nocc_b*nvir_b))
-            a = np.block([[ a_aa  , a_ab],
-                          [ a_ab.T, a_bb]])
-            eref = np.linalg.eigvalsh(a)
-            self.assertAlmostEqual(abs(td.e[:3] - eref[:3]).max(), 0, 8)
+        mf = cell.UKS(xc=xc).run()
+        td = mf.TDA().run(nstates=3, conv_tol=1e-7)
+        a, b = td.get_ab()
+        a = spin_orbital_block(a)
+        eref = np.linalg.eigvalsh(a)
+        self.assertAlmostEqual(abs(td.e[:3] - eref[:3]).max(), 0, place)
+
+    def test_camb3lyp_tda(self):
+        self.check_rsh_tda('camb3lyp')
+
+    def test_wb97_tda(self):
+        self.check_rsh_tda('wb97')
+
+    def test_hse03_tda(self):
+        self.check_rsh_tda('hse03')
 
 
 class WaterBigBoxPBE(unittest.TestCase):
@@ -232,13 +231,7 @@ def test_tda(self):
         ref = [5.37745381, 5.37745449]
         td = self.kernel('TDA', ref)
         a, b = td.get_ab()
-        a_aa, a_ab, a_bb = a
-        nocc_a, nvir_a, nocc_b, nvir_b = a_ab.shape
-        a_aa = a_aa.reshape((nocc_a*nvir_a,nocc_a*nvir_a))
-        a_ab = a_ab.reshape((nocc_a*nvir_a,nocc_b*nvir_b))
-        a_bb = a_bb.reshape((nocc_b*nvir_b,nocc_b*nvir_b))
-        a = np.block([[ a_aa  , a_ab],
-                      [ a_ab.T, a_bb]])
+        a = spin_orbital_block(a)
         eref = np.linalg.eigvalsh(a)
         self.assertAlmostEqual(abs(td.e[:2] - eref[:2]).max(), 0, 8)
 

pyscf/scf/hf.py

@@ -961,9 +961,11 @@ def get_jk(mol, dm, hermi=1, vhfopt=None, with_j=True, with_k=True, omega=None):
             Whether to compute K matrices
 
         omega : float
-            Parameter of range-separated Coulomb operator: erf( omega * r12 ) / r12.
-            If specified, integration are evaluated based on the long-range
-            part of the range-separated Coulomb operator.
+            Parameter of range-separated Coulomb operator.
+            When omega is 0 (or None), integrals are computed with the full-range Coulomb potential.
+            When it is larger than zero, integrals are evaluated with the long-range
+            Coulomb potential erf( omega * r12 ) / r12. When omega is smaller
+            than 0, short-range Coulomb potential erfc( omega * r12 ) / r12 is applied.
 
     Returns:
         Depending on the given dm, the function returns one J and one K matrix,