try to rationalize lassi soc imag sign convention

I have done my best to nail it to the behavior of pyscf.fci.
fci_dhf_slow, which for complex h1 and h2 has the property that
e = np.tensordot (h1,dm1,axes=2) + .5*np.tensordot (h2,dm2,axes=4).

my_pyscf/lassi/op_o0.py

@@ -805,8 +805,8 @@ def root_make_rdm12s (las, ci_fr, nelec_frs, si, ix, orbsym=None, wfnsym=None):
         d1s2, d2s2 = solver.trans_rdm12s (ci_r_real, ci_r_imag, norb, nelec_r)
         d1s2 -= np.asarray (d1s2).transpose (0,2,1)
         d2s2 -= np.asarray (d2s2).transpose (0,2,1,4,3)
-        d1s -= 1j * d1s2 
-        d2s += 1j * d2s2
+        d1s += 1j * d1s2 
+        d2s -= 1j * d2s2
     rdm1s[0,:,:] = d1s[0]
     rdm1s[1,:,:] = d1s[1]
     rdm2s[0,:,:,0,:,:] = d2s[0]

my_pyscf/lassi/op_o1/rdm.py

@@ -683,8 +683,8 @@ def roots_make_rdm12s (las, ci, nelec_frs, si, **kwargs):
 
     # Put rdm1s in PySCF convention: [p,q] -> q'p
     if spin_pure: rdm1s = rdm1s.transpose (0,1,3,2)
-    else: rdm1s = rdm1s[:,0].transpose (0,2,1).conj ()
-    rdm2s = rdm2s.reshape (nroots_si, 2, 2, ncas, ncas, ncas, ncas).transpose (0,1,3,4,2,5,6)
+    else: rdm1s = rdm1s[:,0].transpose (0,2,1)
+    rdm2s = rdm2s.reshape (nroots_si, 2, 2, ncas, ncas, ncas, ncas).transpose (0,1,3,4,2,5,6).conj ()
 
     return rdm1s, rdm2s
 

tests/lassi/test_soc.py

@@ -77,7 +77,7 @@ def case_soc_stdm12s_slow (self, opt=0):
         self.assertAlmostEqual (np.amax(np.abs(dm1s_test[:,:8,8:])), 0)
     dm2_test = lib.einsum ('iabcdi->iabcd', stdm2s_test.sum ((1,4)))
     e0, h1, h2 = ham_2q (las2, las2.mo_coeff, soc=True)
-    e1 = lib.einsum ('pq,iqp->i', h1, dm1s_test)
+    e1 = lib.einsum ('pq,ipq->i', h1, dm1s_test)
     e2 = lib.einsum ('pqrs,ipqrs->i', h2, dm2_test) * .5
     e_test = e0 + e1 + e2
     with self.subTest (sanity='spin-free total energies'):
@@ -134,8 +134,8 @@ def dm_sector (dm, m):
 def case_soc_rdm12s_slow (self, opt=0):
     rdm1s_test, rdm2s_test = roots_make_rdm12s (las2, las2.ci, lsi2.si, opt=opt)
     stdm1s, stdm2s = make_stdm12s (las2, soc=True, opt=opt)    
-    rdm1s_ref = lib.einsum ('ir,jr,jabi->rab', lsi2.si.conj (), lsi2.si, stdm1s)
-    rdm2s_ref = lib.einsum ('ir,jr,isabtcdj->rsabtcd', lsi2.si.conj (), lsi2.si, stdm2s)
+    rdm1s_ref = lib.einsum ('ir,jr,iabj->rab', lsi2.si.conj (), lsi2.si, stdm1s)
+    rdm2s_ref = lib.einsum ('ir,jr,jsabtcdi->rsabtcd', lsi2.si.conj (), lsi2.si, stdm2s)
     with self.subTest (sanity='dm1s'):
         self.assertAlmostEqual (lib.fp (rdm1s_test), lib.fp (rdm1s_ref), 10)
     with self.subTest (sanity='dm2s'):
@@ -154,8 +154,7 @@ def case_soc_rdm12s_slow (self, opt=0):
     with lib.light_speed (5):
         e0, h1, h2 = ham_2q (las2, las2.mo_coeff, soc=True)
     rdm2_test = rdm2s_test.sum ((1,4))
-    # NOTE: dumbass PySCF 1-RDM convention that ket is first
-    e1 = lib.einsum ('pq,iqp->i', h1, rdm1s_test)
+    e1 = lib.einsum ('pq,ipq->i', h1, rdm1s_test)
     e2 = lib.einsum ('pqrs,ipqrs->i', h2, rdm2_test) * .5
     e_test = e0 + e1 + e2 - las2.e_states[0]
     e_ref = lsi2.e_roots - las2.e_states[0]
@@ -243,13 +242,13 @@ def test_hso (hso_test, tag='kernel'):
             stdm2 = stdm2s.sum ((1,4))
             e0eff = h0 - e0
             h0eff = np.eye (7) * e0eff
-            h1eff = lib.einsum ('pq,iqpj->ij', h1, stdm1s)
+            h1eff = lib.einsum ('pq,ipqj->ij', h1, stdm1s)
             h2eff = lib.einsum ('pqrs,ipqrsj->ij', h2, stdm2) * .5
             test_hso (h0eff + h1eff + h2eff, 'make_stdm12s')
             rdm1s, rdm2s = roots_make_rdm12s (las, las.ci, si, soc=True, break_symmetry=True,
                                               opt=opt)
             rdm2 = rdm2s.sum ((1,4))
-            e1eff = lib.einsum ('pq,iqp->i', h1, rdm1s)
+            e1eff = lib.einsum ('pq,ipq->i', h1, rdm1s)
             e2eff = lib.einsum ('pqrs,ipqrs->i', h2, rdm2) * .5
             test_hso ((si * (e0eff+e1eff+e2eff)[None,:]) @ si.conj ().T, 'roots_make_rdm12s')
         with self.subTest ('o0-o1 ham agreement'):
@@ -274,7 +273,7 @@ def test_soc_stdm12s_slow_o0 (self):
         case_soc_stdm12s_slow (self, opt=0)
 
     def test_soc_stdm12s_slow_o1 (self):
-        #case_soc_stdm12s_slow (self, opt=1)
+        case_soc_stdm12s_slow (self, opt=1)
         d_test = make_stdm12s (las2, soc=True, opt=1)
         d_ref = make_stdm12s (las2, soc=True, opt=0)
         for i,j in itertools.product (range (len (d_test[0])), repeat=2): 
@@ -287,7 +286,7 @@ def test_soc_rdm12s_slow_o0 (self):
         case_soc_rdm12s_slow (self, opt=0)
 
     def test_soc_rdm12s_slow_o1 (self):
-        #case_soc_rdm12s_slow (self, opt=1)
+        case_soc_rdm12s_slow (self, opt=1)
         d_test = roots_make_rdm12s (las2, las2.ci, lsi2.si, opt=1)
         d_ref = roots_make_rdm12s (las2, las2.ci, lsi2.si, opt=0)
         for i in range (len (d_test[0])):