lassi op_o1 hci 3frag 2c debug

debug/lassi/debug_opt57_slow.py

@@ -233,7 +233,9 @@ def test_contract_hlas_ci (self):
                         #if opt>0 and not spaces[r].is_single_excitation_of (spaces[s]): continue
                         #elif opt==1: print (r,s, round (lib.fp (hket_pq_s)-lib.fp (hket_ref_s),3))
                         with self.subTest (opt=opt, frag=f, bra_space=r, ket_space=s,
-                                           intyp=interactions[interidx[r,s]]):
+                                           intyp=interactions[interidx[r,s]],
+                                           dneleca=nelec[:,r,0]-nelec[:,s,0],
+                                           dnelecb=nelec[:,r,1]-nelec[:,s,1]):
                             self.assertAlmostEqual (lib.fp (hket_pq_s), lib.fp (hket_ref_s), 8)
 
 

my_pyscf/lassi/op_o1/hci.py

@@ -34,13 +34,7 @@ def __init__(self, ints, nlas, hopping_index, lroots, h1, h2, nbra=1,
     get_ham_2q = hams2ovlp.HamS2Ovlp.get_ham_2q
 
     # Handling for 1s1c: need to do both a'.sm.b and b'.sp.a explicitly
-    def make_exc_tables (self, hopping_index):
-        exc = super ().make_exc_tables (hopping_index)
-        exc_1s1c  = np.pad (exc['1s1c'],   ((0,0),(0,1)), constant_values=0)
-        exc_1s1cT = np.pad (exc['1s1c_T'], ((0,0),(0,1)), constant_values=1)
-        exc['1s1c'] = np.append (exc_1s1c, exc_1s1cT, axis=0)
-        return exc
-
+    all_interactions_full_square = True
     interaction_has_spin = ('_1c_', '_1c1d_', '_1s1c_', '_2c_')
 
     def _init_vecs (self):
@@ -292,8 +286,8 @@ def _crunch_2c_(self, bra, ket, i, j, k, l, s2lt):
         s12 = s2 % 2
         nelec_f_bra = self.nelec_rf[self.rootaddr[bra]]
         nelec_f_ket = self.nelec_rf[self.rootaddr[ket]]
-        # TODO: debug this for 3- and 4-fragment interactions
-        fac = 1 / (1 + int (s11==s12))
+        # TODO: debug this for 4-fragment interactions
+        fac = 1 / (1 + int (s11==s12 and i==k and j==l))
         h_iklj = self.get_ham_2q (i,j,k,l).transpose (0,2,3,1) # Dirac order
         if s11==s12 and i!=k and j!=l: # exchange
             h_iklj -= self.get_ham_2q (i,l,k,j).transpose (0,2,1,3)
@@ -319,19 +313,19 @@ def _crunch_2c_(self, bra, ket, i, j, k, l, s2lt):
             if iad: hci_f_ab[i] += fac * self.ints[i].contract_h20 (s2lt, h20_ik, ket)
         else:
             if iad:
-                h10_i = np.tensordot (h20_ik, self.ints[k].get_1_p (bra, ket, s12), axes=1)
+                h10_i = np.dot (h20_ik, self.ints[k].get_1_p (bra, ket, s12))
                 hci_f_ab[i] += fac * self.ints[i].contract_h10 (s11, h10_i, None, ket)
             if kad:
-                h10_k = np.tensordot (self.ints[i].get_1_p (bra, ket, s11), h20_ik, axes=1)
+                h10_k = np.dot (self.ints[i].get_1_p (bra, ket, s11), h20_ik)
                 hci_f_ab[k] += fac * self.ints[k].contract_h10 (s12, h10_k, None, ket)
         if j == l:
             if jad: hci_f_ab[j] += fac * self.ints[j].contract_h02 (s2lt, h02_lj, ket)
         else:
             if jad:
-                h01_j = np.tensordot (h02_jl, self.ints[l].get_1_h (bra, ket, s12), axes=1)
+                h01_j = np.dot (self.ints[l].get_1_h (bra, ket, s12), h02_lj)
                 hci_f_ab[j] += fac * self.ints[j].contract_h01 (s11, h01_j, None, ket)
             if lad:
-                h01_l = np.tensordot (self.ints[j].get_1_h (bra, ket, s11), h02_jl, axes=1)
+                h01_l = np.dot (h02_lj, self.ints[j].get_1_h (bra, ket, s11))
                 hci_f_ab[l] += fac * self.ints[l].contract_h01 (s12, h01_l, None, ket)
         dt, dw = logger.process_clock () - t0, logger.perf_counter () - w0
         self.dt_2c, self.dw_2c = self.dt_2c + dt, self.dw_2c + dw

my_pyscf/lassi/op_o1/stdm.py

@@ -253,6 +253,13 @@ def make_exc_tables (self, hopping_index):
             ispin[idx]]
         ).T
 
+        # Two-electron interaction: ij -> kk ("coalesce")
+        idx = idx_2e & (ncharge_index == 3) & (np.amax (charge_index, axis=0) == 2)
+        if nfrags > 2: exc_coalesce = np.vstack (
+            list (np.where (idx)) + [findf[-1][idx], findf[0][idx], findf[-1][idx], findf[1][idx],
+            ispin[idx]]
+        ).T
+
         # Two-electron interaction: k(a)j(b) -> i(a)k(b) ("1s1c")
         idx = idx_2e & (nspin_index==3) & (ncharge_index==2) & (np.amin(spin_index,axis=0)==-2)
         if nfrags > 2: exc['1s1c'] = np.vstack (
@@ -262,7 +269,7 @@ def make_exc_tables (self, hopping_index):
         # Two-electron interaction: k(b)j(a) -> i(b)k(a) ("1s1c_T")
         # This will only be used when we are unable to restrict ourselves to the lower triangle
         idx = idx_2e & (nspin_index==3) & (ncharge_index==2) & (np.amax(spin_index,axis=0)==2)
-        if nfrags > 2: exc['1s1c_T'] = np.vstack (
+        if nfrags > 2: exc_1s1cT = np.vstack (
             list (np.where (idx)) + [findf[-2][idx], findf[0][idx], findf[-1][idx]]
         ).T
 
@@ -289,22 +296,28 @@ def make_exc_tables (self, hopping_index):
             ispin[idx]]
         ).T
 
+        if self.all_interactions_full_square and nfrags > 2:
+            exc['1s1c'] = np.append (
+                np.pad (exc['1s1c'], ((0,0),(0,1)), constant_values=0),
+                np.pad (exc_1s1cT,   ((0,0),(0,1)), constant_values=1),
+                axis=0)
+            exc_split = np.append (exc_split, exc_coalesce, axis=0)
+
         # Combine "split", "pair", and "scatter" into "2c"
         if nfrags > 1: exc['2c'] = exc_pair
         if nfrags > 2: exc['2c'] = np.vstack ((exc['2c'], exc_split))
         if nfrags > 3: exc['2c'] = np.vstack ((exc['2c'], exc_scatter))
 
         return exc
 
+    all_interactions_full_square = False
     interaction_has_spin = ('_1c_', '_1c1d_', '_2c_')
 
     def mask_exc_table_(self, exc, lbl, mask_bra_space=None, mask_ket_space=None):
-        # Part 1: restrict to the caller-specified rectangle OR its transpose
-        idx1  = mask_exc_table (exc, col=0, mask_space=mask_bra_space)
-        idx1 &= mask_exc_table (exc, col=1, mask_space=mask_ket_space)
-        #idx2  = mask_exc_table (exc, col=1, mask_space=mask_bra_space)
-        #idx2 &= mask_exc_table (exc, col=0, mask_space=mask_ket_space)
-        exc = exc[idx1]#|idx2]
+        # Part 1: restrict to the caller-specified rectangle
+        idx  = mask_exc_table (exc, col=0, mask_space=mask_bra_space)
+        idx &= mask_exc_table (exc, col=1, mask_space=mask_ket_space)
+        exc = exc[idx]
         # Part 2: identify interactions which are equivalent except for the overlap
         # factor of spectator fragments. Reduce the exc table only to the unique
         # interactions and populate self.nonuniq_exc with the corresponding