In [45]:
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):
#print "Connection " , diagram[2]
#print "Number ", diagram[3]
self.H = diagram[0]
self.T = diagram[1]
self.schema = diagram[2]
self.i = diagram[3]
#later editing, remove this if faulty
self.convert_schema()
def convert_schema(self):
schema = self.schema
for i in schema:
n = len(i)
sc = i[n/2:n]
sc.reverse()
i[n/2:n] = sc
def scan_extract(self,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'] #Target labels
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']
#print "Schema:", self.schema
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 and internal lines
for i in range(len(self.schema)):
n = len(self.schema[i])
for e in range(len(self.schema[i])):
ne = n-e-1
de = e
if self.schema[i][de] == None:
if self.T[i].q_c[de] == 1:
self.T_labels[i][de] = plabels_t.pop()
if self.T[i].q_c[de] == -1:
self.T_labels[i][de] = hlabels_t.pop()
"""
if self.schema[i][ne] == None:
if self.T[i].q_c[ne] == 1:
self.T_labels[i][ne] = plabels_t.pop()
if self.T[i].q_c[ne] == -1:
self.T_labels[i][ne] = hlabels_t.pop()
"""
if self.schema[i][de] != None:
h_index = self.schema[i][de]
if self.T[i].q_c[de] == 1:
particlelabel = plabels.pop()
self.T_labels[i][de] = particlelabel
self.H_labels_a[h_index] = particlelabel
#print self.H.m_qa
if self.H.m_qa[h_index] >= 0:
self.H_labels_c[abs(self.H.m_qa[h_index])-1] = plabels_t.pop()
if self.T[i].q_c[de] == -1:
holelabel = hlabels.pop()
self.T_labels[i][de] = holelabel
self.H_labels_a[h_index] = holelabel
if self.H.m_qa[h_index] >= 0:
self.H_labels_c[abs(self.H.m_qa[h_index])-1] = hlabels_t.pop()
"""
if self.schema[i][ne] != None:
h_index = self.schema[i][ne]
if self.T[i].q_c[ne] == 1:
particlelabel = plabels.pop()
self.T_labels[i][ne] = particlelabel
self.H_labels_a[h_index] = particlelabel
self.H_labels_c[abs(self.H.m_qa[h_index])-1] = plabels_t.pop()
if self.T[i].q_c[ne] == -1:
holelabel = hlabels.pop()
self.T_labels[i][ne] = holelabel
self.H_labels_a[h_index] = holelabel
self.H_labels_c[abs(self.H.m_qa[h_index])-1] = hlabels_t.pop()
"""
"""
#label all external and internal lines
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()
if self.schema[i][e] != None:
h_index = self.schema[i][e]
if self.T[i].q_c[e] == 1:
particlelabel = plabels.pop()
self.T_labels[i][e] = particlelabel
self.H_labels_a[h_index] = particlelabel
self.H_labels_c[abs(self.H.m_qa[h_index])-1] = plabels_t.pop()
if self.T[i].q_c[e] == -1:
holelabel = hlabels.pop()
self.T_labels[i][e] = holelabel
self.H_labels_a[h_index] = holelabel
self.H_labels_c[abs(self.H.m_qa[h_index])-1] = hlabels_t.pop()
"""
#print "---"
#print self.H_labels_c
#print self.H_labels_a
#print self.T_labels
#print "---"
"""
#Label all external lines exiting the cluster operators (targets)
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()
"""
#Ensure all external interaction lines labeled
for i in range(len(self.H_labels_c)):
if self.H_labels_c[i] == None:
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.find_permutations()
self.calc_prefactor()
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(i)])
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]: #equal orientation
if line1[2] != line2[2]: #inequal origin #AND origins differ
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))
self.prefactor *= (-1)**len(self.permutations)
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
#count quasi-loops ? #Important bug fix, 2.nov 2014
#for i in self.schema:
# lp += i.count(None)/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 = []
#print self.T_labels
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]
#a += self.T_labels[i][n-(e+1)]
#j += self.T_labels[i][n+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
ins += self.H_labels_c[i]
out += self.H_labels_a[l]
if self.H.q_a[l] == -1:
#Found a hole passing through
ins += self.H_labels_c[i]
out += self.H_labels_a[l]
######
######
# This needs to happen in the correct ordeR!
######
######
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]
v_out = []
v_ins = []
for i in range(len(self.schema)):
for e in range(len(self.schema[i])):
if self.schema[i][e] != None:
h_index = self.schema[i][e]
if self.T[i].q_c[e] == 1:
#particle in
v_ins.append(self.T_labels[i][e])
coupled = self.H.m_qa[h_index]
"""
if coupled < 0:
#reflecting at interaction
nc = abs(coupled) - 1
v_out.append(self.H_labels_a[nc])
"""
if coupled > 0:
#passing through interaction
nc = abs(coupled) - 1
#print coupled, self.H_labels_c, self.H.m_qa
v_out.append(self.H_labels_c[nc])
if self.T[i].q_c[e] == -1:
#hole out
v_out.append(self.T_labels[i][e])
coupled = self.H.m_qa[h_index]
"""
if coupled < 0:
#reflecting at interaction
nc = abs(coupled) - 1
v_ins.append(self.H_labels_a[nc])
"""
if coupled > 0:
#passing through interaction
nc = abs(coupled) - 1
v_ins.append(self.H_labels_c[nc])
#print "New config:", v_out, v_ins
out = ""
ins = ""
for i in range(len(v_out)):
out += v_out[i]
ins += v_ins[i]
self.h_op = [out, ins]
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
s = "prefactor: %f" % self.prefactor + "\n"
s += "permutations: %s" % self.permutations+ "\n"
s += "sum: %s" % self.summation+ "\n"
s += "h_ops: %s" % self.h_op+ "\n"
s += "t_ops: %s" % self.t_ops+ "\n"
return s
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)" % ("vmin"+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)" % ("tf"+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 +="%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'
plot(self.t_nodes[0][0].x-.01,self.t_nodes[0][0].y-.01,".", color ="White")
plot(self.t_nodes[0][0].x-.01,self.t_nodes[0][0].y+1.1,".", color ="White")
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] != None:
#if lines_t[i][p1_ind][2] == lines_t[i][h1_ind][2] == 0:
#print "Internal loop!"
#Scan for available loop
n = None
for t in range(len(h_nodes)):
#print "t:",t
if h_nodes[t].probe_loop_a():
#print "Found loop!"
#Found loop
h_nodes[t].connect_loop_a()
n = t
break
try:
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)
#Remove if causing trouble
except:
#connect to separate vertices
#connect particle line
for t in range(len(h_nodes)):
if h_nodes[t].probe_line_a(-1):
#print "Found particle!"
h_nodes[t].connect_loop_a()
n = t
break
hn = h_nodes[n]
ncon(t_nodes[i][u],hn,order = 0, p_h = -1)
#connect particle line
for t in range(len(h_nodes)):
if h_nodes[t].probe_line_a(1):
#print "Found particle!"
h_nodes[t].connect_loop_a()
n = t
break
hn = h_nodes[n]
ncon(t_nodes[i][u],hn,order = 0, p_h = 1)
#print "Could not connect loop."
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, nc.y+1*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, nc.y+1*scale],[],[]), order = 0, p_h = -1)
if self.vis:
hold('off')
show()
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"
print ""
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_old(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
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 = []
self.reports = []
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]
self.reports.append(c.report())
#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:
ret = self.cppcodes[i]
except:
print "Invalid index given:", i
print "There are only %i diagrams in this contraction." % self.N
return ret
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 report(self, i):
ret = None
try:
ret = self.reports[i]
except:
print "Invalid index given:", i
print "There are only %i reports in this contraction." % self.N
return ret
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, config = [1,0,0,0]):
#Find combinations of operator elements in H, T that produce excitation level i
skew = 0
tex = []
for h in range(len(H)):
shift = 0
for i in range(len(T)):
#print h, i
#print H
target_e = H[h].excitation()
for t in T[i]:
target_e += t.excitation()
if target_e == E:
#print h,i
stringname = "CC%i_%i" % (h,i)
contr = O(H[h],T[i], name = stringname)
if contr.E() == E:
for e in range(contr.Ndiag):
if config[0] == 1:
tx = contr.latex(e)
#Math(tx)
#print tx
tex.append(tx)
if config[1] == 1:
contr.diagram(e, [4*h,shift], True)
if config[1] == 2:
contr.diagram(e, [4*h,shift], False)
if config[3] == 1:
print contr.code(e)
if config[2] == 1:
print contr.reports[i]
shift += 4
if config[1] == 2:
axes().set_aspect('equal', 'datalim')
show()
return tex
def combine_all(H,T, config = [1,0,0,0]):
#Find combinations of operator elements in H, T that produce excitation level i
tex = []
for h in range(len(H)):
shift = 0
for i in range(len(expT)):
#print h, i
contr = O(H[h],expT[i])
for e in range(contr.Ndiag):
if config[0] == 1:
tx = contr.latex(e)
#print tx
tex.append(tx)
if config[0] == 2:
tx = contr.latex(e)
tex.append(tx)
if config[1] == 1:
contr.diagram(e, [4*h,shift], True)
if config[1] == 2:
contr.diagram(e, [4*h,shift], False)
if config[3] == 1:
contr.code(e)
shift+= 4
del(contr)
if config[1] ==2:
axes().set_aspect('equal', 'datalim')
show()
return tex
"Because of the often large number of terms in all versions of the CCM it is rather hard work to obtain the explicit equations by hand. And after this is done one has to write a program to put them into the computer." - H. G. Kümmel
This program enables the derivation of the coupled cluster equations by the use of computer.
The code in the cell above and this notebook is intended to provide an environment for deriving CC equations, complementary code and diagrams using the diagrammatic approach outlined in Shavitt and Bartlett (S-B) (p. 292). The code is mainly written for educational purposes, and has not been extensively tested, optimized or cleaned.
The core operation in CCAlgebra is the combining of operators. In CC calculations one often employs a diagrammatic approach instead of Wicks theorem to evaluate contractions of second quantized creation/annihilation operators, and this code is an attempt at working the diagrammatic rules into a practical code. All operators defined in this code is assumed to be normal ordered - this is explained in more detail below.
In the following examples, we will explore the standard functionality of CCAlgebra.
The operators we work with here are the cluster operators (T_1, T_2 ...) and the hamiltonian operators. Normal ordering of the hamiltonian for a system of interacting particles will typically produce a number of parts, each corresponding to a diagrammatic filament. This is discussed in more detail in S-B (p. 112, 4.29). (By filament, I mean a single, unconnected operator, that in itself would not produce any contribution when evaluated on the vacuum state.)
In CCalgebra operators are defined as shown below. The hamiltonian filament given here corresponds to the last diagram on page 180 in S-B, (9.105).
In [2]:
#How to define operators
H = Operator([1,-1],[]) # The first list defines q-particle annihilations (below interaction),
# while the second is creations. Particles and holes are indicated by 1 and -1 respectively.
T_1 = Operator([],[1,-1]) #The T_1 cluster operator
T_2 = Operator([],[1,1,-1,-1]) #The T_2 operator; all lists must be normal ordered
Note that the q-particle operators are given in normal order - q-particle creation operators to the left.
Operators may be diagrammatically connected by applying the rules laid out in S-B (p.298-299). The implementation in CCAlgebra uses the permutations library to produce all possible ways of connecting the diagrams, whereby it sorts out the distinct diagrams.
The basic rules are that all lines below the hamiltonian filament must connect to the one of the cluster operators, and all cluster operators must be connected to the hamiltonian filament. When this connection of operators acts on the vacuum state, it may produce an excitation or it may leave the vacuum state at the same energy level.
A class object O(O1, [O2, O3, ...]) is defined and will contain all distinct ways of connecting the operators. (The name should probably be changed to something more explanatory such as "connections()". An example is given by connecting the previously defined operators:
In [3]:
contraction = O(H,[T_1]) #The second argument must be a list of operators (T_1 * T_1 = [T_1, T_1])
A connection of opererators will in principle produce a new set of operators, which share a common property understood as the excitation level of the operators. Any operator in second quantization have this property, and it is easily evaluated by counting the number of lines exiting the operator and subtracting the the number of lines entering the operator. The excitation level is this number divided by 2, as the interaction should leave the number of particles constant.
Considering the connection made above, we may evaluate the excitation level as well as the number of distinct diagrams produced in the connection:
In [4]:
print "Excitation:", contraction.E() #Excitation energy of connected operators
print "Number of distinct diagrams:", contraction.N
In [5]:
c_tex = contraction.latex(0) #The first of possibly multiple distinct diagrams in latex format.
Math(c_tex)
Out[5]:
Compare this sum to the middle upper diagram on page 297 in S-B. As the final excitation of this diagram is 0, it will contribute to the energy equation. Note that while S-B makes an exception in the labelling of internal lines in the energy diagrams, CCAlgebra does not make this distinction. (S-B reserves the labels a,b,i,j for target indices in general, except for the energy contribution.)
In [6]:
contraction.diagram(0, [0,0], True) #Display diagram 0 at position x/y = [0,0] using internal show() function (matplotlib)
The diagram above is identical to the one in S-B, page 297. The dotted line represents the the interaction (H), while the horizontal straight line denotes the $T_1$ cluster operator. All lines are connected in the diagram above.
In general, CCAlgebra should produce diagrams indistinct to those found in S-B. In other words, some discrepancies in features might occur when comparing, but these discrepancies will not affect the outcome of any CC calculation.
There are two different options for displaying diagrams. The last parameter in the function call above is set to "True" to indicate that the diagram is to be both plotted and "show()"n by the O() class. The list containing two zeros is the position of the diagram in the x,y plane. The code is written this way to enable plotting of multiple diagrams in the same figure. (See the example in section 3.5 on how to use this in other ways).
In [7]:
print contraction.code(0)
This code may be implemented in a class in any CC-implementation in C++, but note that it is in itself not a complete class. This functionality is under developement, and will be implemented later on.
In [8]:
H = Operator([1,1,-1],[1])
contraction = O(H,[T_1, T_2]) #The second argument must be a list of operators (T_1 * T_1 = [T_1, T_1])c_tex = contraction.latex(0) #The first of possibly multiple distinct diagrams in latex format.
c_tex = contraction.latex(0) #The first of possibly multiple distinct diagrams in latex format.
Math(c_tex)
Out[8]:
In [9]:
plot([-2,7],[-2,4.5], color = "white") #just making the diagrams scale correctly
contraction.diagram(0, [0,0], False)
contraction.diagram(1, [3,0], False)
contraction.diagram(2, [6,0], False)
axis("off")
show()
In [10]:
print contraction.code(0)
While the examples above deal with single filaments of the hamiltionian, the full normal ordered hamiltionian will be a sum of more than one such filaments. When setting up the CC equations we will need to consider the connection between all this filaments and all terms in the expansion of the exponential ansatz.
To handle this, we will employ some complementary functions below.
The function expand_ansatz(list, int) takes as input a list with each element a list-embraced operator - typically cluster operators in the cluster expansion. The integer value indicates at what order the expansion should be truncated. Note that (1) very high values of int will most likely cause a system overload, and (2) when using a hamiltonian with at most two-body interactions, no higher orders than 4 is needed due to the fact that the hamiltionian needs to connect to all cluster-operators to its right.
In the example below it is shown that any list of lists may be used as input for the expand_ansatz() function:
In [11]:
print expand_ansatz([["a"],["b"]], 3)
This function is called without any arguments and returns a prewritten list of operators corresponding to the normal ordered hamiltonian in S-B. (p.280-281). It may be illuminating to see how it is set up while comparing each element to its diagrammatic counterpart:
In [12]:
def normal_ordered_hamiltonian():
#These elements corresponds to (9.105) in S-B:
F1 = Operator([1],[1]) #Excitation level: 0
F2 = Operator([-1],[-1]) #E:0
F3 = Operator([],[1,-1]) #E:+1
F4 = Operator([1,-1],[]) #E:-1
#These elements correspond to (9.107) in S-B:
V1 = Operator([1,1],[1,1]) #E:0
V2 = Operator([-1,-1],[-1,-1]) #E:0
V3 = Operator([1,-1],[1,-1]) #E:0
V4 = Operator([1],[1,1,-1]) #E:+1
V5 = Operator([1,1,-1],[1]) #E:-1
V7 = Operator([1,-1,-1],[-1]) #E:-1
V6 = Operator([-1],[1,-1,-1]) #E:+1
V9 = Operator([1,1,-1,-1],[]) #E:-2
V8 = Operator([],[1,1,-1,-1]) #E:+2
return [F1,F2,F3,F4,V1,V2,V3,V4,V5,V6,V7,V8,V9]
In [13]:
print cluster_operator([1,2,3]) # returns [T_1, T_2, T_3]
So far we have introduced all the functionality needed to set up the elements in the equations, but we have not set up the equations themselves. To do this, we need yet another pair of functions.
When deriving the equations we need to consider all possible connections between hamiltonian filaments and combinations of cluster operators in the cluster expansion. Currently there are two such higher level functions available.
Combining all elements in the normal ordered hamiltonian with all elements in the cluster expansion will produce a high number of excited slater determinants. To see all of these we may use the combine_all() function. In the following example we set up a normal ordered hamiltonian and a list of list of cluster operators, and then we find all combinations between these operators, representing them mathematically.
In [14]:
H = normal_ordered_hamiltonian()
T_1 = Operator([],[1,-1]) #The T_1 cluster operator
T_2 = Operator([],[1,1,-1,-1]) #The T_2 operator; all lists must be normal ordered
expT = expand_ansatz([[T_1],[T_2]], 1)
tx = combine_all(H,expT, [1,0,0,0]) #The list is explained below
S = "0 = "
for i in tx:
S+= i
Math(S)
Out[14]:
The combine_all() function is called with the two listst of operators as parameters and a maybe cryptic list with instructions on how to represent the combinations. The list of instructions works as follows:
[0,0,0,0] - Produce no representations, only iterate through all connections
[1,0,0,0] - Return a string of latex-formatted expression for each connection.
[0,1,0,0] - Plot all diagrams in same figure.
[0,2,0,0] - Plot each diagram in a separate figure. (May produce a lot of diagrams.)
[0,0,1,0] - Display the report for each diagram. (Explained in section 3 above)
[0,0,0,1] - Print out the code for all diagrams.
The arguments may ofcourse be combined if more than one index is set to non-zero values:
[1,2,1,1] - Return string, plot separate diagrams, show report and print code.
Important: This function is not very practical, and will most likely cause a memory overload if applied to diagrams with many permutations in the sums (due to the length of the C++ string.). For all practical purposes, the combine_to_excitation() function should be applied.
While the combine_all() function might have som educational value, it will show a lot more diagrams than what we in reality seek. When setting up the energy- and amplitide equations, we seek only contributions that correspond to certain excitation levels. This is beause the application of the exponential ansatz connected to the hamiltonian will produce a variety of excited states using the virtual orbitals in the SD, and when solving for the amplitudes we need to project these states down on a bra SD of the target excited state. (this formulation is a bit messy)
To find only the excitation levels of interest, i.e. the contributions to each equation, we may write (using the same functions as above):
In [15]:
tx = combine_to_excitation(H,expT,1, [1,0,0,0]) #All combinations that produce excitation level 1
S = "0 ="
for i in tx:
S+= "+"+ i
Math(S)
Out[15]:
Note the third argument, setting the excitation level of interest to 1. All the sums above have this excitation level. We may also display these 6 contributions to the $t_1$ amplitude equation as diagrams:
In [16]:
tx = combine_to_excitation(H,expT,1, [0,1,0,0]) #All combinations that produce excitation level 1
We are now ready to derive the full CCD equations. We perform the operations in two separate cells, one for the energy and another for the $t_2$ amplitude equation. This example might be considered typical for a session where we want to produce some code for a C++ implementation for CCD.
In [46]:
H = normal_ordered_hamiltonian()
T_2 = Operator([],[1,1,-1,-1])
expT = expand_ansatz([[T_2]], 4) #Natural truncation
tx = combine_to_excitation(H,expT, 0, [1,1,0,0]) #The list is explained below
S = "0 = "
for i in tx:
S+= "+" + i
Math(S)
Out[46]:
As expected, we find only one contribution to the energy expression. Now for the $t_2$ amplitude equation (in the example we return the code):
In [47]:
#tx = combine_to_excitation(H,expT, 2, [0,0,0,1]) #The list is explained below
tx = combine_to_excitation(H,expT, 2, [1,0,0,1]) #The list is explained below
S = "0 = "
for i in tx:
S+= "+" + i
Math(S)
Out[47]:
In summary, we see that CCAlgebra has a potential to simplify both the understanding and implementation of CC-equations. The code above may be copy-pasted into a C++ class file.
Update (2.november 2014): The code above was successfully implemented into my quantum many-body solver Fermion Mingle in C++ with some minor editing, and the code produced the exact same results I have obtained from previous benchmarking. This is a proof-of-concept, showing that CCAlg may be used to write code for solver functions.
The CCD is considered to be the simplest of the CC equations. When we also include single excitations, we may derive the CCSD equations. As the correlation energy is given by completely closed diagrams, we should only find three possible such connections of operators. This is confirmed in the example below:
In [19]:
H = normal_ordered_hamiltonian() #Including one- and two-particle interactions
expT = expand_ansatz([[T_1],[T_2]],3) #Taylor expand a list of lists to the 3rd order
tx = combine_to_excitation(H,expT,0, [1,1,0,0])
S = "0 = "
for i in tx:
S += "+" + i
Math(S)
Out[19]:
The expression above is identical to the energy expression presented in S-B, p.297, meaning that CCAlgebra at least solves the energy equation correctly. We proceed by deriving the $t_2$ amplitude equation:
In [20]:
tx = combine_to_excitation(H,expT,2, [1,0,0,0])
S = "0 = "
for i in tx:
S += "+" + i
Math(S)
Out[20]:
And the $t_1$ amplitude equation:
In [21]:
tx = combine_to_excitation(H,expT,1, [1,0,0,0])
S = "0 = "
for i in tx:
S += "+" + i
Math(S)
Out[21]:
Finally we derive the CCSDT equations. This exercise shows some of the potential of using CCAlgebra, as the equations are notoriously prone to error when derived by hand due to the length and complexity of the equations.
We begin by a confirmation that the energy equation remains the same as before (as no triple amplitudes will enter the energy expression):
In [22]:
H = normal_ordered_hamiltonian() #Including one- and two-particle interactions
T_1 = Operator([],[1,-1]) #The T_1 cluster operator
T_2 = Operator([],[1,1,-1,-1]) #The T_2 operator; all lists must be normal ordered
T_3 = Operator([],[1,1,1,-1,-1,-1])
expT = expand_ansatz([[T_1],[T_2], [T_3]],4) #Taylor expand a list of lists to the 3rd order
tx = combine_to_excitation(H,expT,0, [1,0,0,0])
s = "0 ="
for t in tx:
s += "+" + t + "\n"
Math(s)
Out[22]:
We then derive each amplitude equation in turn. For the $t_1$ amplitude:
In [23]:
tx = combine_to_excitation(H,expT,1, [1,0,0,0])
s = "0 ="
for t in tx:
s += "+" + t + "\n"
Math(s)
Out[23]:
For the $t_2$ amplitudes:
In [24]:
tx = combine_to_excitation(H,expT,2, [1,0,0,0])
s = "0 ="
for t in tx:
s += "+" + t + "\n"
Math(s)
Out[24]:
An finally, the $t_3$ amplitude:
In [25]:
tx = combine_to_excitation(H,expT,3, [1,0,0,0])
s = "0 ="
for t in tx:
s += "+" + t + "\n"
Math(s)
Out[25]:
In [26]:
H = Operator([1,-1,-1],[-1])
contraction = O(H, [T_1, T_2])
for i in range(contraction.N):
print contraction.code(i)
In [26]:
In [27]:
contraction.diagram(0, [0,0], True)
contraction.diagram(1, [0,0], True)
contraction.diagram(2, [0,0], True)
The full CCD equations is easily derived by higher order functions:
In [28]:
H = normal_ordered_hamiltonian() #Including one- and two-particle interactions
expT = expand_ansatz([[T_2]],3) #Taylor expand a list of lists to the 3rd order
Find energy expressions:
In [29]:
tx = combine_to_excitation(H,expT,0, [1,1,0,0]) #find all combinations of elements in H and expT that produce final excitation = 0
In [30]:
s = "0 ="
for t in tx:
s += "+" + t + "\n"
Math(s)
Out[30]:
In [31]:
combine_to_excitation(H,expT,1, [0,1,0,0]) #These will not contribute in CCD due to excitation level = 1
Out[31]:
In [32]:
combine_to_excitation(H,expT,2, [0,1,0,0]) #These will however contribute to the t2 amplitude equation due to excitation level = 2
Out[32]:
In [33]:
#The t2 amplitude equation is
tx = combine_to_excitation(H,expT,2, [1,0,0,0])
s = "0 ="
for t in tx:
s += "+" + t + "\n"
Math(s)
Out[33]:
In [34]:
#This means we have to write the following code for the t2 amplitude eq. of CCsolve in C++:
combine_to_excitation(H,expT,2, [0,0,0,1])
Out[34]:
In [35]:
#While the following code will give the correlation energy:
combine_to_excitation(H,expT,0, [0,0,0,1])
Out[35]:
The CCSD equations may be derived in the same manner:
In [36]:
H = normal_ordered_hamiltonian() #Including one- and two-particle interactions
expT = expand_ansatz([[T_1],[T_2]],2) #Taylor expand a list of lists to the 3rd order
#Note: running the command below will likely cause a memory overload due to the length of some of the distinct diagrams representation as strings of C-code.
#tx = combine_all(H,expT, [2,0,0,0])
s = "0 ="
for t in tx:
s += "+" + t + "\n"
Math(s)
Out[36]:
All diagrams may be plotted in a common display frame:
In [37]:
combine_to_excitation(H,expT,0, [0,1,0,0]) #These will however contribute to the t2 amplitude equation due to excitation level = 2
Out[37]:
In [38]:
combine_to_excitation(H,expT,1, [0,2,0,0])
Out[38]:
In [39]:
combine_to_excitation(H,expT,2, [0,2,0,0])
Out[39]:
In [40]:
H = normal_ordered_hamiltonian() #Including one- and two-particle interactions
T_1 = Operator([],[1,-1]) #The T_1 cluster operator
T_2 = Operator([],[1,1,-1,-1]) #The T_2 operator; all lists must be normal ordered
T_3 = Operator([],[1,1,1,-1,-1,-1])
expT = expand_ansatz([[T_1],[T_2], [T_3]],4) #Taylor expand a list of lists to the 3rd order
#print expT
tx = combine_to_excitation(H,expT,3, [1,0,0,0])
s = "0 ="
for t in tx:
s += "+" + t + "\n"
Math(s)
Out[40]:
In [41]:
print len(expand_ansatz([["a"],["b"],["c"]], 4))
In [41]: