From bd06007dd4129bb633aa71c4b9d25315b5d8adaf Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Mon, 9 Sep 2024 17:34:54 -0500 Subject: [PATCH] Very tight screening in lassi o1 --- debug/lassi/debug_soc.py | 6 +++--- my_pyscf/lassi/lassi.py | 21 +++++++++++++-------- my_pyscf/lassi/op_o0.py | 4 ++-- my_pyscf/lassi/op_o1/frag.py | 4 ++-- my_pyscf/lassi/op_o1/rdm.py | 5 +++-- my_pyscf/lassi/op_o1/stdm.py | 5 +++-- 6 files changed, 26 insertions(+), 19 deletions(-) diff --git a/debug/lassi/debug_soc.py b/debug/lassi/debug_soc.py index 90aa7f34..df017d33 100644 --- a/debug/lassi/debug_soc.py +++ b/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): diff --git a/my_pyscf/lassi/lassi.py b/my_pyscf/lassi/lassi.py index 3dffa4bc..c8ae29f8 100644 --- a/my_pyscf/lassi/lassi.py +++ b/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 and 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 diff --git a/my_pyscf/lassi/op_o0.py b/my_pyscf/lassi/op_o0.py index 460e44b2..b0b388ef 100644 --- a/my_pyscf/lassi/op_o0.py +++ b/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. diff --git a/my_pyscf/lassi/op_o1/frag.py b/my_pyscf/lassi/op_o1/frag.py index 245f6470..fc956a36 100644 --- a/my_pyscf/lassi/op_o1/frag.py +++ b/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) diff --git a/my_pyscf/lassi/op_o1/rdm.py b/my_pyscf/lassi/op_o1/rdm.py index 08900c68..6d41f2f3 100644 --- a/my_pyscf/lassi/op_o1/rdm.py +++ b/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 diff --git a/my_pyscf/lassi/op_o1/stdm.py b/my_pyscf/lassi/op_o1/stdm.py index 7e0d9676..cbd6f129 100644 --- a/my_pyscf/lassi/op_o1/stdm.py +++ b/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