OOP approach to IM-SRG; will develop into IMSRG3

oop_imsrg/flow.py

+from tensornetwork import *
+import numpy as np
+from hamiltonian import PairingHamiltonian
+from occupation_tensors import OccupationTensors
+from generator import WegnerGenerator
+
+class Flow_IMSRG2(object):
+
+    def __init__(self, h, occ_t):
+        assert isinstance(h, PairingHamiltonian), "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._occB = occ_t.occB
+        # self._occC = occ_t.occC
+        # self._occD = occ_t.occD
+
+    # @property
+    # def f(self):
+    #     return self._f
+
+    # @property
+    # def G(self):
+    #     return self._G
+
+    # @f.setter
+    # def f(self, f):
+    #     self._f = f
+
+    # @G.setter
+    # def G(self, G):
+    #     self._G = G
+
+    def flow(self, gen, occ_t):
+        assert isinstance(gen, WegnerGenerator), "Arg 0 must be WegnerGenerator object"
+
+        # f = self.f
+        # G = self.G
+        f = gen.f
+        G = gen.G
+
+        eta1B = gen.calc_eta1B()
+        eta2B = gen.calc_eta2B()
+
+        occA = occ_t.occA
+        occB = occ_t.occB
+        occC = occ_t.occC
+        occD = occ_t.occD
+        
+        # occA = self._occA
+        # occB = self._occB
+        # occC = self._occC
+        # occD = self._occD
+
+        # - Calculate dE/ds
+        # first term
+        sum1_0b_1 = ncon([occA, eta1B], [(0, 1, -1, -2), (0, 1)]).numpy()
+        sum1_0b = ncon([sum1_0b_1, f], [(0, 1), (1, 0)]).numpy()
+        
+        # second term
+        sum2_0b_1 = np.matmul(eta2B, occD)
+        sum2_0b = ncon([sum2_0b_1, G], [(0, 1, 2, 3), (2, 3, 0, 1)]).numpy()
+        
+        dE = sum1_0b + 0.5*sum2_0b
+        
+        
+        # - Calculate df/ds
+        # first term
+        sum1_1b_1 = ncon([eta1B, f], [(-1, 0), (0, -2)]).numpy()
+        sum1_1b_2 = np.transpose(sum1_1b_1)
+        sum1_1b = sum1_1b_1 + sum1_1b_2
+        
+        # second term (might need to fix)
+        sum2_1b_1 = ncon([eta1B, G], [(0, 1), (1, -1, 0, -2)]).numpy()
+        sum2_1b_2 = ncon([f, eta2B], [(0, 1), (1, -1, 0, -2)]).numpy()
+        sum2_1b_3 = sum2_1b_1 - sum2_1b_2
+        sum2_1b = ncon([occA, sum2_1b_3],[(-1, -2, 0, 1), (0,1)]).numpy()
+        
+        # third term
+        sum3_1b_1 = ncon([occC, G], [(-1, -2, -3, 0, 1, 2), (0, 1, 2, -4)]).numpy()
+        sum3_1b_2 = ncon([eta2B, sum3_1b_1], [(2, -1, 0, 1), (0, 1, 2, -2)]).numpy()
+        sum3_1b_3 = np.transpose(sum3_1b_2)
+        sum3_1b = sum3_1b_2 + sum3_1b_3
+        
+        df = sum1_1b + sum2_1b + 0.5*sum3_1b
+        
+        
+        # - Calculate dG/ds
+        # first term (P_ij piece)
+        sum1_2b_1 = ncon([eta1B, G], [(-1, 0), (0, -2, -3, -4)]).numpy()
+        sum1_2b_2 = ncon([f, eta2B], [(-1, 0), (0, -2, -3, -4)]).numpy()
+        sum1_2b_3 = sum1_2b_1 - sum1_2b_2
+        sum1_2b_4 = np.transpose(sum1_2b_3, [1, 0, 2, 3])
+        sum1_2b_5 = sum1_2b_3 - sum1_2b_4
+        
+        # first term (P_kl piece)
+        sum1_2b_6 = ncon([eta1B, G], [(0, -3), (-1, -2, 0, -4)]).numpy()
+        sum1_2b_7 = ncon([f, eta2B], [(0, -3), (-1, -2, 0, -4)]).numpy()
+        sum1_2b_8 = sum1_2b_6 - sum1_2b_7
+        sum1_2b_9 = np.transpose(sum1_2b_8, [0, 1, 3, 2])
+        sum1_2b_10 = sum1_2b_8 - sum1_2b_9
+        
+        sum1_2b = sum1_2b_5 - sum1_2b_10
+
+        # second term
+        sum2_2b_1 = ncon([occB,     G], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
+        sum2_2b_2 = ncon([occB, eta2B], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
+        sum2_2b_3 = ncon([eta2B,  sum2_2b_1], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
+        sum2_2b_4 = ncon([G,      sum2_2b_2], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
+        sum2_2b = sum2_2b_3 - sum2_2b_4
+
+        # third term
+        sum3_2b_1 = ncon([eta2B, G], [(0, -1, 1, -3), (1, -2, 0, -4)]).numpy()
+        sum3_2b_2 = np.transpose(sum3_2b_1, [1, 0, 2, 3])
+        sum3_2b_3 = np.transpose(sum3_2b_1, [0, 1, 3, 2])
+        sum3_2b_4 = np.transpose(sum3_2b_1, [1, 0, 3, 2])
+        sum3_2b_5 = sum3_2b_1 - sum3_2b_2 - sum3_2b_3 + sum3_2b_4
+        sum3_2b = ncon([occA, sum3_2b_5], [(0, 1, -1, -2), (0, 1, -3, -4)]).numpy()
+
+        dG = sum1_2b + 0.5*sum2_2b + sum3_2b    
+        
+        return (dE, df, dG)
\ No newline at end of file

oop_imsrg/generator.py

+from tensornetwork import *
+import numpy as np
+from hamiltonian import PairingHamiltonian
+from occupation_tensors import OccupationTensors
+
+class WegnerGenerator(object):
+
+    def __init__(self, h, occ_t):
+        
+        assert isinstance(h, PairingHamiltonian), "Arg 0 must be Hamiltonian 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._occB = occ_t.occB
+        self._occC = occ_t.occC
+        self._occD = occ_t.occD
+
+
+
+        # self._eta1B = self.eta1B
+        # self._eta2B = self.eta2B
+
+    @property
+    def f(self):
+        return self._f
+
+    @property
+    def G(self):
+        return self._G
+
+    @f.setter
+    def f(self, f):
+        self._f = f
+
+    @G.setter
+    def G(self, G):
+        self._G = G
+    
+    def __decouple_OD(self):
+        f = self.f
+        G = self.G
+        holes = self._holes
+        particles = self._particles
+
+        # - Decouple off-diagonal 1B and 2B pieces
+        fod = np.zeros(f.shape)
+        fod[np.ix_(particles, holes)] += f[np.ix_(particles, holes)]
+        fod[np.ix_(holes, particles)] += f[np.ix_(holes, particles)]
+        fd = f - fod
+
+        God = np.zeros(G.shape)
+        God[np.ix_(particles, particles, holes, holes)] += G[np.ix_(particles, particles, holes, holes)]
+        God[np.ix_(holes, holes, particles, particles)] += G[np.ix_(holes, holes, particles, particles)]
+        Gd = G - God
+
+        return (fd, fod, Gd, God)
+
+    def calc_eta1B(self):
+        fd, fod, Gd, God = self.__decouple_OD()
+
+        holes = self._holes
+        particles = self._particles
+
+        occA = self._occA
+        occC = self._occC
+
+        # - Calculate 1B generator
+        # first term
+        sum1_1b_1 = ncon([fd, fod], [(-1, 0), (0, -2)]).numpy()
+        sum1_1b_2 = np.transpose(sum1_1b_1)
+        sum1_1b = sum1_1b_1 - sum1_1b_2
+        
+        # second term
+        sum2_1b_1 = ncon([fd, God], [(0, 1), (1, -1, 0, -2)]).numpy()
+        sum2_1b_2 = ncon([fod, Gd], [(0, 1), (1, -1, 0, -2)]).numpy()
+        sum2_1b_3 = sum2_1b_1 - sum2_1b_2
+        sum2_1b = ncon([occA, sum2_1b_3],[(-1, -2, 0, 1), (0,1)]).numpy()
+        
+        # third term
+        sum3_1b_1 = ncon([occC, God], [(-1, -2, -3, 0, 1, 2), (0, 1, 2, -4)]).numpy()
+        sum3_1b_2 = ncon([Gd, sum3_1b_1], [(2, -1, 0, 1), (0, 1, 2, -2)]).numpy()
+        sum3_1b_3 = np.transpose(sum3_1b_2)
+        sum3_1b = sum3_1b_2 - sum3_1b_3
+        
+        eta1B = sum1_1b + sum2_1b + 0.5*sum3_1b
+
+        return eta1B
+
+    def calc_eta2B(self):
+        fd, fod, Gd, God = self.__decouple_OD()
+
+        holes = self._holes
+        particles = self._particles
+
+        occA = self._occA
+        occB = self._occB
+
+        # - Calculate 2B generator
+        # first term (P_ij piece)
+        sum1_2b_1 = ncon([fd, God], [(-1, 0), (0, -2, -3, -4)]).numpy()
+        sum1_2b_2 = ncon([fod, Gd], [(-1, 0), (0, -2, -3, -4)]).numpy()
+        sum1_2b_3 = sum1_2b_1 - sum1_2b_2
+        sum1_2b_4 = np.transpose(sum1_2b_3, [1, 0, 2, 3])
+        sum1_2b_5 = sum1_2b_3 - sum1_2b_4
+        
+        # first term (P_kl piece)
+        sum1_2b_6 = ncon([fd, God], [(0, -3), (-1, -2, 0, -4)]).numpy()
+        sum1_2b_7 = ncon([fod, Gd], [(0, -3), (-1, -2, 0, -4)]).numpy()
+        sum1_2b_8 = sum1_2b_6 - sum1_2b_7
+        sum1_2b_9 = np.transpose(sum1_2b_8, [0, 1, 3, 2])
+        sum1_2b_10 = sum1_2b_8 - sum1_2b_9
+        
+        sum1_2b = sum1_2b_5 - sum1_2b_10
+
+        # second term
+        sum2_2b_1 = ncon([occB, God], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
+        sum2_2b_2 = ncon([occB,  Gd], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
+        sum2_2b_3 = ncon([Gd,  sum2_2b_1], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
+        sum2_2b_4 = ncon([God, sum2_2b_2], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
+        sum2_2b = sum2_2b_3 - sum2_2b_4
+
+        # third term
+        sum3_2b_1 = ncon([Gd, God], [(0, -1, 1, -3), (1, -2, 0, -4)]).numpy()
+        sum3_2b_2 = np.transpose(sum3_2b_1, [1, 0, 2, 3])
+        sum3_2b_3 = np.transpose(sum3_2b_1, [0, 1, 3, 2])
+        sum3_2b_4 = np.transpose(sum3_2b_1, [1, 0, 3, 2])
+        sum3_2b_5 = sum3_2b_1 - sum3_2b_2 - sum3_2b_3 + sum3_2b_4
+        sum3_2b = ncon([occA, sum3_2b_5], [(0, 1, -1, -2), (0, 1, -3, -4)]).numpy()
+
+        eta2B = sum1_2b + 0.5*sum2_2b + sum3_2b
+
+        return eta2B
+
+    # @property
+    # def eta1B(self):
+    #     return self._eta1B
+    
+    # @property
+    # def eta2B(self):
+    #     return self._eta2B
+
+    # @eta1B.setter
+    # def eta1B(self, occA):
+    #     fd = self._fd
+    #     Gd = self._Gd
+    #     fod = self._fod
+    #     God = self._God
+
+    #     holes = self._holes
+    #     particles = self._particles
+
+    #     occA = self._occA
+    #     occC = self._occC
+
+    #     # - Calculate 1B generator
+    #     # first term
+    #     sum1_1b_1 = ncon([fd, fod], [(-1, 0), (0, -2)]).numpy()
+    #     sum1_1b_2 = np.transpose(sum1_1b_1)
+    #     sum1_1b = sum1_1b_1 - sum1_1b_2
+        
+    #     # second term
+    #     sum2_1b_1 = ncon([fd, God], [(0, 1), (1, -1, 0, -2)]).numpy()
+    #     sum2_1b_2 = ncon([fod, Gd], [(0, 1), (1, -1, 0, -2)]).numpy()
+    #     sum2_1b_3 = sum2_1b_1 - sum2_1b_2
+    #     sum2_1b = ncon([occA, sum2_1b_3],[(-1, -2, 0, 1), (0,1)]).numpy()
+        
+    #     # third term
+    #     sum3_1b_1 = ncon([occC, God], [(-1, -2, -3, 0, 1, 2), (0, 1, 2, -4)]).numpy()
+    #     sum3_1b_2 = ncon([Gd, sum3_1b_1], [(2, -1, 0, 1), (0, 1, 2, -2)]).numpy()
+    #     sum3_1b_3 = np.transpose(sum3_1b_2)
+    #     sum3_1b = sum3_1b_2 - sum3_1b_3
+        
+    #     eta1B = sum1_1b + sum2_1b + 0.5*sum3_1b
+
+    #     self._eta1B = eta1B
+
+    # @eta2B.setter
+    # def eta2B(self, eta2B):
+    #     fd = self._fd
+    #     Gd = self._Gd
+    #     fod = self._fod
+    #     God = self._God
+
+    #     holes = self._holes
+    #     particles = self._particles
+
+    #     occA = self._occA
+    #     occB = self._occB
+
+    #     # - Calculate 2B generator
+    #     # first term (P_ij piece)
+    #     sum1_2b_1 = ncon([fd, God], [(-1, 0), (0, -2, -3, -4)]).numpy()
+    #     sum1_2b_2 = ncon([fod, Gd], [(-1, 0), (0, -2, -3, -4)]).numpy()
+    #     sum1_2b_3 = sum1_2b_1 - sum1_2b_2
+    #     sum1_2b_4 = np.transpose(sum1_2b_3, [1, 0, 2, 3])
+    #     sum1_2b_5 = sum1_2b_3 - sum1_2b_4
+        
+    #     # first term (P_kl piece)
+    #     sum1_2b_6 = ncon([fd, God], [(0, -3), (-1, -2, 0, -4)]).numpy()
+    #     sum1_2b_7 = ncon([fod, Gd], [(0, -3), (-1, -2, 0, -4)]).numpy()
+    #     sum1_2b_8 = sum1_2b_6 - sum1_2b_7
+    #     sum1_2b_9 = np.transpose(sum1_2b_8, [0, 1, 3, 2])
+    #     sum1_2b_10 = sum1_2b_8 - sum1_2b_9
+        
+    #     sum1_2b = sum1_2b_5 - sum1_2b_10
+
+    #     # second term
+    #     sum2_2b_1 = ncon([occB, God], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
+    #     sum2_2b_2 = ncon([occB,  Gd], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
+    #     sum2_2b_3 = ncon([Gd,  sum2_2b_1], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
+    #     sum2_2b_4 = ncon([God, sum2_2b_2], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
+    #     sum2_2b = sum2_2b_3 - sum2_2b_4
+
+    #     # third term
+    #     sum3_2b_1 = ncon([Gd, God], [(0, -1, 1, -3), (1, -2, 0, -4)]).numpy()
+    #     sum3_2b_2 = np.transpose(sum3_2b_1, [1, 0, 2, 3])
+    #     sum3_2b_3 = np.transpose(sum3_2b_1, [0, 1, 3, 2])
+    #     sum3_2b_4 = np.transpose(sum3_2b_1, [1, 0, 3, 2])
+    #     sum3_2b_5 = sum3_2b_1 - sum3_2b_2 - sum3_2b_3 + sum3_2b_4
+    #     sum3_2b = ncon([occA, sum3_2b_5], [(0, 1, -1, -2), (0, 1, -3, -4)]).numpy()
+
+    #     eta2B = sum1_2b + 0.5*sum2_2b + sum3_2b
+
+    #     self._eta2B = eta2B
\ No newline at end of file

oop_imsrg/hamiltonian.py

+import numpy as np
+from tensornetwork import *
+
+class PairingHamiltonian(object):
+
+    def __init__(self, n_hole_states, n_particle_states, d=1, g=0.5, pb=0.0):
+        self._d = d
+        self._g = g
+        self._pb = pb
+        self._reference = np.append(np.ones(n_hole_states), np.zeros(n_particle_states))
+        self._holes = np.arange(n_hole_states, dtype=np.int64)
+
+        self._n_sp_states = n_hole_states + n_particle_states
+        self._particles = np.arange(n_hole_states,self.n_sp_states, dtype=np.int64)
+        self._sp_basis = np.append(self.holes, self.particles)
+
+        self._H1B, self._H2B = self.__construct()
+        self._E, self._f, self._G = self.__normal_order()
+
+    @property
+    def d(self):
+        return self._d
+
+    @property
+    def g(self):
+        return self._g
+
+    @property
+    def pb(self):
+        return self._pb
+
+    @property
+    def reference(self):
+        return self._reference
+
+    @property
+    def holes(self):
+        return self._holes
+
+    @property
+    def particles(self):
+        return self._particles
+    
+    @property
+    def sp_basis(self):
+        return self._sp_basis
+    
+    @property
+    def n_sp_states(self):
+        return self._n_sp_states
+
+    @property
+    def H1B(self):
+        return self._H1B
+    
+    @property
+    def H2B(self):
+        return self._H2B
+
+    @property
+    def E(self):
+        return self._E
+    
+    @property
+    def f(self):
+        return self._f
+    
+    @property
+    def G(self):
+        return self._G
+    
+    
+    def __delta2B(self, p,q,r,s):
+        pp = np.floor_divide(p,2)
+        qp = np.floor_divide(q,2)
+        rp = np.floor_divide(r,2)
+        sp = np.floor_divide(s,2)
+
+        ps = 1 if p%2==0 else -1
+        qs = 1 if q%2==0 else -1
+        rs = 1 if r%2==0 else -1
+        ss = 1 if s%2==0 else -1
+        
+        if pp != qp or rp != sp:
+            return 0
+        if ps == qs or rs == ss:
+            return 0
+        if ps == rs and qs == ss:
+            return -1
+        if ps == ss and qs == rs:
+            return 1
+        
+        return 0
+    
+    def __deltaPB(self, p,q,r,s):
+        pp = np.floor_divide(p,2)
+        qp = np.floor_divide(q,2)
+        rp = np.floor_divide(r,2)
+        sp = np.floor_divide(s,2)
+
+        ps = 1 if p%2==0 else -1
+        qs = 1 if q%2==0 else -1
+        rs = 1 if r%2==0 else -1
+        ss = 1 if s%2==0 else -1
+        
+        if (pp != qp and rp == sp) or (pp == qp and rp != sp):
+            if ps == qs or rs == ss:
+                return 0
+            if ps == rs and qs == ss:
+                return -1
+            if ps == ss and qs == rs:
+                return 1
+        
+        return 0
+
+    def __construct(self):
+        bas1B = self.sp_basis # get the single particle basis
+
+        # one body part of Hamiltonian is floor-division of basis index
+        # matrix elements are (P-1) where P is energy level
+        H1B = np.diag(self.d*np.floor_divide(bas1B,2))
+
+        # two body part of Hamiltonian constructed from four indices
+        # with non-zero elements defined by pairing term
+        H2B = np.zeros(np.ones(4, dtype=np.int64)*self.n_sp_states)
+        for p in bas1B:
+            for q in bas1B:
+                for r in bas1B:
+                    for s in bas1B:
+                        H2B[p,q,r,s] += self.g*0.5*self.__delta2B(p,q,r,s)
+                        H2B[p,q,r,s] += self.pb*0.5*self.__deltaPB(p,q,r,s)
+                            
+        return (H1B, H2B)
+
+    def __normal_order(self):
+        bas1B = self.sp_basis # get the single particle basis
+        H1B_t = self.H1B   # get the 1B tensor
+        H2B_t = self.H2B   # get the 2B tensor
+        holes = self.holes         # get hole states
+        particles = self.particles # get particle states
+
+        net = TensorNetwork()
+
+        # - Calculate 0B piece
+        H1B_holes = H1B_t[np.ix_(holes,holes)]
+        H2B_holes = H2B_t[np.ix_(holes,holes,holes,holes)]
+
+        ob_node0b = net.add_node(H1B_holes)
+        tb_node0b = net.add_node(H2B_holes)
+        
+        ob_ii = net.connect(ob_node0b[0],ob_node0b[1])
+        tb_ijij1 = net.connect(tb_node0b[0], tb_node0b[2])
+        tb_ijij2 = net.connect(tb_node0b[1], tb_node0b[3])
+        
+        flatten = net.flatten_edges([tb_ijij1, tb_ijij2])
+        ob_contract = net.contract(ob_ii).tensor.numpy()
+        tb_contract = 0.5*net.contract(flatten).tensor.numpy()
+
+        E = ob_contract + tb_contract
+        
+        
+        # - Calculate 1B piece
+        ob_node1b = net.add_node(H1B_t)
+        tb_node1b = net.add_node(H2B_t[np.ix_(bas1B,holes,bas1B,holes)])
+        
+        tb_ihjh = net.connect(tb_node1b[1], tb_node1b[3])
+        tb_contract = net.contract(tb_ihjh)
+        
+        f = ob_node1b.tensor.numpy() + tb_contract.tensor.numpy()
+        
+        G = H2B_t
+        
+        return (E, f, G)
+
+

oop_imsrg/main.py

+from hamiltonian import PairingHamiltonian
+from occupation_tensors import OccupationTensors
+from generator import WegnerGenerator
+from flow import Flow_IMSRG2
+from scipy.integrate import odeint, ode
+import numpy as np
+import time
+
+ha = PairingHamiltonian(4,4)
+ot = OccupationTensors(ha.sp_basis, ha.reference)
+wg = WegnerGenerator(ha, ot)
+fl = Flow_IMSRG2(ha, ot)
+
+def derivative(t, y, hamiltonian, occ_tensors, generator, flow):
+    assert isinstance(hamiltonian, PairingHamiltonian), "Arg 2 must be Hamiltonian object"
+    assert isinstance(occ_tensors, OccupationTensors), "Arg 3 must be OccupationTensors object"
+    assert isinstance(generator, WegnerGenerator), "Arg 4 must be WegnerGenerator object"
+    assert isinstance(flow, Flow_IMSRG2), "Arg 5 must be a Flow_IMSRG2 object"
+
+    E, f, G = ravel(y, hamiltonian.holes, hamiltonian.particles)
+
+    generator.f = f
+    generator.G = G
+    eta1B = generator.calc_eta1B()
+    eta2B = generator.calc_eta2B()
+    
+    dE, df, dG = flow.flow(generator, occ_tensors)
+    
+    dy = unravel(dE, df, dG)
+    
+    return dy
+
+def unravel(E, f, G):
+    unravel_E = np.reshape(E, -1)
+    unravel_f = np.reshape(f, -1)
+    unravel_G = np.reshape(G, -1)
+    
+    return np.concatenate([unravel_E, unravel_f, unravel_G], axis=0)
+
+def ravel(y, holes, particles):
+    
+    bas_len = len(np.append(holes,particles))
+    
+    ravel_E = np.reshape(y[0], ())
+    ravel_f = np.reshape(y[1:bas_len**2+1], (bas_len, bas_len))
+    ravel_G = np.reshape(y[bas_len**2+1:bas_len**2+1+bas_len**4], 
+                         (bas_len, bas_len, bas_len, bas_len))
+    
+    return(ravel_E, ravel_f, ravel_G)
+
+
+def main():
+
+    print("""Pairing model IM-SRG flow: 
+    d              = {:2d}
+    g              = {:2.4f}
+    pb             = {:2.4f}
+    SP basis size  = {:2d}
+    n_holes        = {:2d}