Very tight screening in lassi o1

debug/lassi/debug_soc.py

@@ -41,6 +41,7 @@ def setUpModule():
     las2.state_average_(weights=[1,0,0,0,0],
                             spins=[[0,0],[2,0],[-2,0],[0,2],[0,-2]],
                             smults=[[1,1],[3,1],[3,1],[1,3],[1,3]])
+
     # NOTE: Be careful about state selection. You have to select states that can actually be coupled
     # by a 1-body SOC operator. For instance, spins=[0,0] and spins=[2,2] would need at least a 2-body
     # operator to couple.
@@ -59,7 +60,6 @@ def tearDownModule():
 
 def case_soc_stdm12s_slow (self, opt=0):
     stdm1s_test, stdm2s_test = make_stdm12s (las2, soc=True, opt=opt) 
-    np.save ('stdm1s.{}.npy'.format (opt), stdm1s_test)
     with self.subTest ('2-electron'):
         self.assertAlmostEqual (linalg.norm (stdm2s_test), 12.835690990485933, 6)
     with self.subTest ('1-electron'):
@@ -259,7 +259,7 @@ def test_soc_stdm12s_slow_o0 (self):
 
     def test_soc_stdm12s_slow_o1 (self):
         #case_soc_stdm12s_slow (self, opt=1)
-        d_test = make_stdm12s (las2, soc=True, opt=1)
+        d_test = make_stdm12s (las2, soc=True, opt=1, screen_linequiv=False)
         d_ref = make_stdm12s (las2, soc=True, opt=0)
         for i,j in itertools.product (range (len (d_test[0])), repeat=2): 
             for r in range (2):
@@ -272,7 +272,7 @@ def test_soc_rdm12s_slow_o0 (self):
 
     def test_soc_rdm12s_slow_o1 (self):
         #case_soc_rdm12s_slow (self, opt=1)
-        d_test = roots_make_rdm12s (las2, las2.ci, lsi2.si, opt=1)
+        d_test = roots_make_rdm12s (las2, las2.ci, lsi2.si, opt=1, screen_linequiv=False)
         d_ref = roots_make_rdm12s (las2, las2.ci, lsi2.si, opt=0)
         for i in range (len (d_test[0])):
             for r in range (2):

my_pyscf/lassi/lassi.py

@@ -438,7 +438,8 @@ def _eig_block (las, e0, h1, h2, ci_blk, nelec_blk, rootsym, soc, orbsym, wfnsym
         else: raise (err) from None
     return e, c, s2_blk
 
-def make_stdm12s (las, ci=None, orbsym=None, soc=False, break_symmetry=False, opt=1):
+def make_stdm12s (las, ci=None, orbsym=None, soc=False, break_symmetry=False, screen_linequiv=True,
+                  opt=1):
     ''' Evaluate <I|p'q|J> and <I|p'r'sq|J> where |I>, |J> are LAS states.
 
         Args:
@@ -505,7 +506,8 @@ def make_stdm12s (las, ci=None, orbsym=None, soc=False, break_symmetry=False, op
             d1s, d2s = op_o0.make_stdm12s (las1, ci_blk, nelec_blk, orbsym=orbsym, wfnsym=wfnsym)
             t0 = lib.logger.timer (las, 'LASSI make_stdm12s rootsym {} CI algorithm'.format (
                 sym), *t0)
-            d1s_test, d2s_test = op_o1.make_stdm12s (las1, ci_blk, nelec_blk)
+            d1s_test, d2s_test = op_o1.make_stdm12s (las1, ci_blk, nelec_blk,
+                                                     screen_linequiv=screen_linequiv)
             t0 = lib.logger.timer (las, 'LASSI make_stdm12s rootsym {} TDM algorithm'.format (
                 sym), *t0)
             lib.logger.debug (las,
@@ -523,7 +525,8 @@ def make_stdm12s (las, ci=None, orbsym=None, soc=False, break_symmetry=False, op
         else:
             if not o0_memcheck: lib.logger.debug (
                 las, 'Insufficient memory to test against o0 LASSI algorithm')
-            d1s, d2s = op[opt].make_stdm12s (las1, ci_blk, nelec_blk, orbsym=orbsym, wfnsym=wfnsym)
+            d1s, d2s = op[opt].make_stdm12s (las1, ci_blk, nelec_blk, orbsym=orbsym, wfnsym=wfnsym,
+                                             screen_linequiv=screen_linequiv)
             t0 = lib.logger.timer (las, 'LASSI make_stdm12s rootsym {}'.format (sym), *t0)
         idx_allprods.append (list(np.where(idx_prod)[0]))
         nprods += len (idx_allprods[-1])
@@ -548,7 +551,7 @@ def make_stdm12s (las, ci=None, orbsym=None, soc=False, break_symmetry=False, op
     return stdm1s, stdm2s
 
 def roots_make_rdm12s (las, ci, si, orbsym=None, soc=None, break_symmetry=None, rootsym=None,
-                       opt=1):
+                       screen_linequiv=True, opt=1):
     '''Evaluate 1- and 2-electron reduced density matrices of LASSI states
 
         Args:
@@ -628,7 +631,8 @@ def roots_make_rdm12s (las, ci, si, orbsym=None, soc=None, break_symmetry=None,
                                                 wfnsym=wfnsym)
             t0 = lib.logger.timer (las, 'LASSI make_rdm12s rootsym {} CI algorithm'.format (sym),
                                    *t0)
-            d1s_test, d2s_test = op[opt].roots_make_rdm12s (las1, ci_blk, nelec_blk, si_blk)
+            d1s_test, d2s_test = op[opt].roots_make_rdm12s (las1, ci_blk, nelec_blk, si_blk,
+                                                            screen_linequiv=screen_linequiv)
             t0 = lib.logger.timer (las, 'LASSI make_rdm12s rootsym {} TDM algorithm'.format (sym),
                                    *t0)
             lib.logger.debug (las,
@@ -647,7 +651,7 @@ def roots_make_rdm12s (las, ci, si, orbsym=None, soc=None, break_symmetry=None,
             if not o0_memcheck: lib.logger.debug (las,
                 'Insufficient memory to test against o0 LASSI algorithm')
             d1s, d2s = op[opt].roots_make_rdm12s (las1, ci_blk, nelec_blk, si_blk, orbsym=orbsym,
-                                                  wfnsym=wfnsym)
+                                                  wfnsym=wfnsym, screen_linequiv=screen_linequiv)
             t0 = lib.logger.timer (las, 'LASSI make_rdm12s rootsym {}'.format (sym), *t0)
         idx_int = np.where (idx_si)[0]
         for (i,a) in enumerate (idx_int):
@@ -658,7 +662,7 @@ def roots_make_rdm12s (las, ci, si, orbsym=None, soc=None, break_symmetry=None,
     return rdm1s, rdm2s
 
 def root_make_rdm12s (las, ci, si, state=0, orbsym=None, soc=None, break_symmetry=None,
-                      rootsym=None, opt=1):
+                      rootsym=None, screen_linequiv=True, opt=1):
     '''Evaluate 1- and 2-electron reduced density matrices of one single LASSI state
 
         Args:
@@ -703,7 +707,8 @@ def root_make_rdm12s (las, ci, si, state=0, orbsym=None, soc=None, break_symmetr
         rootsym = getattr (si, 'rootsym', getattr (las, 'rootsym', None))
     rootsym = [rootsym[s] for s in states]
     rdm1s, rdm2s = roots_make_rdm12s (las, ci, si_column, orbsym=orbsym, soc=soc,
-                                      break_symmetry=break_symmetry, rootsym=rootsym, opt=opt)
+                                      break_symmetry=break_symmetry, rootsym=rootsym,
+                                      screen_linequiv=screen_linequiv, opt=opt)
     if len (states) == 1:
         rdm1s, rdm2s = rdm1s[0], rdm2s[0]
     return rdm1s, rdm2s

my_pyscf/lassi/op_o0.py

@@ -627,7 +627,7 @@ def contract_h (c, nel):
             hket_fr_pabq[ifrag][ibra] = np.stack (hket_fr_pabq[ifrag][ibra], axis=-1)
     return hket_fr_pabq
 
-def make_stdm12s (las, ci_fr, nelec_frs, orbsym=None, wfnsym=None):
+def make_stdm12s (las, ci_fr, nelec_frs, orbsym=None, wfnsym=None, **kwargs):
     '''Build LAS state interaction transition density matrices
 
     Args:
@@ -827,7 +827,7 @@ def root_make_rdm12s (las, ci_fr, nelec_frs, si, ix, orbsym=None, wfnsym=None):
 
     return rdm1s, rdm2s
 
-def roots_make_rdm12s (las, ci_fr, nelec_frs, si, orbsym=None, wfnsym=None):
+def roots_make_rdm12s (las, ci_fr, nelec_frs, si, orbsym=None, wfnsym=None, **kwargs):
     '''Build LAS state interaction reduced density matrices for final
     LASSI eigenstates.
 

my_pyscf/lassi/op_o1/frag.py

@@ -309,13 +309,13 @@ def _init_crunch_(self, screen_linequiv):
             isequal = False
             if ci[i] is ci[j]: isequal = True
             elif np.all (ci[i]==ci[j]): isequal = True
-            elif np.all (np.abs (ci[i]-ci[j]) < 1e-8): isequal=True
+            elif np.all (np.abs (ci[i]-ci[j]) < 1e-16): isequal=True
             else:
                 ci_i = ci[i].reshape (lroots[i],-1)
                 ci_j = ci[j].reshape (lroots[j],-1)
                 ovlp = ci_i.conj () @ ci_j.T
                 isequal = np.allclose (ovlp.diagonal (), 1,
-                                       rtol=1e-10, atol=1e-10)
+                                       rtol=1e-16, atol=1e-16)
                                        # need extremely high precision on this one
                 if screen_linequiv and (not isequal):
                     u, svals, vh = linalg.svd (ovlp)

my_pyscf/lassi/op_o1/rdm.py

@@ -600,7 +600,7 @@ def make_fdm1 (iroot, ifrag):
         return fdm
     return make_fdm1
 
-def roots_make_rdm12s (las, ci, nelec_frs, si, **kwargs):
+def roots_make_rdm12s (las, ci, nelec_frs, si, screen_linequiv=True, **kwargs):
     ''' Build spin-separated LASSI 1- and 2-body reduced density matrices
 
     Args:
@@ -644,7 +644,8 @@ def roots_make_rdm12s (las, ci, nelec_frs, si, **kwargs):
 
     # First pass: single-fragment intermediates
     hopping_index, ints, lroots = frag.make_ints (las, ci, nelec_frs, nlas=nlas,
-                                                  _FragTDMInt_class=FragTDMInt)
+                                                  _FragTDMInt_class=FragTDMInt,
+                                                  screen_linequiv=screen_linequiv)
     nstates = np.sum (np.prod (lroots, axis=0))
 
     # Memory check

my_pyscf/lassi/op_o1/stdm.py

@@ -987,7 +987,7 @@ def sprint_profile (self):
         profile += '\n' + fmt_str.format ('putS', self.dt_s, self.dw_s)
         return profile
 
-def make_stdm12s (las, ci, nelec_frs, **kwargs):
+def make_stdm12s (las, ci, nelec_frs, screen_linequiv=True, **kwargs):
     ''' Build spin-separated LAS product-state 1- and 2-body transition density matrices
 
     Args:
@@ -1025,7 +1025,8 @@ def make_stdm12s (las, ci, nelec_frs, **kwargs):
         ncas = ncas * 2
 
     # First pass: single-fragment intermediates
-    hopping_index, ints, lroots = frag.make_ints (las, ci, nelec_frs, nlas=nlas)
+    hopping_index, ints, lroots = frag.make_ints (las, ci, nelec_frs, nlas=nlas,
+                                                  screen_linequiv=screen_linequiv)
     nstates = np.sum (np.prod (lroots, axis=0))
 
     # Memory check