updated with 3B flow equations

main.py

-# Main program for IM-SRG.
+# Main program for IM-SRG(2).
 
 
 # Author: Jacob Davison
@@ -104,7 +104,7 @@ def ravel(y, bas_len):
 
 # @profile
 def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
-    """Main method uses scipy.integrate.ode to solve the IMSRG flow
+    """Main method uses scipy.integrate.ode to solve the IMSRG(2) flow
     equations."""
     start = time.time()
 
@@ -120,7 +120,7 @@ def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
     initf = time.time()
     print("Initialized objects in {:2.4f} seconds\n".format(initf-initi))
 
-    print("""Pairing model IM-SRG flow:
+    print("""Pairing model IM-SRG(2) flow:
     d              = {:2.4f}
     g              = {:2.4f}
     pb             = {:2.4f}
@@ -431,22 +431,25 @@ if __name__ == '__main__':
     # test = main(4,4)
     # test2 = ci_matrix.exact_diagonalization(1.0,0.5,0.001)
     # print(test2)
-
-    h = PairingHamiltonian2B(4,4)
-    occt = OccupationTensors(h.sp_basis, h.reference)
-    wg2b = WegnerGenerator(h, occt)
-    wg3b = WegnerGenerator3B(h, occt)
-    test = wg3b.calc_eta()
-    test2 = wg2b.calc_eta()
-
-    eta1B_test = test2[0]
-    eta2B_test = test2[1]
-
-    eta1B = test[0]
-    eta2B = test[1]
-    eta3B = test[2]
-    print(eta1B)
-    print(eta2B[0,1,4,5])
-    print(eta2B[4,5,0,1])
-    print(eta2B_test[0,1,4,5])
-    print(eta3B.shape)
+    test = main3b(4,4)
+    # h = PairingHamiltonian2B(4,4)
+    # occt = OccupationTensors(h.sp_basis, h.reference)
+    # # wg2b = WegnerGenerator(h, occt)
+    # wg3b = WegnerGenerator3B(h, occt)
+    # fl = Flow_IMSRG3(h, occt)
+    # test = fl.flow(wg3b)
+    # print(test[0])
+    # test = wg3b.calc_eta()
+    # test2 = wg2b.calc_eta()
+    #
+    # eta1B_test = test2[0]
+    # eta2B_test = test2[1]
+    #
+    # eta1B = test[0]
+    # eta2B = test[1]
+    # eta3B = test[2]
+    # print(eta1B)
+    # print(eta2B[0,1,4,5])
+    # print(eta2B[4,5,0,1])
+    # print(eta2B_test[0,1,4,5])
+    # print(eta3B.shape)

oop_imsrg/flow.py

@@ -75,7 +75,9 @@ class Flow_IMSRG2(Flow):
         f = gen.f
         G = gen.G
 
-        eta1B, eta2B = gen.calc_eta()
+        partition = gen.calc_eta()
+        eta1B = partition[0]
+        eta2B = partition[1]
 
         # occA = occ_t.occA
         # occB = occ_t.occB
@@ -169,6 +171,206 @@ class Flow_IMSRG2(Flow):
         return (dE, df, dG)
 
 class Flow_IMSRG3(Flow_IMSRG2):
+    """Calculates the flow equations for the IMSRG(3). Inherits from Flow_IMSRG2."""
 
-    def __init__(self, gen):
-        super().__init__(gen)
+    def __init__(self, h, occ_t):
+        """Class constructor. Instantiates Flow_IMSRG3 object.
+
+        Arguments:
+
+        h -- Hamiltonian object
+        occ_t -- OccupationTensors object"""
+
+        assert isinstance(h, Hamiltonian), "Arg 0 must be PairingHamiltonian object"
+        assert isinstance(occ_t, OccupationTensors), "Arg 1 must be OccupationTensors object"
+
+        # self.f = h.f
+        # self.G = h.G
+
+        self._holes = h.holes
+        self._particles = h.particles
+
+        self._occA = occ_t.occA
+        self._occA2 = occ_t.occA2
+        self._occB = occ_t.occB
+        self._occC = occ_t.occC
+        self._occD = occ_t.occD
+        self._occE = occ_t.occE
+        self._occF = occ_t.occF
+        self._occG = occ_t.occG
+        self._occH = occ_t.occH
+        self._occI = occ_t.occI
+        self._occJ = occ_t.occJ
+
+    def flow(self, gen):
+        """Iterates the IMSRG3 flow equations once. Extends IMSRG2 flow function.
+
+        Arugments:
+
+        gen -- Generator object; generator produces the flow
+
+        Returns:
+
+        (dE, -- zero-body tensor
+         df, -- one-body tensor
+         dG, -- two-body tensor
+         dW) -- three-body tensor"""
+
+        assert isinstance(gen, Generator), "Arg 0 must be WegnerGenerator object"
+
+        # f = self.f
+        # G = self.G
+        f = gen.f
+        G = gen.G
+        W = gen.W
+
+        partition = super().flow(gen)
+        dE = partition[0]
+        df = partition[1]
+        dG = partition[2]
+
+        partition = gen.calc_eta()
+        eta1B = partition[0]
+        eta2B = partition[1]
+        eta3B = partition[2]
+        # print(eta3B.shape)
+        occA = self._occA
+        occA2 = self._occA2
+        occB = self._occB
+        occC = self._occC
+        occD = self._occD
+        occE = self._occE
+        occF = self._occF
+        occG = self._occG
+        occH = self._occH
+        occI = self._occI
+        occJ = self._occJ
+
+        # Calculate 0B flow equation
+        sum3_0b_1 = np.matmul(eta3B, occE)
+        sum3_0b = ncon([sum3_0b_1, W], [(0,1,2,3,4,5), (3,4,5,0,1,2)]).numpy()
+        dE += (1/18)*sum3_0b
+
+        # Calculate 1B flow equation
+        # fourth term
+        sum4_1b_1 = np.matmul(np.transpose(occD,[2,3,0,1]), G)
+        sum4_1b_2 = np.matmul(np.transpose(occD,[2,3,0,1]), eta2B)
+        sum4_1b_3 = ncon([eta3B,  sum4_1b_1], [(0,1,-1,2,3,-2),(2,3,0,1)]).numpy()
+        sum4_1b_4 = ncon([W, sum4_1b_2], [(0,1,-1,2,3,-2),(2,3,0,1)]).numpy()
+        sum4_1b = sum4_1b_3 - sum4_1b_4
+
+        # fifth term
+        sum5_1b_1 = ncon([occF, eta3B],
+                         [(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)]).numpy()
+        sum5_1b_2 = ncon([occF, W],
+                         [(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)]).numpy()
+        sum5_1b_3 = ncon([sum5_1b_1, W],
+                         [(0,1,-1,2,3,4), (2,3,4,0,1,-2)]).numpy()
+        sum5_1b_4 = ncon([sum5_1b_2, eta3B],
+                         [(0,1,-1,2,3,4), (2,3,4,0,1,-2)]).numpy()
+        sum5_1b = sum5_1b_3 - sum5_1b_4
+
+        df += 0.25*sum4_1b + (1/12)*sum5_1b
+
+        # Calculate 2B flow equation
+        # fourth term
+        sum4_2b_1 = np.matmul(-1.0*np.transpose(occA2), f)
+        sum4_2b_2 = np.matmul(-1.0*np.transpose(occA2),  eta1B)
+        sum4_2b_3 = ncon([eta3B,  sum4_2b_1], [(0,-1,-2,1,-3,-4), (1,0)]).numpy()
+        sum4_2b_4 = ncon([W, sum4_2b_2], [(0,-1,-2,1,-3,-4), (1,0)]).numpy()
+        sum4_2b = sum4_2b_3 - sum4_2b_4
+
+        #fifth term
+        sum5_2b_1 = ncon([occG, G], [(-1,-2,-4,0,1,2), (1,2,0,-3)]).numpy()
+        sum5_2b_2 = ncon([occG,  eta2B], [(-1,-2,-4,0,1,2), (1,2,0,-3)]).numpy()
+        sum5_2b_3 = ncon([eta3B,  sum5_2b_1], [(0,-1,-2,1,2,-4), (1,2,0,-3)]).numpy()
+        sum5_2b_4 = ncon([W, sum5_2b_2], [(0,-1,-2,1,2,-4), (1,2,0,-3)]).numpy()
+        sum5_2b_5 = sum5_2b_3 - sum5_2b_4
+        sum5_2b = sum5_2b_5 - np.transpose(sum5_2b_5, [3,2,0,1]) - \
+                    np.transpose(sum5_2b_5, [0,1,3,2]) + \
+                    np.transpose(sum5_2b_5, [2,3,0,1])
+
+        #sixth term
+        sum6_2b_1 = ncon([occH, W], [(-1,-2,-3,-4,0,1,2,3),(1,2,3,0,-5,-6)]).numpy()
+        sum6_2b_2 = ncon([occH, W], [(-3,-4,-5,-6,0,1,2,3),(0,-1,-2,1,2,3)]).numpy()
+        sum6_2b_3 = ncon([eta3B, sum6_2b_1], [(0,-1,-2,1,2,3), (1,2,3,0,-3,-4)]).numpy()
+        sum6_2b_4 = ncon([eta3B, sum6_2b_2], [(1,2,3,0,-3,-4), (0,-1,-2,1,2,3)]).numpy()
+        sum6_2b = sum6_2b_3 - sum6_2b_4
+
+        #seventh term
+        sum7_2b_1 = ncon([occI, W], [(-1,-2,-3,-4,0,1,2,3), (2,3,-5,0,1,-6)]).numpy()
+        sum7_2b_2 = ncon([eta3B, sum7_2b_1], [(0,1,-1,2,3,-4),(2,3,-2,0,1,-3)]).numpy()
+        sum7_2b = sum7_2b_2 - np.transpose(sum7_2b_2,[1,0,2,3]) - \
+                              np.transpose(sum7_2b_2,[0,1,3,2]) + \
+                              np.transpose(sum7_2b_2,[1,0,3,2])
+
+        dG += sum4_2b + 0.5*sum5_2b + (1/6)*sum6_2b + 0.25*sum7_2b
+
+        # Calculate 3B flow equation
+        #first, second, and third terms (one index contraction)
+
+        #terms with P(i/jk) -- line 1 and 2
+        sum1_3b_1 = ncon([eta1B, W], [(-1,0), (0,-2,-3,-4,-5,-6)]).numpy()
+        sum1_3b_2 = ncon([f, eta3B], [(-1,0), (0,-2,-3,-4,-5,-6)]).numpy()
+        sum1_3b_3 = sum1_3b_1 - sum1_3b_2
+        sum1_3b_4 = sum1_3b_3 - np.transpose(sum1_3b_3, [1,0,2,3,4,5]) - \
+                                np.transpose(sum1_3b_3, [2,1,0,3,4,5])
+
+        #terms with P(l/mn) -- line 1 and 2
+        sum1_3b_5 = ncon([eta1B, W], [(0,-4), (-1,-2,-3,0,-5,-6)]).numpy()
+        sum1_3b_6 = ncon([f, eta3B], [(0,-4), (-1,-2,-3,0,-5,-6)]).numpy()
+        sum1_3b_7 = sum1_3b_6 - sum1_3b_5
+        sum1_3b_8 = sum1_3b_7 - np.transpose(sum1_3b_7, [0,1,2,4,3,5]) - \
+                                np.transpose(sum1_3b_7, [0,1,2,5,4,3])
+
+        #terms with P(ij/k)P(l/mn) -- line 3
+        sum1_3b_9  = ncon([eta2B, G], [(-1,-2,-4,0),(0,-3,-5,-6)])
+        sum1_3b_10 = ncon([G, eta2B], [(-1,-2,-4,0),(0,-3,-5,-6)])
+        sum1_3b_11 = sum1_3b_9 - sum1_3b_10
+        sum1_3b_12 = sum1_3b_11 - np.transpose(sum1_3b_11, [0,1,2,4,3,5]) - \
+                                  np.transpose(sum1_3b_11, [0,1,2,5,4,3])
+        sum1_3b_13 = sum1_3b_12 - np.transpose(sum1_3b_12, [2,1,0,3,4,5]) - \
+                                  np.transpose(sum1_3b_12, [0,2,1,3,4,5])
+
+        sum1_3b = sum1_3b_4 + sum1_3b_8 + sum1_3b_13
+
+        #fourth term
+        sum4_3b_1 = ncon([eta2B, occB, W], [(-1,-2,2,3),(2,3,0,1),(0,1,-3,-4,-5,-6)]).numpy()
+        sum4_3b_2 = ncon([G, occB, eta3B], [(-1,-2,2,3),(2,3,0,1),(0,1,-3,-4,-5,-6)]).numpy()
+        sum4_3b_3 = sum4_3b_1 - sum4_3b_2
+        sum4_3b = sum4_3b_3 - np.transpose(sum4_3b_3, [1,0,2,3,4,5]) - \
+                              np.transpose(sum4_3b_3, [2,1,0,3,4,5])
+
+        #fifth term
+        sum5_3b_1 = ncon([eta2B, occB, W], [(2,3,-4,-5),(2,3,0,1),(-1,-2,-3,0,1,-6)]).numpy()
+        sum5_3b_2 = ncon([G, occB, eta3B], [(2,3,-4,-5),(2,3,0,1),(-1,-2,-3,0,1,-6)]).numpy()
+        sum5_3b_3 = sum5_3b_1 - sum5_3b_2
+        sum5_3b = sum5_3b_3 - np.transpose(sum5_3b_3, [0,1,2,5,4,3]) - \
+                              np.transpose(sum5_3b_3, [0,1,2,3,5,4])
+
+        #sixth term
+        sum6_3b_1 = ncon([eta2B, occA, W], [(3,-1,2,-4),(2,3,0,1),(0,-2,-3,1,-5,-6)]).numpy()
+        sum6_3b_2 = ncon([G, occA, eta3B], [(3,-1,2,-4),(2,3,0,1),(0,-2,-3,1,-5,-6)]).numpy()
+        sum6_3b_3 = sum6_3b_1 - sum6_3b_2
+        sum6_3b_4 = sum6_3b_3 - np.transpose(sum6_3b_3, [0,1,2,4,3,5]) - \
+                                np.transpose(sum6_3b_3, [0,1,2,5,4,3])
+        sum6_3b = sum6_3b_4 - np.transpose(sum6_3b_4, [1,0,2,3,4,5]) - \
+                              np.transpose(sum6_3b_4, [2,1,0,3,4,5])
+
+        #seventh term
+        sum7_3b_1 = ncon([eta3B, occJ, W], [(-1,-2,-3,3,4,5), (3,4,5,0,1,2), (0,1,2,-4,-5,-6)]).numpy()
+        sum7_3b_2 = ncon([W, occJ, eta3B], [(-1,-2,-3,3,4,5), (3,4,5,0,1,2), (0,1,2,-4,-5,-6)]).numpy()
+        sum7_3b = sum7_3b_1 - sum7_3b_2
+
+        #eighth term
+        sum8_3b_1 = ncon([eta3B, occC, W], [(3,4,-3,5,-5,-6), (3,4,5,0,1,2), (2,-1,-2,0,1,-4)]).numpy()
+        sum8_3b_2 = ncon([eta3B, occC, W], [(5,-2,-3,3,4,-6), (3,4,5,0,1,2), (-1,0,1,-4,-5,2)]).numpy()
+        sum8_3b_3 = sum8_3b_1 - sum8_3b_2
+        sum8_3b_4 = sum8_3b_3 - np.transpose(sum8_3b_3, [0,1,2,4,3,5]) - \
+                                np.transpose(sum8_3b_3, [0,1,2,5,4,3])
+        sum8_3b = sum8_3b_4 - np.transpose(sum8_3b_4, [2,1,0,3,4,5]) - \
+                              np.transpose(sum8_3b_4, [0,2,1,3,4,5])
+
+        dW = sum1_3b + 0.5*sum4_3b + (-0.5)*sum5_3b + (-1.0)*sum6_3b + (1/6)*sum7_3b + 0.5*sum8_3b
+
+        return (dE, df, dG, dW)

oop_imsrg/generator.py

@@ -185,7 +185,6 @@ class WegnerGenerator(Generator):
         return (eta1B, eta2B)
 
 class WegnerGenerator3B(WegnerGenerator):
-
     """Calculate Wegner's generator for a normal ordered Hamiltonian.
        Truncate at three-body interactions. Inherits from WegnerGenerator."""
 

oop_imsrg/occupation_tensors.py

@@ -260,7 +260,7 @@ class OccupationTensors(object):
         ref = self._reference
         n = len(bas1B)
 
-        occE = np.zeros((n,n,n,n,n,n,n))
+        occE = np.zeros((n,n,n,n,n,n))
 
         for a in bas1B:
             for b in bas1B: