In [32]:
# <!-- collapse=True -->
from IPython.display import display, Math, Latex 
%matplotlib inline

from numpy import *
from itertools import *
from matplotlib.pyplot import *
from copy import copy

class Combiner():
    #Normal ordered operator for cluster algebra (diagrammatic)
    def __init__(self, H, T):       
        self.q_a = H
        self.I = []
        HN = len(self.q_a)
        self.combine(H,T)
        
    def combine(self, H, ops, excitation = None):
        #Assosciate a list of T-operators (ops) to the current operator instance    
        #Find all possible ways to combine the operators using self.distinct_combinations() 
        T = []
        for i in ops:
            T.append(i)

        self.T_operator = T
        
        #Finding acceptable combinations of internal contractions between the operators
        self.I = self.distinct_combinations(self.q_a, T)
               
    def distinct_combinations(self,H,T):
        #Returns all possible combinations of H and T
        #I - all q-particle annihilation operators
        #T list of list with T operators. ex. [[-1,1],[-1,1]]  = T_1 T_1
       lenH  = len(H)
       lenT  = len(T)
       lenTi = []
       for i in range(lenT):
           lenTi.append(len(T[i]))
       #H+=[0 for i in range(lenT-1)] #adding zeros to denote separations in the cluster-operators

       H2 = []
       for i in H:
           H2.append(i)
       H2 += [0 for i in range(lenT-1)] #adding zeros to denote separations in the cluster-operators

       #Creating countlist for T-operator to keep track of q-operators in each cluster operator
       T_budget = self.itemcount(T)
       
       
       #Create all permutations
       H_permuted = permutations(H2)
       
       #Sort out indistinct diagrams and cancelling terms
       excluded = []
       excluded_budgets = []
       accepted = []
       for i in H_permuted:
           if self.acceptable(i, T, excluded, excluded_budgets):
               #print "Accepted"
               accepted.append(self.splitlist(i,0))       
       self.combined = accepted
               
           
       return accepted
       #def evaluate_q_c_total(self):
                
    def assess_excitation(self):
        #Assess current excitation level (also if combined operator)
        self.E = self.q_c.count(-1) + self.q_c.count(1) - (self.q_a.count(1) + self.q_a.count(-1))
        for i in range(len(self.T_operator)):
            self.E += self.T_operator[i].count(1) + self.T_operator[i].count(-1) 
        #fill in contracted elements
        self.E/= 2.0
    def scan_extract(L, e):
        #Exctract element e from list L
        ret = None
        for i in range(len(L)):
            if L[i] == e:
                L = delete(L, i)
                ret = e
    def nloops(self,x,y):
        #returns number of loops in budget
        return (x+y - abs(x-y))//2
    
  
    def nozeroedges(self,i):
        #Assert that there are no borders in the endpoints of the connection pattern i
        ret = True
        try:
            if i[0] == 0 or i[-1] == 0:
                ret = False
        except:
            ret = False
        return ret
    
    def nozerocontact(self,i):
        #Assert that there are no neighbouring borders in the connection pattern i
        ret = True
        for e in range(len(i) - 1):
            if i[e+1] == 0 and i[e] == 0:
                ret = False
        return ret
    
    def splitlist(self,L, d):
        #Returns a split list HL from a list L into constituents, d denotes barrier
        HL = [[]]
        n = 0
        for i in range(len(L)):
            if L[i] != d:
                HL[n].append(L[i])        
            if L[i] == d:
                HL.append([])
                n += 1        
        return HL
    
    def itemcount(self,T):
        #Count number of particle- and hole lines in each constituent part of T
        #Returns a list object of the form [[#particles, #holes], ...]
        itemnumber = []
        for i in range(len(T)):
            itemnumber.append([])
            itemnumber[i].append(T[i].count(1))  #number of q-particles
            itemnumber[i].append(T[i].count(-1)) #number of q-holes
        return itemnumber
    
    def contractable(self,L,T):
        #Asserts that the number of contractions in each p- and hline of L does not superseed # of p-h in T 
        #input two itemcount items lists, returns bool
        ret = True
        for i in range(len(T)):
            for e in range(len(T[i])):
                if T[i][e]<L[i][e]:
                    ret = False
        return ret
        
    def find_identical(self,T):
        #Find identical operators in the list of operators T
        #returns a list of pairs of indices that denotes permutations that does not alter the T operator
        identicals = []
        for i in range(len(T)):
            for e in range(i,len(T)):
                if T[i] == T[e] and i!=e:
                    identicals.append([i,e])
        return identicals
                    
    def permute_elements(self,e1,e2,L):
        #Returns a list where elements at indices e1,e2 in L is permuted
        L_ret = []        
        for i in L:
            L_ret.append(i)
        L_ret[e1] = L[e2]
        L_ret[e2] = L[e1]
        return L_ret
    
    def acceptable(self,i, T, excluded, excluded_budgets):
        #Test if a potential connection pattern is distinct
        #Returns bool
        ret = False
        identicals = self.find_identical(T)
        T_budget   = self.itemcount(T)
        
        if self.nozeroedges(i) and self.nozerocontact(i):
            I = self.splitlist(i, 0)
            I_budget = self.itemcount(I)
            if I_budget not in excluded_budgets:
                excluded_budgets.append(I_budget)
                for e in identicals:
                    excluded_budgets.append(self.permute_elements(e[0], e[1], I_budget))
                if self.contractable(I_budget,T_budget):
                    if I not in excluded:
                        excluded.append(I)
                        for e in identicals:
                            excluded.append(self.permute_elements(e[0], e[1], I))
                        ret = True
        return ret

class qnode():
    #Object for keeping track of vertices
    def __init__(self, qa):
        #self.pos = pos
        self.connected = 0
        self.qa = qa #creation or annihilation +1/-1
        self.group = []
    def connect(i):
        self.connected = i

class qop():
    #a collection of vertices -- in principle an operator
    def __init__(self,qnodes):
        self.qnodes = qnodes
        self.pos = [0,0]
        self.setpos(self.pos)
    def setpos(self, pos):
        self.pos = pos
        for i in range(len(self.qnodes)):
            self.qnodes[i].qpos = self.pos[0] + i
            

class qvert():
    #A vertex in the operator with input, output
    def __init__(self, qin, qout):
        self.qin = qin   #qnode in
        self.qout = qout #qnode out
        self.qpos = 0.0


def qlist(o_cop):
    O = o_cop
    excluded = []
    olist = []
    for i in range(len(O)):
        if O[i] == +1:
            #find -1
            for e in range(len(O)):
                if e!=i:
                    if O[e] == -1:
                        if e not in excluded:
                            if i not in excluded:
                                excluded.append(i)
                                excluded.append(e)
                                olist.append(qvert(qnode(1), qnode(-1)))
    return qop(olist)

    
def connect(H,T,I):
    #connect H and T using I as pattern
    #Return a mapping of the connections (a list of H-connections to T)
    connected_H = []
    connected_T = []    
    map_T = []
    map_H = []
    for i in range(len(H)):
        map_H.append(None)
    for e in range(len(T)):
        map_T.append([])
        for u in range(len(T[e])):
            map_T[e].append(None)
    
    for i in range(len(I)):
        for e in range(len(I[i])):
            t = I[i][e] #connect this quantity to H and T
            scan = True
            for e1 in range(len(T[i])):
                if T[i][e1] == t and map_T[i][e1] == None:
                    for i1 in range(len(H)):
                        if H[i1] == t and map_H[i1] == None:
                            if scan:
                                #Perform connection
                                map_H[i1] = [i,e1]
                                map_T[i][e1] = i1
                                scan = False
    return map_T


class Operator():
    def __init__(self, q_a, q_c, pos = [0,0]):
        self.q_a = q_a
        self.q_c = q_c
        self.internals = len(q_c)
        self.children = []
        self.E = self.excitation()
        self.pos = pos
        self.pair_up_lines()
    def excitation(self):
        #Evaluate operator excitation level
        Ex = (self.q_c.count(-1) + self.q_c.count(1) - (self.q_a.count(1) + self.q_a.count(-1)))/2.0
        for i in range(len(self.children)):
            Ex += self.children[i].excitation()
        return Ex
    def adopt(self, O, I):
        self.children = O
        self.relation = I
        for i in self.children:
            i.set_pos([self.pos[0], self.pos[1]-1])
        #for i in self.relation:
        #    if 
        #append new child operator to this operator, include more outgoing lines, use connection pattern I
        #C = Combiner(self.q_c, O.q_a) #This must be done externally
        return 0
    def get_schematic(self):
        schema = []
        for i in self.children:
            schema.append(i.get_schematic())
        #return a list of connected nodes and type of connections
        
    def set_pos(self, pos):
        self.pos = pos
        for i in self.children:
            i.pos =[pos[0], pos[1]-1]
    def get_c_pos(self,ci):
        ret = None
        if ci<self.internals:
            for i in range(len(self.pos_group)):
                if ci in self.pos_group[i][0]:
                    ret = [self.pos[0]+i, self.pos[1]]
                    break
        else:
            ret = 0
            
        return ret
        
    def pair_up_lines(self):
        #Pair up ***internal*** vertices
        p0 = 0 #Vertical position for iterations
        #Set up a mapping
        pos_groups = [] #Grouping 2 and 2 nodes together
        paired = [[],[]]
        #for i in range((len(self.q_a)+len(self.q_c))/2):
        #    pos_groups.append([])
        
        #Use same principle as below (schematic) to pair them up and map them out
        self.m_qa = []
        self.m_qc = []

        self.inout = [[],[]] #keep track of vertices in operator

        for i in range(len(self.q_c)):
            self.m_qc.append(None)
        for i in range(len(self.q_a)):
            self.m_qa.append(None)

        
        #identify pairs passing through the operator
        for i in range(self.internals):
            for e in range(len(self.q_a)):
                if self.q_a[e] == self.q_c[i] and self.m_qc[i] == None and self.m_qa[e] == None:
                    #paired[0].append(i)
                    #paired[1].append(e)
                    #pos_groups.append([[i],[e]])
                    self.m_qc[i] = e+1
                    self.m_qa[e] = i+1
                    break
                    
        #identify pairs beneath the interaction
        for i in range(len(self.q_a)):
            for e in range(len(self.q_a)):
                if i!=e and self.q_a[e] == -self.q_a[i] and self.m_qa[i] == None and self.m_qa[e] == None:
                    self.m_qa[i] = -e-1 #negative index implies index inside self
                    self.m_qa[e] = -i-1
                    break


        #identify pairs above the interaction
        for i in range(self.internals):
            for e in range(self.internals):
                if i!=e and self.q_c[e] == -self.q_c[i] and self.m_qc[i] == None and self.m_qc[e] == None:
                    self.m_qc[i] = -e-1
                    self.m_qc[e] = -i-1
                    break                                   
        
        

    def get_paired_c(self, p):
        #return the corresponding index of the line below the vertex
        ret = None
        for i in range(len(self.pos_group)):
            if p in self.pos_group[i][0]:
                ret = self.pos_group[i][1][0]
        return ret

class contracted():
    def __init__(self, H, T, pattern):
        self.H = H
        self.T = T
        self.pattern = pattern
        
        

def contract(H,T):
    T_q_c = []
    for i in T:
        i.set_pos([H.pos[0], H.pos[1]-1])
        T_q_c.append(i.q_c)
    C = Combiner(H.q_a,T_q_c) #Combines annihilators from H (Hqa) and creators from T (Tqc)
    diagrams = []
    for i in C.I:
        diagrams.append([copy(H), copy(T), connect(H.q_a, T_q_c, i), i])
    #Pending: optional combinations of operators (T^-1 H T) for property calculations
    return diagrams


class schematic():
    def __init__(self, diagram):
        self.H = diagram[0]            
        self.T = diagram[1] 
        self.schema = diagram[2]
        self.i = diagram[3]
    
    def scan_extract(L, e):
        #Exctract element e from list L
        ret = None
        for i in range(len(L)):
            if L[i] == e:
                L = delete(L, i)
                ret = e
    def nloops(self,x,y):
        #returns number of loops in budget
        return (x+y - abs(x-y))//2
    def itemcount(self,T):
        #Count number of particle- and hole lines in each constituent part of T
        #Returns a list object of the form [[#particles, #holes], ...]
        itemnumber = []
        for i in range(len(T)):
            itemnumber.append([])
            itemnumber[i].append(T[i].count(1))  #number of q-particles
            itemnumber[i].append(T[i].count(-1)) #number of q-holes
        return itemnumber
    
    def contractable(self,L,T):
        #Asserts that the number of contractions in each p- and hline of L does not superseed # of p-h in T 
        #input two itemcount items lists, returns bool
        ret = True
        for i in range(len(T)):
            for e in range(len(T[i])):
                if T[i][e]<L[i][e]:
                    #print T[i][e],L[i][e]
                    ret = False
        return ret
        
    def find_identical(self,T):
        #Find identical operators in the list of operators T
        #returns a list of pairs of indices that denotes permutations that does not alter the T operator
        identicals = []
        for i in range(len(T)):
            for e in range(i,len(T)):
                if T[i] == T[e] and i!=e:
                    identicals.append([i,e])
        return identicals
                    
    def permute_elements(self,e1,e2,L):
        #Returns a list where elements at indices e1,e2 in L is permuted
        L_ret = []        
        for i in L:
            L_ret.append(i)
        L_ret[e1] = L[e2]
        L_ret[e2] = L[e1]
        return L_ret
    
    def acceptable(self,i, T, excluded, excluded_budgets):
        #Test if a potential connection pattern is distinct
        #Returns bool
        ret = False
        identicals = self.find_identical(T)
        T_budget   = self.itemcount(T)
        
        if self.nozeroedges(i) and self.nozerocontact(i):
            I = self.splitlist(i, 0)
            I_budget = self.itemcount(I)
            if I_budget not in excluded_budgets:
                excluded_budgets.append(I_budget)
                for e in identicals:
                    excluded_budgets.append(self.permute_elements(e[0], e[1], I_budget))
                if self.contractable(I_budget,T_budget):
                    if I not in excluded:
                        excluded.append(I)
                        for e in identicals:
                            excluded.append(self.permute_elements(e[0], e[1], I))
                        ret = True
        return ret
    def setup(self):
        plabels_t = ['r','s','x','z','b','a']
        hlabels_t = ['t','u','q', 'w', 'j','i']

        plabels = ['h','g','f','e','d','c']
        hlabels = ['p','o','n','m','l','k']

                        
        #plabels_ex =  ['h','g','f','e','d','c']
        #hlabels_ex =  ['p','o','n','m','l','k']

        self.T_labels = []
        for i in range(len(self.T)):
            self.T_labels.append([])
            for e in range(len(self.T[i].q_c)):
                self.T_labels[i].append(None)

        self.H_labels_c = []
        self.H_labels_a = []        
        for i in range(len(self.H.q_a)):
            self.H_labels_a.append(None)
        for i in range(len(self.H.q_c)):
            self.H_labels_c.append(None)

        #Label all external lines exiting the cluster operators 
        for i in range(len(self.schema)):
            for e in range(len(self.schema[i])):
                if self.schema[i][e] == None:
                    if self.T[i].q_c[e] == 1:
                        self.T_labels[i][e] = plabels_t.pop()
                    if self.T[i].q_c[e] == -1:
                        self.T_labels[i][e] = hlabels_t.pop()
                
        #labels all external lines exiting the interaction
        for i in range(len(self.H.q_c)):
            
            if self.H.q_c[i] == 1:
                self.H_labels_c[i] = plabels_t.pop()
            if self.H.q_c[i] == -1:
                self.H_labels_c[i] = hlabels_t.pop()
                
        
        #label all internal lines
        for i in range(len(self.schema)):
            n = len(self.schema[i])
            for e in range(n):
                if self.schema[i][e] != None:
                    #Current index is connected
                    if self.T[i].q_c[e]==1:
                        self.T_labels[i][e] = plabels.pop()
                    if self.T[i].q_c[e]==-1:
                        self.T_labels[i][e] = hlabels.pop()
                    self.H_labels_a[self.schema[i][e]] = self.T_labels[i][e]
        
        #Identify the cluster-operators
        self.t_operators()
        
        #Identify the interaction operator
        self.v_operator()
        
        
        #Identify summation indices
        self.summation = ""
        for i in range(len(self.schema)):
            for e in range(len(self.schema[i])):
                if self.schema[i][e] != None:
                    self.summation += self.T_labels[i][e]


        #Connect the vertices in H_labels
        #Identify vertices in all operators

        
        
        
        self.count_loops()
        self.count_holes()
        self.count_equivalent_l()
        self.count_equivalent_t()
        self.calc_prefactor()
        self.find_permutations()


    def cancel_permutations(self):
        return 0
        
        
    def find_permutations(self):
        #Identify name, origin and orientation of outgoing lines
        orig = []
        for i in range(len(self.schema)):
            for e in range(len(self.schema[i])):
                if self.schema[i][e] == None:
                    orig.append([self.T_labels[i][e], self.T[i].q_c[e], "t"+str(e)])
        for i in range(len(self.H.q_c)):
            orig.append([self.H_labels_c[i], self.H.q_c[i], "v"])

        permutes = []
        #compare lines
        for i in range(len(orig)):
            for e in range(i,len(orig)):
                if i!=e and self.ineq_external(orig[i],orig[e]):
                    permutes.append(str(orig[i][0])+str(orig[e][0]))
        self.permutations = permutes
                    
                    
        return 0
    def ineq_external(self, line1, line2):
        #compare two external lines, decide wether they are inequivalent (permutations)
        ret = False
        if line1[1] == line2[1]:
            if line1[2] != line2[2]:
                ret = True
        return ret
        
    def calc_prefactor(self):
        self.prefactor = (0.5**self.n_equivalent_Ts)*(0.5**self.equiv_lines)*((-1)**(self.holes-self.loops))
    def count_equivalent_t(self):
        #count equivalent T-vertices
        t_config = []
        for i in range(len(self.T)):
            eq_h = 0
            eq_p = 0            
            for e in range(len(self.T[i].q_c)):
                if self.schema[i][e] != None:
                    if self.T[i].q_c[e] == -1:
                        eq_h += 1
                    if self.T[i].q_c[e] == 1:
                        eq_p += 1

            t_config.append([eq_p, eq_h])
        ignored = []
        t_eqs = 0
        for i in range(len(t_config)):
            for e in range(len(t_config)):
                if i!= e and i not in ignored:

                    if t_config[i] == t_config[e] and len(self.T[i].q_c) == len(self.T[e].q_c):
                        ignored.append(e)
                        t_eqs += 1
        self.n_equivalent_Ts = t_eqs 
        
    def count_equivalent_l(self):
        #Count equivalent lines
        eqs = 0
        for i in range(len(self.T)):
            eq_h = 0
            eq_p = 0            
            for e in range(len(self.T[i].q_c)):
                if self.schema[i][e] != None:
                    if self.T[i].q_c[e] == -1:
                        eq_h += 1
                    if self.T[i].q_c[e] == 1:
                        eq_p += 1
            eqs += eq_h/2 + eq_p/2
        self.equiv_lines = eqs
                
                
        
    def count_loops(self):
        lp = 0
        for i in range(len(self.T)):
            lp += (len(self.T[i].q_c) - sum(self.T[i].q_c))/2
        self.loops = lp
    def count_holes(self):
        nh = 0
        for i in range(len(self.T)):
            for e in range(len(self.T[i].q_c)):
                if self.T[i].q_c[e] == -1:
                    nh += 1
        for i in range(len(self.H.q_c)):
            if self.H.q_c[i] == -1:
                nh += 1
        self.holes = nh
                
    def t_operators(self):
        tt = []
        for i in range(len(self.T_labels)):
            #tt.append([])
            n = len(self.T_labels[i])/2
            a = ""
            j = ""
            for e in range(n):
                a += self.T_labels[i][e] 
                j += self.T_labels[i][n+e]
            tt.append([n,a,j])
        self.t_ops = tt
    def v_operator(self):
        hh = []
        out = ""
        ins = ""
        for i in range(len(self.H.m_qc)):
            l = self.H.m_qc[i]
            if  l>0:
                l-=1

                #found a line passing through the interaction
                if self.H.q_a[l] == 1:
                    #Found a particle passing through
                    out += self.H_labels_c[i]
                    ins += self.H_labels_a[l]
                if self.H.q_a[l] == -1:
                    #Found a hole passing through
                    out += self.H_labels_c[i]
                    ins += self.H_labels_a[l]
        a_ignores = []            
        for i in range(len(self.H.m_qa)):
            l = self.H.m_qa[i]
            if  l<0 and l not in a_ignores:
                #found a line pair annihilating at the interaction
                l = abs(l) - 1
                out += self.H_labels_a[i]
                ins += self.H_labels_a[l]
                a_ignores.append(-(i+1))

        c_ignores = []
        for i in range(len(self.H.m_qc)):
            l = self.H.m_qc[i]
            if  l<0 and l not in c_ignores:
                #found a line pair created at the interaction
                l = abs(l) - 1
                out += self.H_labels_c[i]
                ins += self.H_labels_c[l]  
                c_ignores.append(-(i+1))                  
                                    
  
        self.h_op = [ins, out]
    def report(self):
        #print " Tlabels :", self.T_labels
        #print " Hlabelsa:", self.H_labels_a
        #print " Hlabelsc:", self.H_labels_c 
        #print "Equivavlent lines:", self.equiv_lines
        #print "Loops:", self.loops
        #print "Holes:", self.holes
        #print "Equivavlent Ts:",self.n_equivalent_Ts 
        print "Prefactor   :", self.prefactor        
        print "Permutations:", self.permutations
        print "Summations  :", self.summation
        print "Interaction :", self.h_op
        print "Cluster     :", self.t_ops
                    

class cpp_func():
    def __init__(self, schema, target):
        self.schema = schema
        self.build(target)
    def permute(self, c, s1, s2):     
        c2 = []
        for i in c:
            c2.append(i)
        #print c2
        s1_i = None
        s2_i = None

        for i in range(len(c)):
            if c[i] == s1:
                s1_i = i
                c2[s1_i] = s2
            if c[i] == s2:
                s2_i = i
                c2[s2_i] = s1
            
        c3 = ""
        for i in c2:
            c3 += i

        return c + "-(" + c3 +")"
        
    def build(self, target):
        prefactor = self.schema.prefactor
        permutations = self.schema.permutations
        sums = self.schema.summation
        interaction = self.schema.h_op
        cluster = self.schema.t_ops
        
        #Begin building the function
        s = "double %s = 0.0\n" %target
        c = ""
        #Core function
        c1 = interaction[0][0]
        c2 = interaction[1][0]
        for i in range(len(interaction[0])-1):
            c1 += ","+interaction[0][i+1]
            c2 += ","+interaction[1][i+1]     
        c += "%s(%s)(%s)" % ("v"+str(len(interaction[0])),c1, c2) #Two particle interaction
        for i in range(len(cluster)):
            c1 = cluster[i][1][0]
            c2 = cluster[i][2][0]
            for e in range(len(cluster[i][1])-1):
                c1 += ","+cluster[i][1][e+1]
                c2 += ","+cluster[i][2][e+1]
            c += "*%s(%s)(%s)" % ("t"+str(cluster[i][0]),c1, c2) #Two particle interaction
        
        #Permute all possible
        for i in range(len(permutations)):
                c = self.permute(c, permutations[i][0], permutations[i][1])
                
            
        
        c = "%s += " %target + c
        c+=";\n"
        sms = []
        for i in range(len(sums)):
            indent = 4 * (len(sums)-i-1) * " "
            s,e = self.for_loop(sums[i])
            sms.insert(0,indent + s)
            sms.append(indent + e)
        indent = 4*(len(sums)) * " " 
        sms.insert(len(sms)/2, indent+c)
        
        #merge function
        C=""
        for i in sms:
            C+= i
            
        C = "double %s = 0.0\n"%target+ C 

        C +="double %s *= %f\n" % (target, prefactor)
        #Set up summations
        #print "\n", C
        self.cppfunction = "\n" + C
        self.target = target
        #print sms
            
    
    def for_loop(self, p):
        if p in "abcdefgh":
            #return 
            lps = "for(int %s = nElectrons; %s < nStates; %s ++){\n" % (p,p,p)
            lpe= "}\n"
        if p in "ijklmnop":
            #return 
            lps = "for(int %s = 0; %s < nElectrons; %s ++){\n" % (p,p,p)
            lpe = "}\n"
        return lps, lpe  

class symbolic_diag():
    def __init__(self, schema, target):
        self.schema = schema
        self.target = target
        self.build(target)
    def build(self, target):
        prefactor = self.schema.prefactor
        permutations = self.schema.permutations
        sums = self.schema.summation
        interaction = self.schema.h_op
        cluster = self.schema.t_ops
        
        pre_fact, pre_denom = self.factorize(prefactor)
        s = ""
        for i in permutations:
            s += "P(%s)" % i
        s += "\\frac{%i}{%i}" % (pre_fact, pre_denom)
        s += " \\sum_{%s}" % sums
        s += " \\langle %s || %s \\rangle" % (interaction[0],interaction[1])
        for i in cluster:
            s += " t_{%s}^{%s}" % (i[2], i[1])
        self.stringform = s
    
    def factorize(self, prefactor):
        ret = prefactor
        fract = 1
        while not float(ret) == int(ret):
            ret*=2
            fract*=2
            
        return ret, fract
     
class feynman_diag():
    #a diagrammatic representation of the operators
    def __init__(self, schema, target, pos, vis = True):
        self.schema = schema
        self.target = target
        self.pos = pos
        self.vis = vis
        
        
        #Parameters for plotting
        self.shift = .2 #horizontal shift of interaction vertices
        self.scale = 1.0 #size of diagrams
        self.split = .25 #angle of paired leaving vertices
        
        #building diagrams
        self.Build(target)
    def Build(self,target):
        scale = self.scale
        prefactor = self.schema.prefactor
        permutations = self.schema.permutations
        sums = self.schema.summation
        interaction = self.schema.h_op
        cluster = self.schema.t_ops
        schema = self.schema.schema
        px = self.pos[0]
        py = self.pos[1]
        
        #Create all vertices in the cluster operators
        #Find total number of vertices:
        nt = 0
        for i in range(len(cluster)):
            nt += cluster[i][0]
        #print "NT:", nt
        self.t_nodes = []
        self.exit_nodes = []
        t0 = -scale*nt/2 
        self.NT_tot = nt
        print cluster
        for i in range(len(cluster)):
            self.t_nodes.append([])
            self.exit_nodes.append([])
            for e in range(cluster[i][0]):
                self.t_nodes[i].append(hnode([px+t0, py],  [1,-1],[]))
                print px+t0, py
                self.exit_nodes[i].append([])
                self.exit_nodes[i][-1].append(hnode([px+t0-self.split*scale, py+2*scale],[],[]))
                self.exit_nodes[i][-1].append(hnode([px+t0+self.split*scale, py+2*scale],[],[]))
                t0 += scale

        
        #Set up H-vertices
        #counting hole lines passing through the interaction
        n_h_h = []
        for i in range(10):
            if i in self.schema.H.m_qc:
                n_h_h.append(self.schema.H.q_a[i-1])
        n_passing_holes = n_h_h.count(-1)
        n_passing_parts = n_h_h.count(1)
        #print n_h_h

        #count particle lines passing 
        n_h_h = 0
        for i in range(10):
            n_h_h += self.schema.H.m_qc.count(i)
        
                        
        #counting the line pairs annihilating at the interaction
        n_h_a = 0
        for i in range(10):
            n_h_a += self.schema.H.m_qa.count(-i)

        #counting the line pairs created at the interaction
        n_h_c = 0
        for i in range(10):
            n_h_c += self.schema.H.m_qc.count(-i)
        
        n_h_a /= 2
        n_h_c /= 2
        
        #print "Passing holes:", n_passing_holes
        #print "Passing particles:" , n_passing_parts
        #print "Annihil:", n_h_a
        #print "Created:", n_h_c
        
        self.h_nodes = []
        h0 = - .5*(n_h_a + n_passing_holes +n_passing_parts + n_h_c)*scale +  self.shift
        for i in range(n_h_a):
            #Add vertices that annihilate a pair of lines
            self.h_nodes.append(hnode([px+h0,py+1*scale], [],[1,-1]))
            h0 += scale 
        for i in range(n_h_c):
            #Add vertices that create a pair of lines
            self.h_nodes.append(hnode([px+h0,py+1*scale], [1,-1],[]))
            h0 += scale
        for i in range(n_passing_parts):
            self.h_nodes.append(hnode([px+h0,py+1*scale], [1],[1]))
            h0 += scale            
        for i in range(n_passing_holes):
            self.h_nodes.append(hnode([px+h0,py+1*scale], [-1],[-1]))
            h0 += scale
        #Done setting up the H-vertices
        
        
        #Grouping up lines
        TL = self.schema.T_labels
        HLa = self.schema.H_labels_a
        HLc = self.schema.H_labels_c   

        lines_t = [] #list of all connections in the diagram, both external and internal
        lines_t_internal = [] #list of all connections in the diagram, both external and internal
        #for i in range(len(self.schema)):
        lines_h_out = []
        
        
        HLc = self.schema.H_labels_c 
        #import mapping 
        mqc = self.schema.H.m_qc
        mqa = self.schema.H.m_qa
        
        for i in range(len(schema)):            
            lines_t.append([])
            lines_t_internal.append([])            
            for e in range(len(schema[i])):
                if schema[i][e] == None:
                    lines_t[i].append([TL[i][e], [i,e], None])
                if schema[i][e] != None:
                    lines_t[i].append([TL[i][e], [i,e], i])
        for i in range(len(HLc)):
            lines_h_out.append([HLc[i], i, None])

        self.lines_h_out = lines_h_out
        self.lines_t = lines_t
        self.lines_t_int = lines_t_internal
        
        
        
        
        
        self.Plot(target)
    def Plot(self, target):
        if self.vis:
            figure(figsize = (2, self.NT_tot/1.5), dpi = 80, edgecolor = 'k',facecolor='white')
            axis('off')
            axes().set_aspect('equal', 'datalim')
        hold('on')
        t_nodes = self.t_nodes
        h_nodes = self.h_nodes
        exit_nodes = self.exit_nodes

        
        #draw all nodes and horizontal lines
        pos = self.pos
        scale = self.scale
        cluster = self.schema.t_ops
        lines_t = self.lines_t
        msize =5
        c = 'black'
        for i in range(len(self.t_nodes)):
            if len(self.t_nodes[i])==1:
                nc = self.t_nodes[i][0] #current node
                #Draw 
                plot( [nc.x - 0.25*scale,nc.x + 0.25*scale],[nc.y, nc.y], color = c)
                plot(nc.x, nc.y,"o", color = c, markersize = msize)
            if len(self.t_nodes[i])>1:
                nc_next = self.t_nodes[i][0]                
                plot(nc_next.x, nc_next.y,"o", color = c, markersize = msize)
                for e in range(len(t_nodes[i])-1):
                    nc = self.t_nodes[i][e]
                    nc_next = self.t_nodes[i][e+1]
                    
                    plot( [nc.x,nc_next.x],[nc.y, nc_next.y], color = c)
                    plot(nc_next.x, nc_next.y,"o", color = c, markersize = msize)
                    
        
        #Draw interaction vertices
        if len(self.h_nodes)==1:
            nc = self.h_nodes[0]
            plot(nc.x,nc.y,"o", color = c, markersize = msize)
            plot( [nc.x,nc.x + 0.5*scale],[nc.y,nc.y], color = c, ls = "dotted")
        if len(self.h_nodes)>1:
            #connect horizontal vertices
            nc = self.h_nodes[0]
            plot(nc.x,nc.y,"o", color = c, markersize = msize)
            for e in range(len(self.h_nodes)-1):
                nc = self.h_nodes[e]
                nc_next = self.h_nodes[e+1]
                plot( [nc.x,nc_next.x],[nc.y, nc_next.y], color = c,ls = "dotted")
                plot(nc_next.x, nc_next.y,"o", color = c, markersize = msize)
                
        #Draw all connected and unconnected lines
        for i in range(len(cluster)): 
            for u in range(cluster[i][0]):
                p1_ind = self.get_l_index(lines_t[i],cluster[i][1][u])
                h1_ind = self.get_l_index(lines_t[i],cluster[i][2][u])   
                
                #print lines_t[i][p1_ind] , lines_t[i][h1_ind]                 
                if lines_t[i][p1_ind][2] == lines_t[i][h1_ind][2] == None:
                    #print "Two external unconnected lines!"
                    ncon(t_nodes[i][u],exit_nodes[i][u][0],order = 0, p_h = -1)
                    ncon(t_nodes[i][u],exit_nodes[i][u][1],order = 0, p_h = 1)
                    #print lines_t[i][p1_ind], lines_t[i][h1_ind]           
                if lines_t[i][p1_ind][2] == lines_t[i][h1_ind][2] == 0 or  lines_t[i][p1_ind][2] == lines_t[i][h1_ind][2] == 1:
                    #print "Internal loop!"
                    #Scan for available loop
                    n = None
                    for t in range(len(h_nodes)):
                        if h_nodes[t].probe_loop_a():
                            #Found loop
                            h_nodes[t].connect_loop_a()
                            n = t
                            break
                    hn = h_nodes[n]
                    
                        
                    ncon(t_nodes[i][u],hn,order = -1, p_h = 1)
                    ncon(hn,t_nodes[i][u],order = -1, p_h = 1)
        for i in range(len(cluster)): 
            for u in range(cluster[i][0]):
                p1_ind = self.get_l_index(lines_t[i],cluster[i][1][u])
                h1_ind = self.get_l_index(lines_t[i],cluster[i][2][u])   
                
                #print lines_t[i][p1_ind] , lines_t[i][h1_ind]          
                if lines_t[i][p1_ind][2] != lines_t[i][h1_ind][2]:
                    #print "One line in , one line out"
                    #print lines_t[i][p1_ind], lines_t[i][h1_ind]
                    
                    #draw line lines_t[i][p1_ind]
                    #begin = [self.t_pos[i][u],self.h_pos[hn/2]]
                    #print "HPOS:          ", self.h_pos, hn/2, hn
                    #print "TPOS:          ", self.t_pos,i,u
                    #begin = [self.t_pos[i][u],self.h_pos[hn/2]] #X-vector
                    #plot(begin, "o")
                    
                    if lines_t[i][p1_ind][2] == None:
                        #print "    Draw external particle"
                        ncon(t_nodes[i][u],exit_nodes[i][u][0],order = 0, p_h = 1)
                        #end = [self.pos[1],self.pos[1]+2]
                        #begin = [self.t_pos[i][u],self.t_pos[i][u]] #X-vector
                        #self.draw_lines(1, begin,end, 1)
                        #hn += 1                        
                    if lines_t[i][p1_ind][2] != None:
                        #print "    Draw internal particle"
                        #Scan for available particle annihilator
                        n = None
                        for t in range(len(h_nodes)):
                            if h_nodes[t].probe_line_a(1):
                                #Found loop
                                h_nodes[t].connect_line_a(1)
                                n =t
                                break
                        hn = h_nodes[n]
                        
                        ncon(t_nodes[i][u],hn,order = 0, p_h = 1)
                        
                        
                        #begin = [self.t_pos[i][u],self.h_pos[hn/2]] #X-vector
                        #end = [self.pos[1],self.pos[1]+1]
                        #self.draw_lines(1, begin,end, 1)
                        #hn += 1
                    
                    if lines_t[i][h1_ind][2] == None:
                        #print "    Draw external hole"
                        ncon(t_nodes[i][u],exit_nodes[i][u][1],order = 0, p_h =- 1)
                        #begin = [self.t_pos[i][u],self.t_pos[i][u]] #X-vector
                        #end = [self.pos[1],self.pos[1]+2]
                        #self.draw_lines(1, begin,end, -1)
                        #hn += 1
                    if lines_t[i][h1_ind][2] != None:
                        #print "    Draw internal hole"
                        
                        n = None
                        for t in range(len(h_nodes)):
                            if h_nodes[t].probe_line_a(-1):
                                #Found loop
                                h_nodes[t].connect_line_a(-1)
                                n = t
                                break
                        hn = h_nodes[n]
                        ncon(t_nodes[i][u],hn,order = 0, p_h = -1)
                        
                        
                        #begin = [self.t_pos[i][u],self.h_pos[hn/2]] #X-vector
                        #end = [self.pos[1],self.pos[1]+1]
                        #self.draw_lines(1, begin,end, -1)
                        #hn += 1
                    #hn += 1
        
        #Finally, we need to draw all lines exiting from the interaction
        for i in range(len(h_nodes)):
            if h_nodes[i].probe_line_c(1):
                #create particle leaving the interaction
                hn = h_nodes[i]
                ncon(hn, hnode([hn.x - scale*0.2, 2*scale],[],[]), order = 0, p_h = 1)
            if h_nodes[i].probe_line_c(-1):
                #create particle leaving the interaction
                hn = h_nodes[i]
                ncon(hn, hnode([hn.x + scale*0.2, 2*scale],[],[]), order = 0, p_h = -1)
        if self.vis:
            show()
            hold('off')


    def get_centered_pair(self, L, excluded):
        #pops centered line pairs from list (particle/hole)
        i1,i2 = None, None
        for i in range(len(L)-1):
            if L[i][0] in "abcdefgh" and L[i+1][0] in "ijklmnop":
                i1 = L[i]
                i2 = L[i+1]
                del(L[i+1])
                del(L[i])
                break              
        return i1, i2
        

    def draw_lines(self, order, begin, end, ph = 0):
        n1 = position(begin[0],begin[1])
        n2 = position(end[0], end[1])
        if order == 1:
            plot(begin, end, color = "black")
            half_x = begin[0] + (begin[1]-begin[0])/1.5
            half_y = end[0] + (end[1]-end[0])/1.5
            
            orient_x =  (begin[1]-begin[0])/2.0
            orient_y =  (end[1]-end[0])/2.0
            if ph == 1:
                draw_arrow([half_x, half_y],[orient_x, orient_y])
            if ph == -1:
                draw_arrow([half_x, half_y],[-orient_x, -orient_y])

        if order %2 == 0:
            #draw pairs of lines
            for i in range(order/2):
                ncon(n1,n2,-1,1)
                ncon(n2,n1,-1,-1)
            print 0
        if order %2 != 0:
            print "draw straight line and pairs"
            
    def get_l_index(self, lines, target):
        ret = None
        for i in range(len(lines)):
            if lines[i][0] == target:
                ret = i
        return ret
    def get_t_index(self, To, target):
        ret = None
        for i in range(len(To)):
            for e in range(len(To[i])):
                if To[i][e] == target:
                    ret = [i,e]
        return ret[0], ret[1]
    
    def get_h_index(sefl,Ho,target):
        ret = None
        for i in range(len(Ho)):
            if Ho[i] == target:
                ret = i
        return ret

class position():
    def __init__(self, x, y):
        self.x = x
        self.y = y

class hnode():
    def __init__(self, pos, q_c, q_a):
        self.pos = pos
        self.x = pos[0]
        self.y = pos[1]
        self.q_c = q_c
        self.q_a = q_a
        self.setup()
    def setup(self):
        self.c_connections = [0,0]
        self.a_connections = [0,0]
        #particles, holes
        for i in self.q_c:
            if i == 1:
                self.c_connections[0] = None
            if i == -1:
                self.c_connections[1] = None
        for i in self.q_a:
            if i == 1:
                self.a_connections[0] = None
            if i == -1:
                self.a_connections[1] = None
    def probe_loop_a(self):
        #Probe annihilation for available loop connection
        ret = False
        if self.a_connections[0] == self.a_connections[1] == None:
            ret = True #A connection is possible
        return ret
    def probe_line_a(self, linetype):
        #Probe annihilation for 
        ret = False
        if linetype == -1:
            #scan for holw
            if self.a_connections[1] == None:
                ret = True
        if linetype == 1:
            #scan for particle
            if self.a_connections[0] == None:
                ret = True
        return ret
    def probe_line_c(self, linetype):
        #Probe annihilation for 
        ret = False
        if linetype == -1:
            #scan for holw
            if self.c_connections[1] == None:
                ret = True
        if linetype == 1:
            #scan for particle
            if self.c_connections[0] == None:
                ret = True
        return ret               

    def connect_loop_a(self):
        if self.probe_loop_a():
            self.a_connections[0] = 1
            self.a_connections[1] = 1
    def connect_line_a(self, linetype):
        if self.probe_line_a(linetype):
            if linetype == 1:
                self.a_connections[0] = 1
            if linetype == -1:
                self.a_connections[1] = 1
                
        
            
            
        
        
                
        

def nconnect(n1,n2,S,order="l0", p_h = None):
    N = 60

    Phx = (n1.x+n2.x)/2.0
    Phy = (n1.y+n2.y)/2.0
    
    lP = sqrt((n2.x-n1.x)**2 + (n2.y-n1.y)**2)
    dPx = (n2.x-n1.x)/lP    
    dPy = (n2.y-n1.y)/lP
    
    Cx = Phx - S*dPy
    Cy = Phy + S*dPx  
    
    
    lC = sqrt((S*dPy)**2 + (S*dPx)**2)
    #node(Phx,Phy, c="blue")
    #node(Cx,Cy, c="red")
    R = sqrt((Cx-n1.x)**2 + (Cy-n1.y)**2)

    lPC0 = sqrt(((Cx+R)-n1.x)**2 + (Cy - n1.y)**2)
    lPC1 = sqrt(((Cx+R)-n2.x)**2 + (Cy - n2.y)**2)

    dalpha = arccos((2*R**2 - lP**2)/(2.0*R**2))

    
    CPx = n1.x - Cx
    CPy = n1.y - Cy
    X,Y = 0,0
    if order == "0":
        #X = [n1.x, n2.x]
        #Y = [n1.y, n2.y]
        X = linspace(n1.x, n2.x, 6) #need to do this to get arrows working
        Y = linspace(n1.y, n2.y, 6)
    if order == "l0":
        if S<0:
            dalpha = 2*pi - dalpha
        A = linspace(0,dalpha, N)
        X,Y = rotate_v(CPx,CPy,A)
        X+=Cx
        Y+=Cy
    if order == "r0":
        if S>0:
            dalpha = 2*pi - dalpha      
        A = linspace(0,-dalpha, N)
        X,Y = rotate_v(CPx,CPy,A)
        X+=Cx
        Y+=Cy
        
    msize = 10
    if p_h == 1:
        draw_arrow([X[len(X)/2],Y[len(X)/2]], [-dPx,-dPy])
        #X[len(X)/2],Y[len(X)/2]
        #plot(X[len(X)/2],Y[len(X)/2], "^", color = "black", markersize = msize)
    if p_h == -1:
        draw_arrow([X[len(X)/2],Y[len(X)/2]], [dPx,dPy])
        #plot(X[len(X)/2],Y[len(X)/2], "v", color = "black", markersize = msize)
    plot(X,Y, color = "black")
    
def ncon(n1,n2,order = 0, p_h = None):
    if order == 0:
        nconnect(n1,n2,1,"0", p_h)
    if order > 0:
        nconnect(n1,n2,(-2+order),"l0", p_h)
    if order < 0:
        nconnect(n1,n2,(-2-order),"r0", p_h)

def draw_arrow(pos, point, s = .2, h = .1):
    #Draw an arrow at pos, pointing in the direction of point.
    #normalize direction
    p2 = sqrt(point[0]**2 + point[1]**2)
    point[0] /= p2
    point[1] /= p2
    
    #pi/2 degree rotation
    p_rotx, p_roty = point[1], -point[0]

    x0, y0 = pos[0], pos[1]
    x1, y1 = pos[0] - s*point[0], pos[1] - s*point[1]
    
    #plot the arrow
    plot([x0, x1+h*p_rotx],[y0, y1+h*p_roty], color = "black")
    plot([x0, x1-h*p_rotx],[y0, y1-h*p_roty], color = "black")
                
def rotate_v(x,y,alpha):
    ca = cos(alpha)
    sa = sin(alpha)
    return ca*x - sa*y, sa*x + ca*y            
    
        
def normal_ordered_hamiltonian():
    #Two particle hamiltonian
    F1 = Operator([1],[1])
    F2 = Operator([-1],[-1])
    F3 = Operator([],[1,-1])
    F4 = Operator([1,-1],[])
    
    V1 = Operator([1,1],[1,1])
    V2 = Operator([-1,-1],[-1,-1])
    V3 = Operator([1,-1],[1,-1])
    
    V4 = Operator([1],[1,1,-1])
    V5 = Operator([1,1,-1],[1])
    
    V7 = Operator([1,-1,-1],[-1])
    V6 = Operator([-1],[1,-1,-1])
    V9 = Operator([1,1,-1,-1],[])
    V8 = Operator([],[1,1,-1,-1])
    return [F1,F2,F3,F4,V1,V2,V3,V4,V5,V6,V7,V8,V9]                     
        
def cluster_operator(configuration):
    #configuration =[]
    T1= Operator([],[1,-1])
    T2= Operator([],[1,1,-1,-1])
    T3= Operator([],[1,1,1,-1,-1,-1])
    T4= Operator([],[1,1,1,1,-1,-1,-1,-1])
    T5= Operator([],[1,1,1,1,1,-1,-1,-1,-1,-1])
    T_list = [T1, T2, T3, T4, T5]
    ret = []
    for i in configuration:
        ret.append(T_list[i-1])
    return ret
                                    

class O():
    def __init__(self, H, T, name = "CC"):
        self.name = name
        self.H = H
        self.T = T
        self.diagrams = contract(H, T)
        self.Ndiag = len(self.diagrams)
        self.schematics = []
        self.stringforms = []
        self.cppcodes = []
        subindex = ["a","b","c","d","e","f","g","h","i","j","k","l","m","n","o","p","q","r","s","t","u","v","x","y","z"]
        for i in self.diagrams:
            self.schematics.append(schematic(i))
            self.schematics[-1].setup()
            self.N = len(self.schematics)
        for e in range(len(self.schematics)):
            
            c = self.schematics[e]
            #name = self.name+"%i" %e
            sym = symbolic_diag(c,"%s" % (self.name+subindex[e]))
            self.stringforms.append(sym.stringform)
            CPPc = cpp_func(c,"%s" % (self.name+subindex[e]))
            self.cppcodes.append(CPPc.cppfunction)
    def E(self, i = None):
        E = self.H.E
        for i in range(len(self.T)):
            E += self.T[i].E
        return E
            
    def latex(self, i):
        ret = None
        try:
            ret =  self.stringforms[i]
        except:
            print "Invalid index given:", i
            print "There are only %i diagrams in this contraction." % self.N
        return ret
    def code(self, i):
        try:
            print self.cppcodes[i]
        except:
            print "Invalid index given:", i
            print "There are only %i diagrams in this contraction." % self.N
    def diagram(self,i, pos = [0,0], vis = True):
        try:
            dia = feynman_diag(self.schematics[i], "D9", pos, vis)
        except:
            print "Invalid index given:", i
            print "There are only %i diagrams in this contraction." % self.N

def list_mult(A, B):
    #multiply two lists
    ret = []
    for i in A:
        for e in B:
            ret.append(i+e)
    return ret


def list_pow(T, n):
    #Return list to the power of n
    th = T
    g = []
    for i in range(n-1):
        th = list_mult(th,T)
        g=th
    if n==1:
        g = T
    return g
            
def factorial(n):
    N = 1
    for i in range(1,n+1):
        N*=i
    return N

def expand_ansatz(T, n):
    #Taylor expand the exponential ansazt (exp(list(t-operatrs))) to a given truncation (n+1)
    expanded_eT = []
    prefactors  = []
    for i in range(1,n+1): #ignore i=0 -> uncorrelated energy
        elem = []
        prefactors.append(factorial(n)) # (**-1) inverted
        expT = list_pow(T,i)
        expanded_eT+=expT
    return expanded_eT

def combine_to_excitation(H,T,E):
    #Find combinations of operator elements in H, T that produce excitation level i
    for h in range(len(H)):
        for i in range(len(expT)):
            #print h, i
            contr = O(H[h],expT[i])
            if contr.E() == E:
                for e in range(contr.Ndiag):
                    #contr.latex(e)
                    contr.latex(e)
                    #print "Excitation:", contr.E() #Excitation
                    #contr.diagram(e, [4*h,4*i + 4*e], False)
                    #show()

In [34]:
H = Operator([1,-1],[])
T = [Operator([],[1,1,-1,-1]),Operator([],[1,-1])]

contr = O(H,T)
print "Excitation:", contr.E()
Math(contr.latex(1))
contr.diagram(1)


Excitation: 2.0
[[2, 'ab', 'ki'], [1, 'c', 'j']]
-1.5 0
-0.5 0
0.5 0

In [16]:
contr.code(1)


double CCb = 0.0
for(int c = nElectrons; c < nStates; c ++){
    for(int k = 0; k < nElectrons; k ++){
        CCb += v1(k)(c)*t2(a,b)(k,i)*t1(c)(j)-(v1(k)(c)*t2(b,a)(k,i)*t1(c)(j))-(v1(k)(c)*t2(a,b)(k,j)*t1(c)(i)-(v1(k)(c)*t2(b,a)(k,j)*t1(c)(i)));
    }
}
double CCb *= 1.000000


In [ ]:
contr.diagram(1)

In [4]:
H = Operator([1,1,-1],[1])
T = [Operator([],[1,1,-1,-1]),Operator([],[1,-1])]
contr = O(H,T)
contr.diagram(0)
contr.diagram(1)



In [6]:
print contr.code(1)
print contr.code(2)


double D9 = 0.0
for(int d = nElectrons; d < nStates; d ++){
    for(int k = 0; k < nElectrons; k ++){
        for(int c = nElectrons; c < nStates; c ++){
            D9 += v2(c,k)(b,d)*t2(c,a)(k,i)*t1(d)(j)-(v2(c,k)(a,d)*t2(c,b)(k,i)*t1(d)(j))-(v2(c,k)(b,d)*t2(c,a)(k,j)*t1(d)(i)-(v2(c,k)(a,d)*t2(c,b)(k,j)*t1(d)(i)));
        }
    }
}
double D9 *= 1.000000


double D9 = 0.0
for(int k = 0; k < nElectrons; k ++){
    for(int d = nElectrons; d < nStates; d ++){
        for(int c = nElectrons; c < nStates; c ++){
            D9 += v2(c,k)(b,d)*t2(c,a)(i,j)*t1(d)(k)-(v2(c,k)(a,d)*t2(c,b)(i,j)*t1(d)(k))-(v2(c,k)(b,d)*t2(c,a)(j,i)*t1(d)(k)-(v2(c,k)(a,d)*t2(c,b)(j,i)*t1(d)(k)));
        }
    }
}
double D9 *= 1.000000


In [ ]:
#find energy
H = Operator([1,1,-1,-1],[])
T = [Operator([],[1,1,-1,-1])]
contr1 = O(H,T)
contr1.diagram(0)
contr1.latex(0)

In [25]:
#find energy
H = Operator([1,1,-1,-1],[])
T = [Operator([],[1,-1]),Operator([],[1,-1])]
contr2 = O(H,T)
contr2.diagram(0)
contr2.latex(0)


Out[25]:
'\\frac{1}{2} \\sum_{ckdl} \\langle kl || cd \\rangle t_{k}^{c} t_{l}^{d}'

In [26]:
#find energy
H = Operator([1,-1],[])
T = [Operator([],[1,-1])]
contr3 = O(H,T)
contr3.diagram(0)
contr3.latex(0)


Out[26]:
'\\frac{1}{1} \\sum_{ck} \\langle k || c \\rangle t_{k}^{c}'

In [29]:
print contr1.code(0)
print contr2.code(0)
print contr3.code(0)


double D9 = 0.0
for(int l = 0; l < nElectrons; l ++){
    for(int k = 0; k < nElectrons; k ++){
        for(int d = nElectrons; d < nStates; d ++){
            for(int c = nElectrons; c < nStates; c ++){
                D9 += v2(k,l)(c,d)*t2(c,d)(k,l);
            }
        }
    }
}
double D9 *= 0.250000


double D9 = 0.0
for(int l = 0; l < nElectrons; l ++){
    for(int d = nElectrons; d < nStates; d ++){
        for(int k = 0; k < nElectrons; k ++){
            for(int c = nElectrons; c < nStates; c ++){
                D9 += v2(k,l)(c,d)*t1(c)(k)*t1(d)(l);
            }
        }
    }
}
double D9 *= 0.500000


double D9 = 0.0
for(int k = 0; k < nElectrons; k ++){
    for(int c = nElectrons; c < nStates; c ++){
        D9 += v1(k)(c)*t1(c)(k);
    }
}
double D9 *= 1.000000


In [2]:
T = [[Operator([],[1,1,-1,-1])]] #CCD
expT = expand_ansatz(T,2)
H = normal_ordered_hamiltonian()
combine_to_excitation(H,expT, 2)


P(ab)P(ij)\fract{1}{1} \sum_{c} \langle c || b \rangle t_{ij}^{ca}
P(ab)P(ij)\fract{-1}{1} \sum_{k} \langle k || j \rangle t_{ki}^{ab}
P(ij)\fract{1}{2} \sum_{cd} \langle cd || ab \rangle t_{ij}^{cd}
P(ab)\fract{1}{2} \sum_{kl} \langle kl || ij \rangle t_{kl}^{ab}
P(ab)P(ij)\fract{-1}{1} \sum_{ck} \langle ck || bj \rangle t_{ki}^{ca}
P(ab)\fract{1}{2} \sum_{cdkl} \langle kl || cd \rangle t_{ki}^{cd} t_{lj}^{ab}
P(ij)P(ab)\fract{1}{4} \sum_{cdkl} \langle kl || cd \rangle t_{ij}^{cd} t_{kl}^{ab}
P(ij)\fract{1}{2} \sum_{ckld} \langle kl || cd \rangle t_{kl}^{ca} t_{ij}^{db}
\fract{1}{2} \sum_{ckdl} \langle kl || cd \rangle t_{ki}^{ca} t_{lj}^{db}

In [ ]: