Commit e61bdfb4 authored by Davison, Jacob's avatar Davison, Jacob
Browse files

resolve merge conflict in main

parents 66141575 3bfe1f0d
......@@ -10,6 +10,8 @@ from scipy.integrate import odeint, ode
import numpy as np
import time
import pickle
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import tracemalloc
......@@ -20,7 +22,8 @@ import random
import tensornetwork as tn
#tn.set_default_backend("tensorflow")
# user files
# USER MODULES
from oop_imsrg.spin_sq import *
from oop_imsrg.hamiltonian import *
from oop_imsrg.occupation_tensors import *
from oop_imsrg.generator import *
......@@ -30,6 +33,11 @@ from oop_imsrg.plot_data import *
import oop_imsrg.ci_pairing.cipy_pairing_plus_ph as ci_matrix
from oop_imsrg.tests2B import *
sys.path.append('/mnt/home/daviso53/Research/')
from pyci.density_matrix.density_matrix import density_1b, density_2b
import pyci.imsrg_ci.pyci_p3h as pyci
#import reference_state_ensemble.reference_ensemble as re
def get_vacuum_coeffs(E, f, G, basis, holes):
H2B = G
......@@ -43,12 +51,24 @@ def get_vacuum_coeffs(E, f, G, basis, holes):
H0B = E - np.trace(H1B[np.ix_(holes,holes)]) - 0.5*result_ij.tensor
#print(result_ij.tensor, result_ji.tensor)
return (H0B, H1B, H2B)
def decouple_od(vo1b, vo2b, particles, holes):
vo1bod = np.zeros_like(vo1b)
vo1bod[np.ix_(particles, holes)] += vo1b[np.ix_(particles, holes)]
vo1bod[np.ix_(holes, particles)] += vo1b[np.ix_(holes, particles)]
vo1bd = vo1b - vo1bod
vo2bod = np.zeros_like(vo2b)
vo2bod[np.ix_(particles, particles, holes, holes)] += vo2b[np.ix_(particles, particles, holes, holes)]
vo2bod[np.ix_(holes, holes, particles, particles)] += vo2b[np.ix_(holes, holes, particles, particles)]
vo2bd = vo2b - vo2bod
return (vo1bd, vo1bod, vo2bd, vo2bod)
# @profile
def derivative(t, y, hamiltonian, occ_tensors, generator, flow):
def derivative(t, y, inputs):
"""Defines the derivative to pass into ode object.
Arguments:
......@@ -66,32 +86,28 @@ def derivative(t, y, hamiltonian, occ_tensors, generator, flow):
dy -- next step in flow"""
assert isinstance(hamiltonian, Hamiltonian), "Arg 2 must be Hamiltonian object"
hamiltonian, occ_tensors, generator, flow = inputs
#assert isinstance(hamiltonian, Hamiltonian), "Arg 2 must be Hamiltonian object"
assert isinstance(occ_tensors, OccupationTensors), "Arg 3 must be OccupationTensors object"
assert isinstance(generator, Generator), "Arg 4 must be Generator object"
assert isinstance(flow, Flow), "Arg 5 must be a Flow object"
E, f, G = ravel(y, hamiltonian.n_sp_states)
# half = int(len(y)/2)
# timefi = time.time()
# Hamiltonian flow
E, f, G = ravel(y, hamiltonian.n_sp_states)
generator.f = f
# timeff = time.time()
# timegi = time.time()
generator.G = G
# timegf = time.time()
# timefli = time.time()
dE, df, dG = flow.flow(generator)
# timeflf = time.time()
dE, df, dG = flow.flow(f, G, generator)
# print("""
# f time: {:2.4f} s
# G time: {:2.4f} s
# Flow time: {:2.4f} s
# """.format(timeff-timefi, timegf-timegi, timeflf-timefli))
# Spin-squared flow
#E_spin, f_spin, G_spin = ravel(y[half::], tspinsq.n_sp_states)
# generator_spin.f = f_spin
# generator_spin.G = G_spin
#dE_spin, df_spin, dG_spin = flow_spin.flow(f_spin, G_spin, generator)
dy = unravel(dE, df, dG)
#dy = np.concatenate([unravel(dE, df, dG), unravel(dE_spin, df_spin, dG_spin)], axis=0)
dy = unravel(dE,df,dG)
return dy
......@@ -139,7 +155,7 @@ def ravel(y, bas_len):
return(ravel_E, ravel_f, ravel_G)
# @profile
def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1, flow_data_log=0, generator='wegner', output_root='.'):
def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0, verbose=1, flow_data_log=0, generator='wegner', output_root='.'):
"""Main method uses scipy.integrate.ode to solve the IMSRG(2) flow
equations.
......@@ -179,11 +195,13 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1, flow_dat
if not os.path.exists(output_root):
os.mkdir(output_root)
if ref == []:
if ref is None:
ha = PairingHamiltonian2B(n_holes, n_particles, d=d, g=g, pb=pb)
ref = ha.reference # this is just for printing
# ss = TSpinSq(n_holes, n_particles)
else:
ha = PairingHamiltonian2B(n_holes, n_particles, ref=ref, d=d, g=g, pb=pb)
# ss = TSpinSq(n_holes, n_particles)
ot = OccupationTensors(ha.sp_basis, ha.reference)
......@@ -196,6 +214,9 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1, flow_dat
wg = generator_dict[generator] #WegnerGenerator(ha, ot)
fl = Flow_IMSRG2(ha, ot)
#wg_spin = WegnerGenerator(ss, ot)
# fl_spin = Flow_IMSRG2(ss, ot)
initf = time.time() # finish instantiation timer
if verbose:
......@@ -218,44 +239,168 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1, flow_dat
flowi = time.time()
# --- Solve the IM-SRG flow
#y0 = np.concatenate([unravel(ha.E, ha.f, ha.G), unravel(ss.E, ss.f, ss.G)], axis=0)
y0 = unravel(ha.E, ha.f, ha.G)
H0B, H1B, H2B = get_vacuum_coeffs(ha.E, ha.f, ha.G, ha.sp_basis, ha.holes)
zero, eta1B_vac, eta2B_vac = get_vacuum_coeffs(0.0, wg.eta1B, wg.eta2B, ha.sp_basis, ha.holes)
pickle.dump( (H0B, H1B, H2B, eta1B_vac, eta2B_vac), open( output_root+"/vac_coeffs_unevolved.p", "wb" ) )
solver = ode(derivative,jac=None)
solver.set_integrator('vode', method='bdf', order=5, nsteps=1000)
solver.set_f_params(ha, ot, wg, fl)
solver.set_f_params([ha, ot, wg, fl])
solver.set_initial_value(y0, 0.)
sfinal = 50
# ds = 1
ds = 0.01
s_vals = []
E_vals = []
E_gs_lst = []
s_expect_lst = []
iters = 0
convergence = 0
fci_hme = pyci.matrix(n_holes, n_particles, 0.0, d, g, pb)
dens_weights,v = np.linalg.eigh(fci_hme)
if verbose:
print("iter, \t s, \t E, \t ||eta1B||, \t ||eta2B||")
print_columns = ['iter',
's',
'E',
'CI_gs',
'||eta1b||',
'||eta2b||']
column_string = '{: >6}, '
for i, string in enumerate(print_columns[1::]):
if i != len(print_columns)-2:
column_string += '{: >'+str(11)+'}, '
else:
column_string += '{: >'+str(11)+'}'
print(column_string.format(*print_columns))
eig_idx = 0
while solver.successful() and solver.t < sfinal:
#print('solver success,', solver.successful())
ys = solver.integrate(sfinal, step=True)
ys = solver.integrate(solver.t+ds, step=True)
Es, fs, Gs = ravel(ys, ha.n_sp_states)
# E_spins, f_spins, G_spins = ravel(ys[half::], ss.n_sp_states)
s_vals.append(solver.t)
E_vals.append(Es)
# compute density matrices, get expectation values
# H0B, H1B, H2B = get_vacuum_coeffs(Es, fs, Gs, ha.sp_basis, ha.holes)
# SS0B, SS1B, SS2B = get_vacuum_coeffs(E_spins, f_spins, G_spins, ss.sp_basis, ss.holes)
# hme = pyci.matrix(n_holes, n_particles, H0B, H1B, H2B, H2B, imsrg=True)
# w,v = np.linalg.eigh(hme)
# v0 = v[:,eig_idx]
# rho1b = density_1b(n_holes,n_particles, weights=v0)
# rho2b = density_2b(n_holes,n_particles, weights=v0)
# contract_1b = np.einsum('ij,ij', rho1b, SS1B)
# rho_reshape_2b = np.reshape(rho2b, (ss.n_sp_states**2,ss.n_sp_states**2))
# h2b_reshape_2b = np.reshape(SS2B, (ss.n_sp_states**2,ss.n_sp_states**2))
# contract_2b = np.einsum('ij,ij', rho_reshape_2b, h2b_reshape_2b)
# s_expect = SS0B + contract_1b + 0.25*contract_2b
# E_gs_lst.append(w[eig_idx])
# s_expect_lst.append(s_expect)
if iters %10 == 0 and verbose:
zero, eta1Bv, eta2Bv = get_vacuum_coeffs(0.0,wg.eta1B,wg.eta2B,ha.sp_basis,ha.holes)
# ggme = pyci.matrix(n_holes, n_particles, zero, eta1Bv, eta2Bv, eta2Bv, imsrg=True)
# ssme = pyci.matrix(n_holes, n_particles, SS0B, SS1B, SS2B, SS2B, imsrg=True)
norm_eta1B = np.linalg.norm(np.ravel(wg.eta1B))
norm_eta2B = np.linalg.norm(np.ravel(wg.eta2B))
print("{:>6d}, \t {:0.4f}, \t {:0.8f}, \t {:0.8f}, \t {:0.8f}".format(iters,
solver.t,
Es,
norm_eta1B,
norm_eta2B))
num_sp = n_holes+n_particles
axes = (num_sp**2,num_sp**2)
H0B, H1B, H2B = get_vacuum_coeffs(Es, fs, Gs, ha.sp_basis, ha.holes)
hme = pyci.matrix(n_holes, n_particles, H0B, H1B, H2B, H2B, imsrg=True)
w,v = np.linalg.eigh(hme)
# SS2B_r = np.reshape(SS2B, (num_sp**2,num_sp**2))
# H2B_r = np.reshape(H2B, (num_sp**2, num_sp**2))
# commute1b = np.linalg.norm(SS1B.dot(H1B) - H1B.dot(SS1B))
# commute2b = np.linalg.norm(SS2B_r.dot(H2B_r) - H2B_r.dot(SS2B_r))
# hfd,hfod,hGd,hGod = wg.decouple_OD()
# sfd,sfod,sGd,sGod = wg_spin.decouple_OD()
# hfd,hfod,hGd,hGod = decouple_od(fs, Gs, ha.particles, ha.holes)
# sfd,sfod,sGd,sGod = decouple_od(f_spins, G_spins, ha.particles, ha.holes)
# hGd,hGod,sGd,sGod = np.reshape(hGd,axes),np.reshape(hGod,axes),np.reshape(sGd,axes),np.reshape(sGod,axes)
# commute1bd = np.linalg.norm(hfd.dot(sfd) - sfd.dot(hfd))
# commute1bod = np.linalg.norm(hfod.dot(sfod) - sfod.dot(hfod))
# commute2bd = np.linalg.norm(hGd.dot(sGd) - sGd.dot(hGd))
# commute2bod = np.linalg.norm(hGod.dot(sGod) - sGod.dot(hGod))
# commute = np.linalg.norm(hme.dot(ssme) - ssme.dot(hme))
# commutePairBlock = np.linalg.norm(hme[0:2,0:2].dot(ssme[0:2,0:2]) - ssme[0:2,0:2].dot(hme[0:2,0:2]))
data_columns = [iters,
solver.t,
Es,
w[eig_idx],
norm_eta1B,
norm_eta2B]
column_string = '{:>6d}, '
for i, string in enumerate(print_columns[1::]):
if i != len(print_columns)-2:
column_string += '{: .8f}, '
else:
column_string += '{: .8f}'
# print(column_string.format(*data_columns))
print(column_string.format(*data_columns))
#print(v0*v0)
# ws,vs = np.linalg.eigh(ssme)
# with np.printoptions(linewidth=999,precision=4):
# print(ssme)
# print(ggme)
# for i in range(wg.eta1B.shape[0]):
# print("{d}".format(d=wg.eta1B[i,:]))
# for i in range(wg.eta1B.shape[0]):
# print("{d}".format(d=f_spins[i,:]))
# for i in range(ssme[0:6,0:6].shape[0]):
# print("{d}".format(d=ssme[0:6,0:6][i,:]))
#print(np.reshape(wg.eta2B,axes))
#print(ssme[0:6,0:6])A
#"{:>6d}, {: .8f}, {: .8f}, {: .8f}, {: .8f}, {: .8f}, {: .8f}, {: .8f}, {: .8f}, {: .8f}, {: .8f}, {: .8f}, {: .8f}, {: .8f}, {: .8f}, {: .8f}".format(*data_columns))
# solver.t,
# Es,
# w[0],
# E_spins,
# s_expect,
# SS0B,
# contract_1b,
# contract_2b,
# norm_eta1B,
# norm_eta2B,
# 0.0,
# 0.0,
# 0.0,
# 0.0))
# commute1bd,
# commute1bod,
# commute2bd,
# commute2bod))
if flow_data_log and iters %10 == 0:
H0B, H1B, H2B = get_vacuum_coeffs(Es, fs, Gs, ha.sp_basis, ha.holes)
......@@ -277,7 +422,7 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1, flow_dat
# if verbose: print("---- Energy diverged at iter {:>06d} with energy {:3.8f}\n".format(iters,E_vals[-1]))
# break
if iters > 10000:
if iters > 300:
if verbose: print("---- Energy diverged at iter {:>06d} with energy {:3.8f}\n".format(iters,E_vals[-1]))
break
......@@ -300,7 +445,7 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1, flow_dat
num_sp = n_holes+n_particles
del ha, ot, wg, fl, solver, y0, sfinal
return (convergence, iters, d, g, pb, num_sp, s_vals, E_vals, time_str)
......@@ -352,14 +497,63 @@ if __name__ == '__main__':
# ----------------------------------------------------------------
refs = [[1,1,1,1,0,0,0,0],[1,1,0,0,1,1,0,0],[1,1,0,0,0,0,1,1],
[0,0,1,1,1,1,0,0],[0,0,1,1,0,0,1,1]]
# main(4,4, generator='white')
# data = pickle.load(open('expect_flow.p', 'rb'))
# fig = plt.figure(figsize=(8,4))
# sns.lineplot(x='s', y=data['E_gs']/data['E_gs'][0], data=data)
# sns.lineplot(x='s', y=data['s_expect']/data['s_expect'][0], data=data)
# plt.legend(['E(s)/E(s=0)', 'SS(s)/SS(s=0)'])
# plt.savefig('flow_conservation.png')
#0,3,14,15,28,35
ref = 0.5*np.asarray(refs[0]) + (1.0-0.5)/4.0*(np.asarray(refs[1]) + np.asarray(refs[2]) + np.asarray(refs[3]) + np.asarray(refs[4]))
# main(4,4, g=5, ref=[1,1,1,1,0,0,0,0])
g = 0.5
pb = 0.1
hme = pyci.matrix(4,4,0.0,1.0,g,pb)
w,v = np.linalg.eigh(hme)
v0 = v[:,5]
#ref = 0.7*np.array([1,1,1,1,0,0,0,0])+0.3*np.array([1,1,0,0,1,1,0,0])
basis = pyci.gen_basis(4,4)[:,1::]
#ref = 0.2*basis[0,:]+ 0.8*basis[6, :]
#ref = basis.T.dot(v0**2)
#ref = 0.8*basis[0,:] + 0.2*basis[1,:]
#ref = basis.T.dot(v0*v0)
ref = basis[0,:]
# ref_input = sys.argv[1]
# ref = [int(x) for x in list(ref_input)]
#ref = pickle.load(open('reference_g2.00_pb0.01_4-4.p', 'rb'))
main(4,4, g=g, ref=ref, pb=pb, generator='wegner')
print('FCI ENERGY = {: .8f}'.format(w[0]))
data = pickle.load(open('expect_flow.p', 'rb'))
fig = plt.figure(figsize=(8,4))
sns.lineplot(x='s', y=data['E_gs']-data['E_gs'][0], data=data)
sns.lineplot(x='s', y=data['s_expect']-data['s_expect'][0], data=data)
plt.legend(['E(s)/E(s=0)', 'SS(s)/SS(s=0)'])
plt.savefig('flow_conservation.png')
hme = pyci.matrix(4,4,0.0,1.0,0.5,0.0)
w,v = np.linalg.eigh(hme)
fig = plt.figure(figsize=(8,4))
sns.lineplot(x='s', y=data['E_gs'], data=data)
sns.lineplot(x='s', y=w[0], data=data)
plt.legend(['E_gs evolution', 'FCI gs'])
plt.savefig('E_gs_error_g2.png')
# refs = [[1,1,1,1,0,0,0,0],[1,1,0,0,1,1,0,0],[1,1,0,0,0,0,1,1],
# [0,0,1,1,1,1,0,0],[0,0,1,1,0,0,1,1]]
ref = 0.8*np.array([1,1,1,1,0,0,0,0])+0.2*np.array([1,1,0,0,1,1,0,0])
main(4,4,g=2, ref=ref, pb=0.0, generator='white')
# ref = 0.5*np.asarray(refs[0]) + (1.0-0.5)/4.0*(np.asarray(refs[1]) + np.asarray(refs[2]) + np.asarray(refs[3]) + np.asarray(refs[4]))
# # main(4,4, g=5, ref=[1,1,1,1,0,0,0,0])
# H1B_true, H2B_true = pickle.load(open('comparison.p','rb'))
# H1B, H2B = pickle.load(open('vac_coeffs_unevolved.p', 'rb'))
......
......@@ -115,9 +115,9 @@ def main3b(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
initi = time.time() # start instantiation timer
if ref == None:
if ref is None:
ha = PairingHamiltonian2B(n_holes, n_particles, d=d, g=g, pb=pb)
ref = [1,1,1,1,0,0,0,0] # this is just for printing
ref = ha.reference # this is just for printing
else:
ha = PairingHamiltonian2B(n_holes, n_particles, ref=ref, d=d, g=g, pb=pb)
print("Built Hamiltonian")
......@@ -138,7 +138,7 @@ def main3b(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
snapshot = tracemalloc.take_snapshot()
top_stats = snapshot.statistics('lineno')
total = sum(stat.size for stat in top_stats)
print("Total allocated size: %.1f GB" % (total / 1024**3))
print("Total allocated size: %.1f MB" % (total / 1024**2))
print("""Pairing model IM-SRG(3) flow:
......@@ -163,15 +163,31 @@ def main3b(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
solver.set_initial_value(y0, 0.)
sfinal = 50
ds = 0.1
ds = 0.01
s_vals = []
E_vals = []
iters = 0
convergence = 0
print_columns = ['iter',
's',
'E',
'||eta1b||',
'||eta2b||',
'||eta3b||']
column_string = '{: >6}, '
for i, string in enumerate(print_columns[1::]):
if i != len(print_columns)-2:
column_string += '{: >'+str(11)+'}, '
else:
column_string += '{: >'+str(11)+'}'
print(column_string.format(*print_columns))
while solver.successful() and solver.t < sfinal:
ys = solver.integrate(sfinal, step=True)
ys = solver.integrate(solver.t+ds, step=True)
Es, fs, Gs, Ws = ravel(ys, ha.n_sp_states)
s_vals.append(solver.t)
E_vals.append(Es)
......@@ -179,9 +195,34 @@ def main3b(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
iters += 1
# if iters == 1:
# break
if iters %10 == 0: print("iter: {:>6d} \t scale param: {:0.4f} \t E = {:0.9f}".format(iters, solver.t, Es))
if len(E_vals) > 100 and abs(E_vals[-1] - E_vals[-2]) < 10**-8 and E_vals[-1] != E_vals[0]:
if iters %10 == 0:
# print("iter: {:>6d} \t scale param: {:0.4f} \t E = {:0.9f}".format(iters, solver.t, Es))
norm_eta1B = np.linalg.norm(np.ravel(wg.eta1B))
norm_eta2B = np.linalg.norm(np.ravel(wg.eta2B))
norm_eta3B = np.linalg.norm(np.ravel(wg.eta3B))
data_columns = [iters,
solver.t,
Es,
norm_eta1B,
norm_eta2B,
norm_eta3B]
column_string = '{:>6d}, '
for i, string in enumerate(print_columns[1::]):
if i != len(print_columns)-2:
column_string += '{: .8f}, '
else:
column_string += '{: .8f}'
# print(column_string.format(*data_columns))
print(column_string.format(*data_columns))
if len(E_vals) > 100 and abs(E_vals[-1] - E_vals[-2]) < 10**-8:
print("---- Energy converged at iter {:>06d} with energy {:1.8f}\n".format(iters,E_vals[-1]))
convergence = 1
break
......@@ -259,7 +300,7 @@ def test_exact(plots_dir):
if __name__ == '__main__':
#test_exact("plots3b/")
main3b(4,4)
main3b(2,2,g=0.5,pb=0.1)
#test = main3b(4,4)
#tracemalloc.start()
......
This diff is collapsed.
......@@ -26,7 +26,7 @@ class Flow_IMSRG2(Flow):
h -- Hamiltonian object
occ_t -- OccupationTensors object"""
assert isinstance(h, Hamiltonian), "Arg 0 must be PairingHamiltonian object"
#assert isinstance(h, Hamiltonian), "Arg 0 must be PairingHamiltonian object"
assert isinstance(occ_t, OccupationTensors), "Arg 1 must be OccupationTensors object"
# self.f = h.f
......@@ -58,7 +58,7 @@ class Flow_IMSRG2(Flow):
# def G(self, G):
# self._G = G
def flow(self, gen):
def flow(self, f, G, gen):
"""Iterates the IMSRG2 flow equations once.
Arugments:
......@@ -75,8 +75,8 @@ class Flow_IMSRG2(Flow):
# f = self.f
# G = self.G
f = gen.f.astype(np.float32)
G = gen.G.astype(np.float32)
# f = gen.f.astype(np.float32)
# G = gen.G.astype(np.float32)
eta1B, eta2B = gen.calc_eta()
# eta1B = partition[0]#.astype(np.float32)
......@@ -233,9 +233,11 @@ class Flow_IMSRG3(Flow_IMSRG2):
self._particles = h.particles
self._occA = occ_t.occA
self._occA2 = occ_t.occA2
self._occA4 = occ_t.occA4
self._occB = occ_t.occB
self._occB4 = occ_t.occB4
self._occC = occ_t.occC
self._occC6 = occ_t.occC6
self._occD = occ_t.occD
self._occE = occ_t.occE
self._occF = occ_t.occF
......@@ -244,6 +246,9 @@ class Flow_IMSRG3(Flow_IMSRG2):
self._occI = occ_t.occI
self._occJ = occ_t.occJ
self._h = h
self._occ_t = occ_t
def flow(self, gen):
"""Iterates the IMSRG3 flow equations once. Extends IMSRG2 flow function.
......@@ -266,10 +271,15 @@ class Flow_IMSRG3(Flow_IMSRG2):
G = gen.G
W = gen.W
partition = super().flow(gen)
gen2b = WegnerGenerator(self._h, self._occ_t)
gen2b.f = f
gen2b.G = G
partition = super().flow(f,G,gen2b)
dE = partition[0]
df = partition[1]
dG = partition[2]
del gen2b
partition = gen.calc_eta()
eta1B = partition[0]
......@@ -277,9 +287,11 @@ class Flow_IMSRG3(Flow_IMSRG2):
eta3B = partition[2]
occA = self._occA
occA2 = self._occA2
occA4 = self._occA4
occB = self._occB
occB4 = self._occB4
occC = self._occC
occC6 = self._occC6
occD = self._occD
occE = self._occE
occF = self._occF
......@@ -289,26 +301,34 @@ class Flow_IMSRG3(Flow_IMSRG2):
occJ = self._occJ
# Calculate 0B flow equation
sum3_0b_1 = np.matmul(eta3B, occE)
sum3_0b_1 = np.multiply(eta3B, occE)
sum3_0b = tn.ncon([sum3_0b_1, W], [(1,2,3,4,5,6), (4,5,6,1,2,3)])#.numpy()
dE += (1/18)*sum3_0b
# Calculate 1B flow equation
# fourth term
sum4_1b_1 = np.matmul(np.transpose(occD,[2,3,0,1]), G)
sum4_1b_2 = np.matmul(np.transpose(occD,[2,3,0,1]), eta2B)
sum4_1b_1 = np.multiply(np.transpose(occD.tensor,[2,3,0,1]), G)
sum4_1b_2 = np.multiply(np.transpose(occD.tensor,[2,3,0,1]), eta2B)
sum4_1b_3 = tn.ncon([eta3B, sum4_1b_1], [(1,2,-1,3,4,-2),(3,4,1,2)])#.numpy()
sum4_1b_4 = tn.ncon([W, sum4_1b_2], [(1,2,-1,3,4,-2),(3,4,1,2)])#.numpy()
sum4_1b_4 = tn.ncon([W, sum4_1b_2], [(1,2,-1,3,4,-2),(3,4,1,2)])#.numpy()
sum4_1b = sum4_1b_3 - sum4_1b_4
# fifth term
sum5_1b_1 = tn.ncon([eta3B, occF, W], [(6,7,-1,8,9,10),
(6,7,8,9,10,1,2,3,4,5),
(3,4,5,1,2,-2)])#.numpy()
sum5_1b_2 = tn.ncon([W, occF, eta3B], [(6,7,-1,8,9,10),
(6,7,8,9,10,1,2,3,4,5),
(3,4,5,1,2,-2)])#.numpy()
sum5_1b = sum5_1b_1 - sum5_1b_2
# sum5_1b_1 = tn.ncon([eta3B, occF, W], [(6,7,-1,8,9,10),
# (6,7,8,9,10,1,2,3,4,5),
# (3,4,5,1,2,-2)])#.numpy()
# sum5_1b_2 = tn.ncon([W, occF, eta3B], [(6,7,-1,8,9,10),
# (6,7,8,9,10,1,2,3,4,5),
# (3,4,5,1,2,-2)])#.numpy()
# sum5_1b = sum5_1b_1 - sum5_1b_2