updated me2B for hopefully correct values

oop_imsrg/spin_sq.py

@@ -152,15 +152,26 @@ class TSpinSq(object):
 
         spin = self.s
         sisj_matrix = np.array([[spin**2,   0,    0,     0],
-                                [0,  -spin**2,  spin,     0],
-                                [0,    spin, -spin**2,     0],
+                                [0,  -(spin**2),  (0.5*spin**2),     0],
+                                [0,    (0.5*spin**2), -(spin**2),     0],
                                 [0,      0,    0, spin**2]])
 
-        sym1 = sisj_matrix[int(ps==qs),int(rs==ss)]
-        sym2 = sisj_matrix[int(qs==ps),int(ss==rs)]
-        asym1 = sisj_matrix[int(ps==qs),int(ss==rs)]
-        asym2 = sisj_matrix[int(qs==ps),int(rs==ss)]
+        bra_idx = (1-int(ps==qs))+int(qs>ps)
+        bra_idx_swap = (1-int(ps==qs))+int(ps>qs)
 
+        ket_idx = (1-int(rs==ss))+int(ss>rs)
+        ket_idx_swap = (1-int(rs==ss))+int(rs>ss)
+
+        sym1 = sisj_matrix[bra_idx, ket_idx]
+        sym2 = sisj_matrix[bra_idx_swap, ket_idx_swap]
+        asym1 = sisj_matrix[bra_idx, ket_idx_swap]
+        asym2 = sisj_matrix[bra_idx_swap, ket_idx]
+
+        # sym1 = int(ps==qs)*int(rs==ss)*sisj_matrix[0,0]
+        # sym2 = int(qs==ps)*int(ss==rs)*sisj_matrix[4,4]
+        # asym1 = sisj_matrix[int(ps==qs),int(ss==rs)]
+        # asym2 = sisj_matrix[int(qs==ps),int(rs==ss)]
+        
         me = 0.5*(int(pp==rp)*int(qp==sp)*(sym1+sym2) - int(pp==sp)*int(qp==rp)*(asym1+asym2))
 
 
@@ -193,7 +204,6 @@ class TSpinSq(object):
                 for r in bas1B:
                     for s in bas1B:
                         SS2B[p,q,r,s] += self.me2B(p,q,r,s)
-                        SS2B[p,q,r,s] += self.me2B(p,q,r,s)
 
         return (SS1B, SS2B)