Reduce numerical errors in TDHF diagonalization

sunqm · Dec 17, 2024 · 159430d · 159430d
1 parent 59a8cb5
commit 159430d
Show file tree

Hide file tree

Showing 3 changed files with 11 additions and 6 deletions.
diff --git a/pyscf/grad/test/test_tdrks_grad.py b/pyscf/grad/test/test_tdrks_grad.py
@@ -60,7 +60,7 @@ def tearDownClass(cls):
         dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = cls.original_grids
 
     def test_tda_singlet_lda(self):
-        td = tdscf.TDA(mf_lda).run(nstates=nstates)
+        td = tdscf.TDA(mf_lda).run(conv_tol=1e-6, nstates=nstates)
         tdg = td.nuc_grad_method()
         g1 = tdg.kernel(td.xy[2])
         self.assertAlmostEqual(g1[0,2], -9.23916667e-02, 6)
@@ -78,7 +78,7 @@ def test_tda_triplet_lda(self):
         self.assertAlmostEqual(abs((e1[2]-e2[2])/.002 - g1[0,2]).max(), 0, 4)
 
     def test_tda_singlet_b88(self):
-        td = tdscf.TDA(mf_gga).run(nstates=nstates)
+        td = tdscf.TDA(mf_gga).run(conv_tol=1e-6, nstates=nstates)
         tdg = td.nuc_grad_method()
         g1 = tdg.kernel(state=3)
         self.assertAlmostEqual(g1[0,2], -9.32506535e-02, 6)

diff --git a/pyscf/pbc/tdscf/test/test_krhf.py b/pyscf/pbc/tdscf/test/test_krhf.py
@@ -163,10 +163,7 @@ def setUpClass(cls):
         cls.cell = cell
         cls.mf = mf
 
-        mol = molgto.Mole()
-        for key in ['verbose','output','atom','basis']:
-            setattr(mol, key, getattr(cell, key))
-        mol.build()
+        mol = cell.to_mol()
         molmf = molscf.RHF(mol).density_fit(auxbasis=mf.with_df.auxbasis).run()
         cls.mol = mol
         cls.molmf = molmf

diff --git a/pyscf/tdscf/_lr_eig.py b/pyscf/tdscf/_lr_eig.py
@@ -440,9 +440,17 @@ def eig(aop, x0, precond, tol_residual=1e-5, nroots=1, x0sym=None, pick=None,
         xt[:,:half_size] -= c.T.dot(xs[:,half_size:].conj())
         xt[:,half_size:] -= c.T.dot(xs[:,:half_size].conj())
 
+        # Remove quasi linearly dependent bases, as they cause more numerical
+        # errors in _symmetric_orth
+        xt_norm = np.linalg.norm(xt, axis=1)
+        xt_to_keep = (dx_norm > tol_residual) & (xt_norm > max(lindep**.5, tol_residual))
+        xt = xt[xt_to_keep]
+        xt /= xt_norm[xt_to_keep, None]
+
         if x0sym is None:
             xt = _symmetric_orth(xt)
         else:
+            xt_ir = xt_ir[xt_to_keep]
             xt_orth = []
             xt_orth_ir = []
             for ir in set(xt_ir):