logoff commit

lassi op_o1 prescreen linear dependencies not working; wrong
energies

my_pyscf/lassi/excitations.py

@@ -312,7 +312,8 @@ def get_ham_pq (self, h0, h1, h2, ci_p):
         nelec_frs = np.append (nelec_fs_p[:,None,:], nelec_frs_q, axis=1)
         with temporary_env (self, ncas_sub=norb_ref, mol=fcisolvers[0].mol):
             ham_pq, _, ovlp_pq = op[self.opt].ham (self, h1, h2, ci_fr, nelec_frs, soc=0,
-                                                   orbsym=self.orbsym_ref, wfnsym=self.wfnsym_ref)
+                                                   orbsym=self.orbsym_ref,
+                                                   wfnsym=self.wfnsym_ref)[:3]
         t1 = self.log.timer ('get_ham_pq', *t0)
         return ham_pq + (h0*ovlp_pq)
 

my_pyscf/lassi/lassi.py

@@ -356,16 +356,16 @@ def _eig_block (las, e0, h1, h2, ci_blk, nelec_blk, rootsym, soc, orbsym, wfnsym
     t0 = (lib.logger.process_clock (), lib.logger.perf_counter ())
     if (las.verbose > lib.logger.INFO) and (o0_memcheck):
         ham_ref, s2_ref, ovlp_ref = op_o0.ham (las, h1, h2, ci_blk, nelec_blk, soc=soc,
-                                               orbsym=orbsym, wfnsym=wfnsym)
+                                               orbsym=orbsym, wfnsym=wfnsym)[:3]
         t0 = lib.logger.timer (las, 'LASSI diagonalizer rootsym {} CI algorithm'.format (
             rootsym), *t0)
 
         h1_sf = h1
         if soc:
             h1_sf = (h1[0:las.ncas,0:las.ncas]
                      - h1[las.ncas:2*las.ncas,las.ncas:2*las.ncas]).real/2
-        ham_blk, s2_blk, ovlp_blk = op[opt].ham (las, h1_sf, h2, ci_blk, nelec_blk, orbsym=orbsym,
-                                               wfnsym=wfnsym)
+        ham_blk, s2_blk, ovlp_blk, get_ovlp = op[opt].ham (las, h1_sf, h2, ci_blk, nelec_blk,
+                                                           orbsym=orbsym, wfnsym=wfnsym)
         t0 = lib.logger.timer (las, 'LASSI diagonalizer rootsym {} TDM algorithm'.format (
             rootsym), *t0)
         lib.logger.debug (las,
@@ -388,8 +388,8 @@ def _eig_block (las, e0, h1, h2, ci_blk, nelec_blk, rootsym, soc, orbsym, wfnsym
     else:
         if (las.verbose > lib.logger.INFO): lib.logger.debug (
             las, 'Insufficient memory to test against o0 LASSI algorithm')
-        ham_blk, s2_blk, ovlp_blk = op[opt].ham (las, h1, h2, ci_blk, nelec_blk, soc=soc,
-                                                 orbsym=orbsym, wfnsym=wfnsym)
+        ham_blk, s2_blk, ovlp_blk, get_ovlp = op[opt].ham (las, h1, h2, ci_blk, nelec_blk, soc=soc,
+                                                           orbsym=orbsym, wfnsym=wfnsym)
         t0 = lib.logger.timer (las, 'LASSI H build rootsym {}'.format (rootsym), *t0)
     log_debug = lib.logger.debug2 if las.nroots>10 else lib.logger.debug
     if np.iscomplexobj (ham_blk):
@@ -430,7 +430,8 @@ def _eig_block (las, e0, h1, h2, ci_blk, nelec_blk, rootsym, soc, orbsym, wfnsym
         lc = 'checking if LASSI basis has lindeps: |ovlp| = {:.6e}'.format (ovlp_det)
         lib.logger.info (las, 'Caught error %s, %s', str (err), lc)
         if ovlp_det < LINDEP_THRESH:
-            raw2orth, orth2raw = citools.get_orth_basis (ci_blk, las.ncas_sub, nelec_blk)
+            raw2orth, orth2raw = citools.get_orth_basis (ci_blk, las.ncas_sub, nelec_blk,
+                                                         _get_ovlp=get_ovlp)
             xhx = raw2orth (ham_blk).conj ().T
             #x = canonical_orth_(ovlp_blk, thr=LINDEP_THRESH)
             lib.logger.info (las, '%d/%d linearly independent model states',

my_pyscf/lassi/op_o0.py

@@ -1,5 +1,6 @@
 import sys
 import numpy as np
+import functools
 from scipy import linalg
 from pyscf.fci import cistring
 from pyscf import fci, lib
@@ -585,7 +586,8 @@ def contract_h (c, nel):
         ket = None
         ham_eff[i,:] = dotter (hket, nelec_ket, spinless2ss=spinless2ss, iket=i, oporder=2)
 
-    return ham_eff, s2_eff, ovlp_eff
+    _get_ovlp = functools.partial (get_ovlp, ci_fr, norb_f, nelec_frs)
+    return ham_eff, s2_eff, ovlp_eff, _get_ovlp
 
 def contract_ham_ci (las, h1, h2, ci_fr_ket, nelec_frs_ket, ci_fr_bra, nelec_frs_bra, 
                      soc=0, orbsym=None, wfnsym=None):

my_pyscf/lassi/op_o1/hams2ovlp.py

@@ -472,6 +472,6 @@ def ham (las, h1, h2, ci, nelec_frs, soc=0, nlas=None, _HamS2Ovlp_class=HamS2Ovl
     if las.verbose >= lib.logger.TIMER_LEVEL:
         lib.logger.info (las, 'LASSI Hamiltonian crunching profile:\n%s', outerprod.sprint_profile ())
 
-    return ham, s2, ovlp
+    return ham, s2, ovlp, outerprod.get_ovlp
 
 

my_pyscf/lassi/op_o1/stdm.py

@@ -425,19 +425,19 @@ def crunch_ovlp (self, bra, ket):
 
     def get_ovlp (self, rootidx=None):
         exc_null = self.exc_null
-        offs_lroots = self.offs_lroots
+        offs_lroots = self.offs_lroots.copy ()
         nstates = self.nstates
         if rootidx is not None:
             rootidx = np.atleast_1d (rootidx)
             bra_null = np.isin (self.exc_null[:,0], rootidx)
             ket_null = np.isin (self.exc_null[:,1], rootidx)
-            idx_null = bra_null & ket_null
-            exc_null = exc_null[idx_null,:]
-            lroots = self.lroots[:,idx_null]
+            exc_null = exc_null[bra_null&ket_null,:]
+            lroots = self.lroots[:,rootidx]
             nprods = np.prod (lroots, axis=0)
             offs1 = np.cumsum (nprods)
-            offs0 = offs1 - nprods 
-            offs_lroots = np.stack ([offs0, offs1], axis=1)
+            offs0 = offs1 - nprods
+            for i, iroot in enumerate (rootidx):
+                offs_lroots[iroot,:] = [offs0[i], offs1[i]]
             nstates = offs1[-1]
         ovlp = np.zeros ([nstates,]*2, dtype=self.dtype)
         for bra, ket in exc_null:

tests/lassi/test_4frag.py

@@ -124,9 +124,9 @@ def test_stdm12s (self):
     def test_ham_s2_ovlp (self):
         h1, h2 = ham_2q (las, las.mo_coeff, veff_c=None, h2eff_sub=None)[1:]
         lbls = ('ham','s2','ovlp')
-        mats_o0 = op_o0.ham (las, h1, h2, las.ci, nelec_frs)#, orbsym=orbsym, wfnsym=wfnsym)
+        mats_o0 = op_o0.ham (las, h1, h2, las.ci, nelec_frs)[:3]#, orbsym=orbsym, wfnsym=wfnsym)
         fps_o0 = [lib.fp (mat) for mat in mats_o0]
-        mats_o1 = op_o1.ham (las, h1, h2, las.ci, nelec_frs)#, orbsym=orbsym, wfnsym=wfnsym)
+        mats_o1 = op_o1.ham (las, h1, h2, las.ci, nelec_frs)[:3]#, orbsym=orbsym, wfnsym=wfnsym)
         for lbl, mat, fp in zip (lbls, mats_o1, fps_o0):
             with self.subTest(opt=1, matrix=lbl):
                 self.assertAlmostEqual (lib.fp (mat), fp, 9)

tests/lassi/test_opt57_slow.py

@@ -158,9 +158,9 @@ def test_ham_s2_ovlp (self):
         h1, h2 = ham_2q (las, las.mo_coeff, veff_c=None, h2eff_sub=None)[1:]
         lbls = ('ham','s2','ovlp')
         t0, w0 = lib.logger.process_clock (), lib.logger.perf_counter ()
-        mats_o0 = op_o0.ham (las, h1, h2, las.ci, nelec_frs, orbsym=orbsym, wfnsym=wfnsym)
+        mats_o0 = op_o0.ham (las, h1, h2, las.ci, nelec_frs, orbsym=orbsym, wfnsym=wfnsym)[:3]
         t1, w1 = lib.logger.process_clock (), lib.logger.perf_counter ()
-        mats_o1 = op_o1.ham (las, h1, h2, las.ci, nelec_frs, orbsym=orbsym, wfnsym=wfnsym)
+        mats_o1 = op_o1.ham (las, h1, h2, las.ci, nelec_frs, orbsym=orbsym, wfnsym=wfnsym)[:3]
         t2, w2 = lib.logger.process_clock (), lib.logger.perf_counter ()
         #print (t1-t0, t2-t1)
         #print (w1-w0, w2-w1)