In [39]:

from IPython.display import display, Math, Latex, Image, HTML
%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
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)" % ("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)" % ("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


Reveal.js 3 Slide Demo

# The Coupled Cluster method:

## Motivation, derivation and applications

Exam presentation | Audun Skau Hansen | UiO |  December 18th 2014

# Motivating the Coupled Cluster method

## The problem

• We have a system of interacting fermions, possibly in a common potential.
• We seek a wavefunction that describes and gives us the properties of this system:
• Energy:
• Ground state
• Ionization
• Excited states
• Fermionic density distribution
• Observables

## The strategy

• We formulate a strategy for solving this problem:
• We seek a general wavefunction that is a solution to the time-independent nonrelativistic Schrödinger equation.
• The wavefunction must then
• behave according to experiment.
• behave according to the postulates of QM
• We will employ the variational principle to optimize this general wavefunction.
• We will then seek methods that improve upon this first approximation.
• In the process, a number of assumptions about the problem will be made.

## Many particle wavefunctions

• We have a wavefunction assosciated with the system: $$|\Psi \rangle$$
• We also have wavefunctions assosciated with each fermion $$|\psi_i(x_i) \rangle$$ for the $i$th fermion, centered at $x_i$.

It is reasonable to assume that the systems wavefunction will be a function of each fermionic function: $$|\Psi \rangle = |\Psi(\psi_0(x_0), \psi_1(x_1), ..., \psi_N(x_N))\rangle$$

## Composite systems I

• Each constituent fermionic wavefunction exists in Hilbert space.
• Each constituent fermionic wavefunction is a solution to the corresponding one-body problem (where all other fermions are removed).
• The composite system will have a state space existing in the space spanned by the tensor product of each constituent space
$$\mathcal{H}_1 \otimes \mathcal{H}_2 \otimes \mathcal{H}_3 \otimes ...$$
• A subspace of the above tensor product might also suffice for the composite system.

## Composite systems II

When setting up a composite system for a number of subsystems, we then naturally assume

$$\Psi_h = \psi_0 \otimes \psi_1 \otimes ... = \psi_0(x_0) \psi_1(x_1) ...$$

This is called the Hartree function or -product. (All individual quantum numbers are contained in $x_n$)

The Hartree product is completely uncorrelated, meaning that the probability of finding fermions simultaneously at locations $x_0,x_1,...$ is simply

$$|\psi_i(x_0)\psi_j(x_1)...|^2 dx_0 dx_1 ... = |\psi_i(x_0)|^2 dx_0 |\psi_j(x_1)|^2 dx_2 ...$$

i.e. the product of each constituent particle wavefunction squared. The motion of these particles is independent of each other.

## Indistinguishability

• Fermions are experimentally known to be identical particles with half integer spin.
• Being identical is equivalent to being experimentally indistinguishable. [1]
• To take this into account, we must assume that the composite wavefunction is the same for all possible permutations of particles
$$\Psi_p = \frac{1}{N} \sum_{i<j} P_{ij} \Psi_h$$

Where P is the permutation operator and N is some normalizing factor.

(If we consider identical particle wavefunctions with identical quantum numbers, this property is already inherent in the Hartree function, since any permutation will just yield an equivalent expression of the Hartree function.)

Consequence:

• The systems wavefunction may be considered a linear combination of spatially permuted states; measurement will be unaffected by permutations.

[1] An exception is when fermions are being prepared prior to experiment, for example by being given opposite momentum.

## The Pauli Exclusion Principle

We know from experiment that identical fermions obey the Pauli principle:

• Assumtion: No two identical fermions may occupy the same QM state.

At least no exception for this principle has been found. [Citation]

## Taking the Pauli Principle into account I

The composite function should fulfill

$$\langle \Psi | \Psi \rangle = 0$$

when constituent particles have identical quantum numbers. The permuted Hartree function does not fulfill this condition, as can be seen by considering two particles with the same quantum numbers:

$$\Psi(a,b) = \frac{1}{\sqrt{2}} (\psi_a(a) \psi_b(b) + \psi_a(b)\psi_b(a))$$$$\rightarrow \langle \Psi(a,a) | \Psi(a,a) \rangle \neq 0$$

## Taking the Pauli Principle into account II

If we instead consider the antisymmetric linear combination of permuted states:

$$\Psi(a, b) = \frac{1}{\sqrt{2}}(\psi_a(a)\psi_b(b) - \psi_a(b)\psi_b(a))$$

we will find that the wavefunction behaves as we want it to:

$$\langle \Psi(a,a) | \Psi(a,a) \rangle = 0$$$$|\Psi(a,b)|^2 dx_a dx_b = |-\Psi(b,a)|^2 dx_a dx_b$$

## The Slater Determinant (SD)

The results so far is generalized in the Slater determinant:

$$\Psi_{SD} (\mathbf{x_1},\mathbf{x_2},...,\mathbf{x_N}) = \ \frac{1}{\sqrt{N!}} \begin{vmatrix} \psi_1(\mathbf{x_1}) & \psi_2(\mathbf{x_1}) & \cdots & \psi_N(\mathbf{x_1}) \\ \psi_1(\mathbf{x_2}) & \psi_2(\mathbf{x_2}) & \cdots & \psi_N(\mathbf{x_2}) \\ \vdots & \vdots & \ddots & \vdots \\ \psi_1(\mathbf{x_N}) & \psi_2(\mathbf{x_N}) & \cdots & \psi_N(\mathbf{x_N}) \\ \end{vmatrix}$$

Another equivalent notation for the SD is to use the antisymmetrizer

$$\Psi_{SD} = \frac{1}{\sqrt{N!}}(-)^p \hat{P} \prod_{i=1}^N \psi_i(x_i) \equiv \sqrt{N!} \mathcal{A} \Psi_h$$

This notation makes it easier to manipulate the expressions involved.

Some observations:

• Given that the constituent wavefunctions fulfills the Scrödinger equation for the system, so will the SD built by these functions.
• For an infinite single-particle basis, we may construct an infinite number of SD's.
• In the limit of an infinite number of SD's, any N-fermion function can be represented exactly.

## The Fermi Vacuum

Each single particle wavefunction will typically have a groundstate, and a possibly infinite number of excited states.

The system as a whole will also have energy levels ("orbitals") that may be occupied by particles.

In the case where the N particle SD is built from the N lowest wavefunctions in the many body problem, we may define the Fermi level to be the energy level above the highest energy state in the SD.

In this context, we often refer to the SD as the reference state, and denote it by $|\Psi_0 \rangle$

An excited SD may then be constructed by replacing one of the states below the Fermi level with one above. In the particle-hole formalism, this means to create particle-hole state.

A notation that is often used for such excitations is

$$a^{\dagger}_a a_i |\Psi_0\rangle = |\Psi_i^a \rangle$$

Where a fermion is removed from state $i$ (below F) and placed in state $a$ (above F).

For states above the Fermi level; abcd...

For states below the Fermi level; ijkl...

For all states; pqrs...

## The exact wavefunction

We may assume that there exists some operator that transforms one single SD built of an infinite single particle basis into the exact wavefunction:

$$\Phi = \hat{\Omega} \Psi_{SD}$$

One possibility for an infinite basis set, is to find and expansion of excited SD's that exactly represents the system

$$\Phi = |\Psi_0\rangle + \sum_{ai} c_{i}^a | \Psi_i^a \rangle + \sum_{abij} c_{ij}^{ab} | \Psi_{ij}^{ab} \rangle + ...$$

In this case

$$\Omega = (1+ \hat{C} + \hat{C}_2 + ...)$$

This approach is called Full CI, and is unattainable whitout any truncation in the basis set or number of excitations.

## The Schrödinger Equation

The nonrelativistic time-independent Schrödinger equation is:

$$H|\Psi_{SD} \rangle = E |\Psi_{SD} \rangle$$
• We may have a number of contributions to the hamiltonian:
• Fermion-fermion interaction
• Fermion-potential interaction
• Internal contributions from the common potential
• Kinetic energy

## The Variational Principle

The standard procedure when evaluating the energy of a solution to the Schrödinger equation is to multiply by $\langle \Psi_{0}|$ from the left:

$$E_0 \leq \langle \Psi_{0} |H|\Psi_{0} \rangle$$

Note that because the single SD is only an approximation to the solution, we only find an upper limit to the ground state energy.

The variational principle states that any function will be an upper limit to the ground state energy.

## Considering Spin

We have previously mentioned that fermions in general have half-integer spin. This will affect our system:

• It doubles the number of states accessible to the fermions.
• Opposite-spin states are orthogonal to each other.

## Interactions

The particles in the system might be subject to a number of interactions that must be included in the energy expression.

These are generally referred to as N-body operators, where N relate to the number of particles involved in the interaction.

We may formulate these as operators acting on only certain particles in the product states:

One-body operator

$$\hat{F} = \sum_p \hat{f}(x_p)$$

Meaning that it only affects particle with quantum numbers $x_p$.

Two-body operator

$$\hat{V} = \sum_{pq} \hat{v}(x_p,x_q)$$

And so on.

The interaction with a external potential field[3] is an example of a one-body operator, while an interaction between the fermions is an example of a two-body operator.

## Interactions

The particles in the system might be subject to a number of interactions that must be included in the energy expression.

These are generally referred to as N-body operators, where N relate to the number of particles involved in the interaction.

We may formulate these as operators acting on only certain particles in the product:

One-body operator

On second quantized form, a general one body operator is written

$$\sum_{pq} f^p_q \hat{a}_p^\dagger \hat{a}_q$$

Two-body operator

On second quantized form, a general two body operator is

$$\hat{V} = \sum_{pqrs} v^{pq}_{rs} \hat{a}_p^\dagger \hat{a}_q^\dagger \hat{a}_r \hat{a}_s$$

And so on.

The interaction with a external potential field[3] is an example of a one-body operator, while an interaction between the fermions is an example of a two-body operator.

## Calculating the ground state energy I

Using properties of the antisymmetrizer, we find that

$$E_0 = \langle \sqrt{N!} \mathcal{A} \Psi_h |\hat{H} | \sqrt{N!} \mathcal{A} \Psi_h \rangle = N! \langle \Psi_h | \hat{H}|\mathcal{A}^2 \Psi_h \rangle = N! \langle \Psi_h | \hat{H}|\mathcal{A} \Psi_h \rangle = \langle \Psi_h | \hat{H}| \Psi_{0} \rangle$$

When evaluating the one-body part of the energy of the slater determinant, we therefore must evaluate expressions of the form

$$\langle \Psi_h|\hat{F} | \Psi_{0}\rangle = \sum_{p} \langle \Psi_h| f^p_p f(x_p) | \Psi_{0}\rangle = \sum_i f^i_i = \sum_i^N \epsilon_i$$

This is just the sum over each single particle energy, and follows from the orthogonality of the products in the reference state.

## Calculating the ground state energy II

In the case of two-body interactions, we must also consider cases of one permutation in the reference state:

$$\langle \Psi_h|\hat{V} | \Psi_{0} \rangle = \sum_{pq} \langle \Psi_h | v^{pq}_{rs} \hat{v}(x_p,x_q) | \Psi_0 \rangle = \langle i j | ( 1 - P_{ij}) | i j \rangle = v^{ij}_{ij} - v^{ij}_{ji} \equiv \langle ij | ij \rangle_{AS}$$

This is caused by the two-body operator's ability to act on terms in the SD that differ from the hartree product by only one permutation, making it identical to the hartree product.

The two terms arising from the two particle interaction have been given the names direct and exchange term, respectively.

When we work with electronic systems, the first term is often refered to as the coulomb term.

## Fermionic correlation

• The fact that the motion of each fermion should be dependent upon the other fermions in the system is commonly referred to as correlation.

• Some of this correlation is attributed to the Pauli exlusion principle, and is inherent in the SD.

• Other correlations will be caused by forces acting amongst the fermions, such as coulomb repulsion between electrons or nuclear forces between nucleons.

• The correlation energy is a quantity commonly used in quantum chemistry to describe the difference in energy between the Hartree Fock energy and the exact solution of the nonrelativistic Schrödinger equation. [S-B, p.8]

## The Hartree Fock method

In many-fermion systems, it is the dependence

$$\frac{1}{|r_i - r_j|}$$

in $\hat{v}(x_i, x_j)$ that poses the difficulty in solving the equations.

In the Hartree Fock method for electron systems, this problem is solved by two equivalent assumptions:

• The system can be represented by a single Slater determinant.
• Each electron interacts with a mean field set up by the other electrons.

Both of these assumtions is wrong [citation needed], but it turns out that this will give us an initial approximative solution that often accounts for more than 99% of the energy and 95% of the wavefunction. [S-B, p6].

This approximation may also serve as a starting point for post Hartree Fock calculations.

## The Fock matrix

To obtain the HF reference, the Fock matrix may be derived by first expanding each constituent function in the SD in a single particle basis.

We then find the expression for the energy

$$E[{C}] = \sum_{\alpha \beta}^A \sum_i^N C_{\alpha i}^* C_{\beta i} \langle \phi_{\alpha i} | \hat{f}(x_i) | \phi_{\beta i} + \frac{1}{2} \sum_{\alpha \beta \gamma \delta}^A \sum_{ij, i \neq j} \langle \phi_{\alpha i} \phi_{\beta j} |\hat{v}(x_i, x_j) | \phi_{\gamma, i} \phi_{\delta j} \rangle_{AS}$$

We want to minimize the energy as a function of the coefficients, while constraining the orthogonality of the basis functions. This means we have to minimize the functional

$$F[C, \lambda] = E[C] - \sum_i \lambda_i \sum_{\alpha}^A C_{i \alpha}^* C_{i \alpha}$$

So that

$$\frac{\partial}{\partial C_{i \alpha}} F[C, \lambda] = 0$$

By performing the differentiation, we find the Fock matrix

$$h^{HF}_{\alpha \beta} = \langle \alpha | f | \beta \rangle + \sum_{\gamma \delta}^A \sum_j^N C_{\beta j}^* C_{\gamma i} C_{\delta j} \langle \alpha \beta |\hat{v} | \gamma \delta \rangle_{AS}$$

## Hartree Fock - the Mean Field approximation

Hartree Fock has also been called the Self-Consistent Field method, and the Mean Field approximation.

When minimizing the energy, we basically replace the interaction between each fermion with a mean field, so that each fermion interacts with this field. We may interpret the Fock matrix as

$$h^{HF} = \hat{t}_{kin} + \hat{u}_{ext} + \hat{u}_{hf}$$

Where $\hat{u}_{hf}$ is the mean field.

[Slides, FYS-KJM4480]

Consequence

• Even for a complete basis, Hartree Fock will only provide an approximation.

By neglecting instantaneous interactions between the fermions, the Hartree Fock method will fail to account for parts of the correlation energy.

## The HF solution process

The typical solution procedure is to

1. Guess som initial coefficients.
2. Set up the Fock matrix.
3. Diagonalize the Fock matrix to obtain new coefficients.
4. Check for convergence, if not: calculate new Fock matrix, repeat 3.

When it converges, we will have obtained a set of eigenvectors of the Fock matrix, representing the single particle states entering the SD.

## Correlations in the HF method

Only part of the correlations occuring in the exact system will be accounted for in the HF approximation:

• The Fermi Correlation (from the Pauli principle).
• To some extent the coulomb interaction. (Mean field)

## Restricted Hartree Fock

For closed shell atoms and molecules, such as He, $H_2$ and Be, the fermions (in this case electrons), may be assumed to be paired with opposite spin particles.

In the restricted Hartree Fock method, we consider only systems where all fermions are paired with opposite spin particles.

## Unrestricted Hartree Fock

For systems with singly oocupied states (unpaired fermions, typically odd numbered) the restriction of spin-pairing is no longer valid.

To take this into account, we set up two Fock matrices; one for each spin orientation. The resulting system is

$$F^{\alpha}(C^{\alpha}, C^{\beta}) C^{\alpha} = S C^{\alpha} \epsilon_{\alpha}$$$$F^{\beta}(C^{\alpha}, C^{\beta}) C^{\beta} = S C^{\beta} \epsilon_{\beta}$$

The ground state energy will be a function of the eigenvectors of these two matrices.

(The dependencie of opposite spin coefficients in the Fock matrices is due to the direct- or coulomb term from the two particle integrals)

[Modern Quantum Chemistry, Szabo, p214]

## Spin contamination

The $\hat{S}^2$ operator commutes with the hamiltonian $\hat{H}$, meaning that they have a common set of eigenfunctions.

When performing UHF this might not always be the case, since we treat fermions that might otherwise occupy the same orbital differently.

[Modern Quantum Chemistry, Szabo]

• Thouless theorem
• Brillouins theorem

• The disassosciation problem

• Avoids spin contamination
• Is size extensive
• Takes into account instantaneous interactions

## Motivation for post HF methods, Summary

• A single determinant is not necessarily a good approximation to the true wavefunction.

• The lack of instantaneous coulumb repulsion in HF

• Spin contamination might occur in UHF and needs correction.

# Deriving The Coupled Cluster method

## The exponential ansatz

The exponential ansatz is an operator acting on the reference state:

$$e^{\hat{T}} |\Psi_o \rangle$$

Where the cluster operators

$$\hat{T} = 1 + \hat{T}_1 + \hat{T_2} + ... = 1 + \sum_{ai} t_i^a a_{a}^\dagger a_i + \sum_{ai} t_{ij}^{ab} a_{a}^\dagger a_{b}^\dagger a_i a_j + ...$$

The coefficients $t_i^a$ are called amplitudes, and will need to be solved for. As for any exponential function, we may expand it as

$$e^{\hat{T}} = 1 + \hat{T} + \frac{1}{2!} \hat{T}^2 + \frac{1}{3!} \hat{T}^3 + ...$$

So that

$$e^{\hat{T}} |\Psi_0 \rangle = (1 + \hat{T} + \frac{1}{2!} \hat{T}^2 + \frac{1}{3!} \hat{T}^3 + ...)|\Psi_0 \rangle$$

## The cluster operator

As can be seen from the second quantized form of the cluster operators, when acting on the reference they will cause an excitation to n-particle n-hole state according to their order.

$$\hat{T}_1|\Psi_0\rangle = \sum_{ai} t_i^a a_{a}^\dagger a_i|\Psi_0\rangle = \sum_{ai} t_i^a|\Psi_i^a\rangle$$$$\hat{T}_2|\Psi_0\rangle = \sum_{abij} t_{ij}^{ab} a_{a}^\dagger a_{b}^\dagger a_i a_j|\Psi_0\rangle = \sum_{abij} t_{ij}^{ab}|\Psi_{ij}^{ab}\rangle$$

## Extensivity and consistency

A quantum mechanical model may be called extensive if the energy of this model scales correctly with the size of the system. [S-B, p11]

Size consistency means that the energy of two subsystems combined or apart should, for large separations, fulfill

$$E(AB) = E(A) + E(B)$$

[S-B, p12]

## Extensivity

The extensivity property is built into the exponential ansatz.

IN S-B, p255, it is shown that if the zero-order wavefunction of a noninteracting system is separable

$$\phi_0(A,B) = \phi_0(A)\psi_0(B)$$

And $\hat{T}$ for the system is additive

$$\hat{T}(A,B) = \hat{T}(A) + \hat{T}(B)$$

It follows that the total wavefunction is multiplicatively separable:

$$\Psi(A,B) = e^{T(A) + T(B)}\phi_0(A,B) = e^{T(A)}\phi_0(A)e^{T(B)}\phi_0(B)$$

And so the energy is additive

$$\hat{H}(A,B)\Psi(A,B) = [E(A) + E(B)]\Psi(A,B)$$

## Size consistency

When comparing CC to CI, we see that the cluster operator corresponding to the 2p2h SD's in CI is

$$\hat{C}_2 = \hat{T}_2 + \frac{1}{2}(\hat{T}_1)^2$$

The last term is commonly referred to as a disconnected cluster.

In a system for two particles A and B, but separated so that practically no interaction occurs, the coupled cluster energies for this system will yield

$$E_{CC} = E_{A,CC} + E_{B,CC}$$

This is not the case for CI

$$E_{CI} \neq E_{A,CI} + E_{B,CI}$$

[Crawford and Schaefer, p18]

• Comparison with CI
• Disconnected vs. connected cluster

## The CC equations I

We now want to derive the equations that lets us find the amplitudes and the energy resulting from the exponential ansatz.

We begin by observing that the normal product form of the Schrödinger equation for the exponential ansatz can be written

$$(\hat{H}_N - \Delta e) e^{\hat{T}} |\Psi_0 \rangle = 0$$

We multiply by the inverse exponential to find

$$(e^{\hat{-T}} \hat{H}_N e^{\hat{T}} - \Delta e) |\Psi_0 \rangle = 0$$

The non-Hermitian operator

$$e^{\hat{-T}} \hat{H}_N e^{\hat{T}} \equiv \mathcal{H}$$

is a similarity transformed Hamiltonian, often called the CC effective Hamiltonian.

## The CC equations II

If we apply the Baker-Campbell-Hausdorff expansion [S-B, p293] to the CC effective hamiltionian, we may rewrite it as

$$\mathcal{H} = \hat{H}_N + [\hat{H}_N, \hat{T}] + \frac{1}{2} [[\hat{H}_N, \hat{T}], \hat{T}] + \frac{1}{3!} [[[\hat{H}_N, \hat{T}], \hat{T}], \hat{T}] + \frac{1}{4!} [[[[\hat{H}_N, \hat{T}], \hat{T}], \hat{T}], \hat{T}].$$

The reason that this sum of nested operators truncates at the four-fold commutator, becomes apparent if we use Wicks generalized theorem to evaluate the expansion.

Following Shavitt and Bartlett (p.294) we have for the commutator of two operators with even number of creation and annihilation operators;

$$[A,B] = AB - BA = np(A,B) + (\dot{A}\dot{B}) - np(B,A) - (\dot{B}\dot{A})$$

Where $np$ indicates a normal ordered product of operators and dotted operators indicate a sum of all possible contractions between the normal ordered operators A and B.

This is simplified further since both A and B (H and T) contain an even number of creation and annihilation operators:

$$np(A,B) = np(B,a)$$

## The CC equations III

We end up with

$$[A,B] = (\dot{A}\dot{B}) - (\dot{B}\dot{A})$$

Since the cluster operators $T_m$ commute, the nested operators will only have surviving terms where $\hat{H}_N$ contracts with one or more cluster operators.

The reason for the truncation of the Hausdorff expansion is then apparent; $\hat{H}_N$ is a sum of operators which at most contains four creation or annihilation operators. It is therefore not possible to contract any term in $\hat{H}_N$ to more than 4 cluster operators.

## The CC equations IV

We also have to take into account that particle creation operators can only produce a nonzero contraction with a particle annihilation operator to its left, and a hole annihilation operator may only produce a nonzero contraction with a hole creation operator to its left, so that the only contributing terms when unraveling the nested operator is terms that begin with the normal ordered hamiltonian:

$$\mathcal{H} = e^{\hat{-T}} \hat{H}_N e^{\hat{T}} = \hat{H}_N + (\dot{\hat{H}_N}\dot{\hat{T}}) + \frac{1}{2}(\dot{\hat{H}_N}\dot{\hat{T}}\dot{\hat{T}}) + \frac{1}{3!}(\dot{\hat{H}_N}\dot{\hat{T}}\dot{\hat{T}}\dot{\hat{T}}) + \frac{1}{4!}(\dot{\hat{H}_N}\dot{\hat{T}}\dot{\hat{T}}\dot{\hat{T}}\dot{\hat{T}})\equiv (\hat{H}_Ne^{\hat{T}})_C$$

The dots indicate a sum over all terms in which the hamiltonian connect by at least one contraction to each of the cluster operators to its right. The subscript C indicates connected terms.

## The CC equations V

We are finally in a position to derive explicit expressions for the energy and amplitudes.

We have the general equation

$$(\hat{H}_Ne^{\hat{T}})_C |\Phi_0 \rangle = \Delta e |\Phi_0 \rangle$$

Finding the energy is now simply

$$\langle \Phi_0 |(\hat{H}_Ne^{\hat{T}})_C |\Phi_0 \rangle = \Delta e$$

The amplitudes may be found from the equations

$$\langle \Phi_i^a |(\hat{H}_Ne^{\hat{T}})_C |\Phi_0 \rangle = 0$$$$\langle \Phi_{ij}^{ab} |(\hat{H}_Ne^{\hat{T}})_C |\Phi_0 \rangle = 0$$

...

What remains now is the tedious process of working out the explicit equations for a give truncation in the exponential ansatz.

## Truncating the ansatz

As previously stated, we have

$$\hat{T} = \hat{T}_1 + \hat{T}_2 + \hat{T}_3 + ...$$

We truncate this ansatz by limiting the number of cluster operators entering the exponent. The standard truncations is

CCD: $\hat{T} = \hat{T}_2$ (doubles)

CCSD: $\hat{T} = \hat{T}_1 + \hat{T}_2$ (singles, doubles)

CCSDT: $\hat{T} = \hat{T}_1+ \hat{T}_2 + \hat{T}_3$ (singles, doubles, triples)

Thouless Theorem states that a CCS approach would only transform a single determinant into another single determinant, so this truncation does not occur. (In analogy to Brillouin's theorem in CI).

## The explicit form of the equation

We could now proceed by working out the explicit equations for amplitudes and energy using Wick's theorem.

Another approach is however outlined in S-B (p.297) that makes the process much simpler and interesting, and it is based on the application of diagrams.

## The diagrammatic normal ordered hamiltonian

(Illustration from S-B, p. 297)

The normal ordered hamiltonian is a sum of "filaments" that may be contracted with cluster operators to form diagrams.

The + and - signs below the interaction line indicates virtual particle annihilation operators. They will be helpful when we evaluate the different ways of contracting each diagrams.

## The diagrammatic cluster operators

(Illustration from S-B, p. 298)

The diagrammatic cluster operators of different excitation levels.

Vertical lines between the + and - symbols may be used to denote repeated operators.

## Connecting operators

(Illustration from S-B, p.299)

In the example above, the contraction

$$\langle \Phi_0 | \hat{H} \hat{T_2} \hat{T_2} | \Phi_0 \rangle$$

is evaluated.

## Interpreting diagrams

In Shavitt Bartlett, 10 rules for intepreting coupled-cluster diagrams are given:

(1)  Label internal and external lines with particle and hole labels
(2)  Assosciate f(i,a) with every 1-p operator vertex
(3)  Assosciate <lout rout || lin rin> with every 2-p operator ver
(4)  Assosciate t(ab,ij) for every amplitude
(5)  Sum over all internal lines
(6)  Multiply by 1/2 for each pair of equivalent internal lines
(7)  Multiply by 1/2 for each pair of equivalent T-vertices
(8)  Include a phase factor (-)**(holes-loops)
(9)  Sum over all distinct permutations of inequivalent external pairs
(10) Cancel each factor .5 in open diagrams with equivalent vertives arising from rule 7 with a permutation of the labels of a pair of external lines connected to the equivalent vertices.

## Implementing diagrams

Although the diagrammatic rules greatly simplifies the process of deriving the equations, there is still a high possibility for errors due to human inaccuracy.

The CCSD T1 amplitude have for example 14 terms, making even the diagrammatic process quite tedious.

## In the words of H.G.Kümmel

"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

[A biography of the coupled cluster method, June 2001]

## One more simplification

To simplify the process of deriving explicit CC equations, I have developed a python script called CCAlgebra (Coupled Cluster Algebra).

The intention is to use the diagrammatical approach in an object oriented context to automate the derivation of the equations involved. It is basically an attempt to teach the computer to derive and interpret diagrams, circumventing the use of wicks theorem or second quantization manipulations.

In the remaining slides I will give an overview of this script, and show how to use it in code developement.

## Defining operators



In [40]:

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



The following operators are defined above:

## Find all possible contractions



In [41]:

contraction = O(H,[T_1]) #The second argument must be a list of operators (T_1 * T_1 = [T_1, T_1])



## Show excitation and number of resulting diagrams



In [42]:

print "Excitation:", contraction.E() #Excitation energy of connected operators
print "Number of distinct diagrams:", contraction.N




Excitation: 0.0
Number of distinct diagrams: 1



## Show a mathematical representation (latex)



In [43]:

c_tex = contraction.latex(0) #The first of possibly multiple distinct diagrams in latex format.
Math(c_tex)




Out[43]:

$$\frac{1}{1} \sum_{ck} \langle k || c \rangle t_{k}^{c}$$



## Display a diagrammatic representation



In [44]:

contraction.diagram(0, [0,0], True) #Display diagram 0 at position x/y = [0,0] using internal show() function (matplotlib)






## A more complex example



In [45]:

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[45]:

$$P(ba)\frac{-1}{1} \sum_{ckd} \langle ak || cd \rangle t_{k}^{c} t_{ij}^{db}$$




In [46]:

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()






## Higher level functions

It is easy to develop higher level operations on top of this framework.

Considering that the exponent in the ansatz will typically be a sum of cluster operators, a neat function is to expand a list of operator objects:



In [47]:

print expand_ansatz([["a"],["b"]], 3)




[['a'], ['b'], ['a', 'a'], ['a', 'b'], ['b', 'a'], ['b', 'b'], ['a', 'a', 'a'], ['a', 'a', 'b'], ['a', 'b', 'a'], ['a', 'b', 'b'], ['b', 'a', 'a'], ['b', 'a', 'b'], ['b', 'b', 'a'], ['b', 'b', 'b']]



A predefined normal ordered hamiltonian will also come in handy



In [48]:

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]



## Cluster operators

A cluster operator is also predefined



In [49]:

print cluster_operator([1,2,3]) # returns [T_1, T_2, T_3]




[[<__main__.Operator instance at 0x10c015e60>], [<__main__.Operator instance at 0x10c015440>], [<__main__.Operator instance at 0x10bd9b8c0>]]



## Evaluating sums of operators

The normal ordered hamiltonian and the cluster expansion is both sums of operators. We will need to combine these in all possible ways, at least find all combinations that correspond to a given energy level.

Below we derive the CC energy



In [50]:

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]], 4)

tx = combine_to_excitation(H,expT,0, [1,0,0,0]) #All combinations that produce excitation level 1
S = "0 ="
for i in tx:
S+= "+"+ i
Math(S)




Out[50]:

$$0 =+\frac{1}{1} \sum_{ck} \langle k || c \rangle t_{k}^{c}+\frac{1}{4} \sum_{cdkl} \langle kl || cd \rangle t_{kl}^{cd}+\frac{1}{2} \sum_{ckdl} \langle kl || cd \rangle t_{k}^{c} t_{l}^{d}$$



## Evaluating the energy

It is easy to recognize the energy when viewed as diagrams



In [51]:

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]], 4)

tx = combine_to_excitation(H,expT,0, [0,1,0,0]) #All combinations that produce excitation level 0






## The CCD equations

To show some of the potential in such a code, we will now develop the full CCD equation.



In [52]:

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[52]:

$$0 = +\frac{1}{4} \sum_{cdkl} \langle kl || cd \rangle t_{kl}^{cd}$$




In [53]:

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, 2, [1,1,0,0]) #The list is explained below
S = "0 = "
for i in tx:
S+= "+" + i + "\n"
Math(S)




Out[53]:

$$0 = +P(ba)\frac{-1}{1} \sum_{c} \langle a || c \rangle t_{ij}^{cb} +P(ij)\frac{1}{1} \sum_{k} \langle k || j \rangle t_{ik}^{ab} +\frac{1}{2} \sum_{cd} \langle ab || cd \rangle t_{ij}^{cd} +\frac{1}{2} \sum_{kl} \langle kl || ij \rangle t_{kl}^{ab} +P(ba)P(ij)\frac{-1}{1} \sum_{ck} \langle ak || cj \rangle t_{ik}^{cb} +P(ij)\frac{-1}{2} \sum_{cdkl} \langle kl || cd \rangle t_{ik}^{cd} t_{jl}^{ab} +\frac{1}{4} \sum_{cdkl} \langle kl || cd \rangle t_{ij}^{cd} t_{kl}^{ab} +P(ab)\frac{-1}{2} \sum_{ckld} \langle kl || cd \rangle t_{kl}^{ca} t_{ij}^{db} +P(ab)P(ij)\frac{1}{2} \sum_{ckdl} \langle kl || cd \rangle t_{ik}^{ca} t_{jl}^{db}$$



## Generating code

When the diagrammatic representation is understood by the computer, it is easy to translate it to latex, diagrams, text or even code.

In the example below, we generate a C++ code directly for the T2 CCD amplitude instead of diagrams or latex:



In [54]:

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, 2, [1,0,0,1]) #The list is explained below




double CC0_0a = 0.0;
for(int c = nElectrons; c < nStates; c ++){
CC0_0a += v1(a)(c)*tf2(c,b)(i,j)-(v1(b)(c)*tf2(c,a)(i,j));
}
CC0_0a *= -1.000000;

double CC1_0a = 0.0;
for(int k = 0; k < nElectrons; k ++){
CC1_0a += v1(k)(j)*tf2(a,b)(i,k)-(v1(k)(i)*tf2(a,b)(j,k));
}
CC1_0a *= 1.000000;

double CC4_0a = 0.0;
for(int d = nElectrons; d < nStates; d ++){
for(int c = nElectrons; c < nStates; c ++){
CC4_0a += v2(a,b)(c,d)*tf2(c,d)(i,j);
}
}
CC4_0a *= 0.500000;

double CC5_0a = 0.0;
for(int l = 0; l < nElectrons; l ++){
for(int k = 0; k < nElectrons; k ++){
CC5_0a += v2(k,l)(i,j)*tf2(a,b)(k,l);
}
}
CC5_0a *= 0.500000;

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

double CC12_1a = 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 ++){
CC12_1a += v2(k,l)(c,d)*tf2(c,d)(i,k)*tf2(a,b)(j,l)-(v2(k,l)(c,d)*tf2(c,d)(j,k)*tf2(a,b)(i,l));
}
}
}
}
CC12_1a *= -0.500000;

double CC12_1b = 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 ++){
CC12_1b += v2(k,l)(c,d)*tf2(c,d)(i,j)*tf2(a,b)(k,l);
}
}
}
}
CC12_1b *= 0.250000;

double CC12_1c = 0.0;
for(int d = nElectrons; d < nStates; d ++){
for(int l = 0; l < nElectrons; l ++){
for(int k = 0; k < nElectrons; k ++){
for(int c = nElectrons; c < nStates; c ++){
CC12_1c += v2(k,l)(c,d)*tf2(c,a)(k,l)*tf2(d,b)(i,j)-(v2(k,l)(c,d)*tf2(c,b)(k,l)*tf2(d,a)(i,j));
}
}
}
}
CC12_1c *= -0.500000;

double CC12_1d = 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 ++){
CC12_1d += v2(k,l)(c,d)*tf2(c,a)(i,k)*tf2(d,b)(j,l)-(v2(k,l)(c,d)*tf2(c,b)(i,k)*tf2(d,a)(j,l))-(v2(k,l)(c,d)*tf2(c,a)(j,k)*tf2(d,b)(i,l)-(v2(k,l)(c,d)*tf2(c,b)(j,k)*tf2(d,a)(i,l)));
}
}
}
}
CC12_1d *= 0.500000;



## Three body interactions?

In principle, extending the code to include three body interactions is trivial.

This functionality is however not yet implemented.



In [60]:

H   = Operator([1,1,1,-1,-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
expT = expand_ansatz([[T_2]], 4) #Natural truncation
tx = combine_to_excitation([H],expT,0, [0,1,0,1]) #All combinations that produce excitation level 1
S = "0 ="
for i in tx:
S+= "+"+ i
Math(S)




Out[60]:

$$0 =$$



## How CCAlgebra works

• Operators are created as classes
• The O() class will create a combined operator form two ore more operators.
• Generating all distinct combinations is done by implementing itertools:
• Set up strings of symbols representing the cluster operators:
• | + + - -
• Create all permutations of this string
• Sort out all non-contributing terms
• Cases where no contractions happen with certain clusters | + - + + - -
• | | + + - - -
• Cases that do not connect to the hamiltonian
• Indistinct permutations
• And so on
• When only distinct diagrams remain, these are returned as a list of "schemas" (or blueprints).
• These schemas may be passed to different classes, that sets up the desired representation:
• Diagrams
• Latex
• Code
• ... Other possibilities?

# Fermion Mingle

www.github.com/audunsh/FermionMingle

During this semester I have also written an CC extension to a Hartree Fock solver I wrote as part of the course FYS4411.

The solver is named "Fermion Mingle", and currently supports

• Gaussian basis sets and integral calculations
• Restricted Hartree-Fock
• Unrestricted Hartree-Fock
• Coupled Cluster (CCD, CCSD)
• Electron density evaluation

Some functions are implemented in python:

• CCalgebra - An environment for coupled cluster calculations and code generation
• TurbomoleConverter - A script for generating C++ code from turbomole files

While the solver in its current state is designed to do calculations on atoms and molecules, the intention is to make it more general and to include fermionic gas caluclations. (electron and neutron).

The code generated for CCD above was successfully copy-pasted into a function in this solver and reproduced exactly the handwritten code, which in turn have been benchmarked against the literature.

## Fermion Mingle environment

A number of classes is written for flexible calculations on a variety of atomic and molecular systems:

    cout << "Performing UHF on H2 using STO-6G set" << endl;
vec3 corePos1 = {0,0,0};                       //setting up some position vectors for the cores
vec3 corePos2 = {0,0,6.0};
fmingle myparty;                               //creating the system
myparty.add_nucleus(corePos1, 1);              //creating a nucleus with charge 1
myparty.add_nucleus(corePos2, 1);               //creating a nucleus with charge 1
myparty.fminglebasisbank.add_STO_3G_h(corePos1);//creating an electron at core 1 using STO-3G basis set
myparty.fminglebasisbank.add_STO_6G_h(corePos2);//creating an electron at core 2 using STO-6G basis set
//myparty.rhf_solve(2);                         //perform a rhf procedure for 2 electrons

myparty.uhf_solve(1,1);                         //perform a uhf procedure for config 1 up, 1 down
//myparty.ccsd_solve(2);                        //perform a CCSD procedure on the obtained reference state
cout << "UHF energy:" << myparty.uhf_energy << endl;
cout << "Correlation energy:" << myparty.correlation_energy << endl;

# A model experiment: H2

## Dissosciation energy

Studying the behaviour of the $H_2$ ground state energy with increasing bond lengths, will reveal some features of the methods discussed.

When the atoms (and electrons) are close together this system may be accurately described by RHF, since the two electrons will pair up with opposite spin in a single orbital.

Increasing the distance between the nucleuses will however make this description less accurate. For large distances, the artifical constraint of spin pairing will no longer be a good representation of the system, and we should expect the energy to diverge above the expected energy.

In reality the bond is expected to be broken for large distances, so that the energy should just be the sum of to separate systems.

I have done some calculations on this system using Fermion Mingle.

[Szabo, Modern Quantum Chemistry, p165]

## The amplitude coefficients

(Observation, 18.12.14: In regions with no T2 contribution, T1 amplitudes peaks to an order of $10^{-5}$ (from $10^{-15}$.)

## General result comparison

System Total energy (a.u.) Correlation (a.u.) Parameters Comparison
H2 -1.13728 -0.0205616 Bond length: 1.4, STO-3G, RHF+CCSD HF-limit: -1.132 (2), Exact correlation -0.03969 (1)
H2 -1.14593 -0.0206048 Bond length: 1.4, STO-6G, RHF+CCSD HF-limit: -1.132 (2), Exact correlation -0.03969 (1)
O2 -147.696 -0.145003 Bond lenght: 2.287, STO-3G, RHF+CCSD HF-limit: -, Exact correlation -0.37 (1)
O -73.7092 -0.047385 STO-3G, RHF+CCSD HF-limit-74.729 (2), Correlation: -0.262 (2)
Be -14.4037 -0.0517703 STO-3G, RHF+CCSD General HF: -14.67 (3)
Be -14.5199 -0.00484722 Hydrogenlike (4 functions), RHF+CCSD General HF: -14.67 (3)
He -2.84228 -0.00868762 Hydrogenlike (4 functions), RHF+CCSD General HF: -2.904 (3)

(1) Results from Thijssen, p84, obtained from variational calculus

(2) results from 460performance.pdf (see text above)

(3) Source given in "Report from FYS4411 - Gøran Brekke Svaland and Audun Skau Hansen"

## Future Prospects

A number of future prospects for Fermion Mingle:

• Improving convergence, possibly using relaxation.
• Full integration of CCAlgebra into Fermion Mingle (or the other way around)
• Periodic fermion gas calculations (following the work of G.Baardsen)
• Make calculation of properties using CC possible.
• Support for parallellized calculations
• Optimizing atomic-to-molecular transformation (following the work of O.T.Norli)
• Support for UHF-CCSD
• Optimizing the CC method using intermediates

# End of presentation

## The water molecule

We have previously done some calculations on H2O using RHF and STO-3G to determine geometry and ground state energy. The results corresponded rather nice to the ones presented in the literature (see Report from FYS4411 - Gøran Brekke Svaland and Audun Skau Hansen).

We now use the geometry and energy presented in the following document to compare our Coupled Cluster Singles Doubles (CCSD) solver to the corresponding calculations performed by NWChem: http://institute.loni.org/NWChem2012/documents/tce-session.pdf.

Prior to the CCSD calculation the different algorithms starts out with the following parameters:

HF-energy (NWChem) : -74.962663062148

HF-energy (CCSolve): -74.962677575226

HF-limit : -76.067 (http://chemistry.illinoisstate.edu/standard/che460/handouts/460performance.pdf)

The results from each iteration is presented in the table below. The rightmost coulomn contains the results obtained by CCSolve.

Iter Correlation(NWChem) Correlation (CCsolve) Rel.Error
1 -0.0358672469179 -0.035897626185118 0.0008469919
2 -0.0454068882657 -0.04544804888198 0.000906484
3 -0.0483870059027 -0.048432634274069 0.0009429881
4 -0.0494370597647 -0.049484831256215 0.0009663093
5 -0.0498391184890 -0.049888079084351 0.0009823728
6 -0.0500021724029 -0.050051825093841 0.0009930107
7 -0.0500711904756 -0.050121251581979 0.0009997986
8 -0.0501014381364 -0.050151741527804 0.0010040309
9 -0.0501150974135 -0.05016554465157 0.0010066276
10 -0.0501214303300 -0.050171962746426 0.0010081998
11 -0.0501244348663 -0.050175017473014 0.0010091407
12 -0.0501258887096 -0.050176500709333 0.0010096978
13 -0.0501266039080 -0.050177233008851 0.0010100246
14 -0.0501269605251 -0.050177599511074 0.0010102146
15 -0.0501271402835 -0.050177784948453 0.0010103242
16 -0.0501272316751 -0.050177879584002 0.0010103871
17 -0.0501272784536 -0.05017792820593 0.0010104229
18 -0.0501273025229 -0.050177953317881 0.0010104433
19 -0.0501273149581 -0.050177966340254 0.0010104547
20 -0.0501273214031 (converged) (none)


In [ ]: