Commit 64583bae authored by Jacob Davison's avatar Jacob Davison
Browse files

updated code logic w/r to IMSRG3

parent cf40e5a8
......@@ -28,6 +28,7 @@ from oop_imsrg.flow import *
from oop_imsrg.plot_data import *
# from oop_imsrg.display_memory import *
import oop_imsrg.ci_pairing.cipy_pairing_plus_ph as ci_matrix
from oop_imsrg.tests2B import *
# @profile
def derivative(t, y, hamiltonian, occ_tensors, generator, flow):
......@@ -161,29 +162,29 @@ def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
ys = solver.integrate(sfinal, step=True)
Es, fs, Gs = ravel(ys, ha.n_sp_states)
#s_vals.append(solver.t)
#E_vals.append(Es)
s_vals.append(solver.t)
E_vals.append(Es)
iters += 1
if iters == 176:
break
# if iters %10 == 0: print("iter: {:>6d} \t scale param: {:0.4f} \t E = {:0.9f}".format(iters, solver.t, Es))
# if iters == 176:
# break
if iters %10 == 0: print("iter: {:>6d} \t scale param: {:0.4f} \t E = {:0.8f}".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]:
# print("---- Energy converged at iter {:>06d} with energy {:1.8f}\n".format(iters,E_vals[-1]))
# convergence = 1
# break
if len(E_vals) > 100 and abs(E_vals[-1] - E_vals[-2]) < 10**-8 and E_vals[-1] != E_vals[0]:
print("---- Energy converged at iter {:>06d} with energy {:1.8f}\n".format(iters,E_vals[-1]))
convergence = 1
break
# if len(E_vals) > 100 and abs(E_vals[-1] - E_vals[-2]) > 1:
# print("---- Energy diverged at iter {:>06d} with energy {:3.8f}\n".format(iters,E_vals[-1]))
# break
if len(E_vals) > 100 and abs(E_vals[-1] - E_vals[-2]) > 1:
print("---- Energy diverged at iter {:>06d} with energy {:3.8f}\n".format(iters,E_vals[-1]))
break
# if iters > 20000:
# print("---- Energy diverged at iter {:>06d} with energy {:3.8f}\n".format(iters,E_vals[-1]))
# break
if iters > 20000:
print("---- Energy diverged at iter {:>06d} with energy {:3.8f}\n".format(iters,E_vals[-1]))
break
# if iters % 1000 == 0:
# print('Iteration {:>06d}'.format(iters))
if iters % 1000 == 0:
print('Iteration {:>06d}'.format(iters))
flowf = time.time()
end = time.time()
time_str = "{:2.5f}\n".format(end-start)
......@@ -219,226 +220,12 @@ def exact_diagonalization(d, g):
return E
# @profile
def scan_params():
tracemalloc.start()
log_dir = "logs\\"
plot_dir = "plots\\"
if not os.path.exists(log_dir):
os.mkdir(log_dir)
print("Created log directory at "+log_dir)
if not os.path.exists(plot_dir):
os.mkdir(plot_dir)
print("Created plots directory at "+plot_dir+"\n")
if os.path.isfile(log_dir+'total_mem.txt'):
os.remove(log_dir+'total_mem.txt')
start = 0.0001
stop = 5
num = 200
gsv = np.linspace(start, stop, num)
pbv = np.copy(gsv)
# gsv = np.append(np.linspace(-stop,-start,num), np.linspace(start, stop,num))
# pbs = np.copy(gsv)
# data_container = np.array([])
for g in gsv:
pb_list = np.array(['convergence', 'iters', 'd', 'g', 'pb', 'n_sp_states', 's_vals', 'E_vals', 'time_str'])
for pb in pbv:
data = main(4,4, d=stop, g=g, pb=pb) # (convergence, s_vals, E_vals)
pb_list = np.vstack([pb_list, data])
snapshot = tracemalloc.take_snapshot()
top_stats = snapshot.statistics('lineno')
total_mem = sum(stat.size for stat in top_stats)
with open(log_dir+'total_mem.txt', 'a') as f:
f.write('{:.1f}\n'.format(total_mem))
if data[0] == 0:
print("Energy diverged. Continuing to next g value...\n")
break
del data, snapshot, top_stats, total_mem
with open('{:s}g-{:2.4f}.pickle'.format(log_dir,g), 'wb') as f:
pickle.dump(pb_list, f, pickle.HIGHEST_PROTOCOL)
plot_data(log_dir, plot_dir)
del pb_list # delete resources that have been written
# data_container = np.append(data_container, pb_list)
def test_exact(plots_dir):
assert isinstance(plots_dir, str), "Enter plots directory as string"
assert (plots_dir[-1] == '\\' or
plots_dir[-1] == '/'), "Directory must end in slash (\\ or /)"
if not os.path.exists(plots_dir):
os.mkdir(plots_dir)
start = -1.0
stop = 1.0
num = 21
g_vals = np.linspace(start, stop, num)
for pb in g_vals:
E_corrs = []
E_exacts = []
for g in g_vals:
data = main(4,4, d=1.0, g=g, pb=0.0)
E_vals = data[7]
E_corr = E_vals[-1]
E_exact = exact_diagonalization(1.0, g)
E_corrs.append(E_corr - (2-g))
E_exacts.append(E_exact - (2-g))
plt.figure(figsize=[12,8])
plt.plot(data[6], data[7])
plt.ylabel('Energy')
plt.xlabel('scale parameter')
plt.title('Convergence for \n g={:2.4f}, pb={:2.4f}'.format(g,pb))
pb_plots_dir = plots_dir+'pb{:2.4f}\\'.format(pb)
if not os.path.exists(pb_plots_dir):
os.mkdir(pb_plots_dir)
plt.savefig(pb_plots_dir+'g{:2.4f}_pb{:2.4f}.png'.format(g,pb))
plt.close()
plt.figure(figsize=[12,8])
plt.plot(g_vals, E_exacts, marker='s')
plt.plot(g_vals, E_corrs, marker='v')
plt.ylabel('E$_{corr}$')
plt.xlabel('g')
plt.legend(['exact', 'IMSRG(2)'])
plt.title('Correlation energy with pb = {:2.4f}'.format(pb))
plt.savefig(plots_dir+'pb{:2.4f}.png'.format(pb))
plt.close()
print(E_exacts)
break
def test_refs(plots_dir):
assert isinstance(plots_dir, str), "Enter plots directory as string"
assert (plots_dir[-1] == '\\' or
plots_dir[-1] == '/'), "Directory must end in slash (\\ or /)"
if not os.path.exists(plots_dir):
os.mkdir(plots_dir)
start = -1.0
stop = 1.0
num = 5
g_vals = np.linspace(start, stop, num)
refs = list(map("".join, itertools.permutations('11110000')))
refs = list(dict.fromkeys(refs)) # remove duplicates
refs = [list(map(int, list(ref))) for ref in refs]
fig_corr = plt.figure(figsize=[12,8])
fig_conv = plt.figure(figsize=[12,8])
ax_corr = fig_corr.add_subplot(1,1,1)
ax_conv = fig_conv.add_subplot(1,1,1)
for pb in g_vals:
E_exacts = []
for g in g_vals:
E_corrs = []
# plt.figure(figsize=[12,8])
ref_rand = random.sample(refs, 20)
E_exact = ci_matrix.exact_diagonalization(1.0, g, pb)
# E_exacts.append(E_exact - (2-g))
E_exacts.append(E_exact)
refs_conv = []
for ref in refs:
data = main(4,4, ref=ref, d=1.0, g=g, pb=pb)
if data[0] == 1:
E_vals = data[7]
E_corr = E_vals[-1]
# E_corrs.append(E_corr - (2-g))
E_corrs.append(E_corr)
refs_conv.append(ref)
ax_conv.plot(data[6], data[7])
# ymin, ymax = ax_conv.get_ylim()
# ax_conv.set_ylim(bottom=0.7*ymin, top=0.7*ymax)
ax_conv.set_ylabel('Energy')
ax_conv.set_xlabel('scale parameter')
ax_conv.set_title('Convergence for \n g={:2.4f}, pb={:2.4f}'.format(g,pb))
ax_conv.legend(refs_conv)
pb_plots_dir = plots_dir+'pb{:2.4f}\\'.format(pb)
if not os.path.exists(pb_plots_dir):
os.mkdir(pb_plots_dir)
fig_conv.savefig(pb_plots_dir+'g{:2.4f}_pb{:2.4f}.png'.format(g,pb))
ax_conv.clear()
with open(plots_dir+'g{:2.4f}_pb{:2.4f}.txt'.format(g,pb), 'w') as f:
f.write('Pairing model: d = 1.0, g = {:2.4f}, pb = {:2.4f}\n'.format(g, pb))
f.write('Full CI diagonalization -- correlation energy: {:2.8f}\n'.format(E_exacts[-1]))
f.write('IMSRG(2) variable reference state -- correlation energies:\n')
if len(refs_conv) == 0:
f.write('No convergence for range of reference states tested\n')
else:
for i in range(len(refs_conv)):
f.write('{:2.8f} | {d}\n'.format(E_corrs[i], d=refs_conv[i]))
f.write('Ground state from IMSRG(2):\n')
e_sort_ind = np.argsort(E_corrs)
print(e_sort_ind)
print(E_corrs)
f.write('{:2.8f} | {d}\n'.format(E_corrs[e_sort_ind[0]], d=refs_conv[e_sort_ind[0]]))
# corr_data = np.reshape(E_corrs, (len(g_vals), 2))
# ax_corr.plot(g_vals, E_exacts, marker='s')
# for i in range(10):
# ax_corr.plot(g_vals, corr_data[:,i], marker='v')
# ymin, ymax = ax_corr.get_ylim()
# ax_corr.set_ylim(bottom=0.7*ymin, top=0.7*ymax)
# ax_corr.set_ylabel('E$_{corr}$')
# ax_corr.set_xlabel('g')
# ax_corr.legend(['exact', 'IMSRG(2)'])
# ax_corr.set_title('Correlation energy with pb = {:2.4f}'.format(pb))
# fig_corr.savefig(plots_dir+'pb{:2.4f}.png'.format(pb))
# ax_corr.clear()
# def search_refs(log_dir):
# refs = list(map("".join, itertools.permutations('11110000')))
# refs = list(dict.fromkeys(refs)) # remove duplicates
# refs = [list(map(int, list(ref))) for ref in refs]
# # refs = random.sample(refs, 2)
# E_conv = []
# refs_conv = []
# for ref in refs:
# data = main(4,4, ref=ref, g=0.5, pb=0.1)
#
# if data[0] == 1: # energy converged
# Es = data[7]
# E_conv.append(Es[-1])
# refs_conv.append(ref)
#
# for i in range(len(E_conv)):
# print('{:2.4f} | {d}'.format(E_conv[i], d=refs_conv[i]))
if __name__ == '__main__':
#test_refs('logs_refs\\')
# test = main(4,4)
print(ci_matrix.exact_diagonalization(1.0, 0.5,0.0))
test = main(4,4)
# # test_exact('plots_exact\\')
# # print(ci_matrix.exact_diagonalization(1.0, 0.5,0.0))
# print(exact_diagonalization(1.0,0.5))
......@@ -158,7 +158,7 @@ def exact_diagonalization(d, g, pb):
w,v = np.linalg.eig(hme)
w_sort = w[np.argsort(w)]
print("--- Ground state energy: {:2.4f}\n".format(w_sort[0]))
print("--- Ground state energy: {:2.8f}\n".format(w_sort[0]))
return w_sort[0] # ground state energy
#
# hme = matrix(states, 1.0, 0.5, 0.1) # delta, g, f
......
......@@ -72,12 +72,12 @@ class Flow_IMSRG2(Flow):
# f = self.f
# G = self.G
f = gen.f
G = gen.G
f = gen.f.astype(np.float32)
G = gen.G.astype(np.float32)
partition = gen.calc_eta()
eta1B = partition[0]
eta2B = partition[1]
eta1B = partition[0].astype(np.float32)
eta2B = partition[1].astype(np.float32)
# occA = occ_t.occA
# occB = occ_t.occB
......@@ -260,15 +260,23 @@ class Flow_IMSRG3(Flow_IMSRG2):
sum4_1b = sum4_1b_3 - sum4_1b_4
# fifth term
sum5_1b_1 = ncon([occF, eta3B.astype(np.float32)],
[(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)]).numpy()
sum5_1b_2 = ncon([occF, W.astype(np.float32)],
[(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)]).numpy()
sum5_1b_3 = ncon([sum5_1b_1, W.astype(np.float32)],
[(0,1,-1,2,3,4), (2,3,4,0,1,-2)]).numpy()
sum5_1b_4 = ncon([sum5_1b_2, eta3B.astype(np.float32)],
[(0,1,-1,2,3,4), (2,3,4,0,1,-2)]).numpy()
sum5_1b = sum5_1b_3 - sum5_1b_4
sum5_1b_1 = ncon([eta3B, occF, W], [(5,6,-1,7,8,9),
(5,6,7,8,9,0,1,2,3,4),
(2,3,4,0,1,-2)]).numpy()
sum5_1b_2 = ncon([W, occF, eta3B], [(5,6,-1,7,8,9),
(5,6,7,8,9,0,1,2,3,4),
(2,3,4,0,1,-2)]).numpy()
sum5_1b = sum5_1b_1 - sum5_1b_2
# sum5_1b_1 = ncon([occF, eta3B.astype(np.float32)],
# [(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)]).numpy()
# sum5_1b_2 = ncon([occF, W.astype(np.float32)],
# [(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)]).numpy()
# sum5_1b_3 = ncon([sum5_1b_1, W.astype(np.float32)],
# [(0,1,-1,2,3,4), (2,3,4,0,1,-2)]).numpy()
# sum5_1b_4 = ncon([sum5_1b_2, eta3B.astype(np.float32)],
# [(0,1,-1,2,3,4), (2,3,4,0,1,-2)]).numpy()
# sum5_1b = sum5_1b_3 - sum5_1b_4
df += 0.25*sum4_1b + (1/12)*sum5_1b
......@@ -281,28 +289,53 @@ class Flow_IMSRG3(Flow_IMSRG2):
sum4_2b = sum4_2b_3 - sum4_2b_4
#fifth term
sum5_2b_1 = ncon([occG, G], [(-1,-2,-4,0,1,2), (1,2,0,-3)]).numpy()
sum5_2b_2 = ncon([occG, eta2B], [(-1,-2,-4,0,1,2), (1,2,0,-3)]).numpy()
sum5_2b_3 = ncon([eta3B, sum5_2b_1], [(0,-1,-2,1,2,-4), (1,2,0,-3)]).numpy()
sum5_2b_4 = ncon([W, sum5_2b_2], [(0,-1,-2,1,2,-4), (1,2,0,-3)]).numpy()
sum5_2b_5 = sum5_2b_3 - sum5_2b_4
sum5_2b = sum5_2b_5 - np.transpose(sum5_2b_5, [3,2,0,1]) - \
np.transpose(sum5_2b_5, [0,1,3,2]) + \
np.transpose(sum5_2b_5, [2,3,0,1])
sum5_2b_1 = ncon([eta3B, occG, G], [(3,-1,-2,4,5,-4),
(3,4,5,0,1,2),
(1,2,0,-3)]).numpy()
sum5_2b_2 = ncon([W, occG, eta2B], [(3,-1,-2,4,5,-4),
(3,4,5,0,1,2),
(1,2,0,-3)]).numpy()
sum5_2b = sum5_2b_2 - np.transpose(sum5_2b_2, [3,2,0,1]) - \
np.transpose(sum5_2b_2, [0,1,3,2]) + \
np.transpose(sum5_2b_2, [2,3,0,1])
# sum5_2b_1 = ncon([occG, G], [(-1,-2,-4,0,1,2), (1,2,0,-3)]).numpy()
# sum5_2b_2 = ncon([occG, eta2B], [(-1,-2,-4,0,1,2), (1,2,0,-3)]).numpy()
# sum5_2b_3 = ncon([eta3B, sum5_2b_1], [(0,-1,-2,1,2,-4), (1,2,0,-3)]).numpy()
# sum5_2b_4 = ncon([W, sum5_2b_2], [(0,-1,-2,1,2,-4), (1,2,0,-3)]).numpy()
# sum5_2b_5 = sum5_2b_3 - sum5_2b_4
# sum5_2b = sum5_2b_5 - np.transpose(sum5_2b_5, [3,2,0,1]) - \
# np.transpose(sum5_2b_5, [0,1,3,2]) + \
# np.transpose(sum5_2b_5, [2,3,0,1])
#sixth term
sum6_2b_1 = ncon([occH, W], [(-1,-2,-3,-4,0,1,2,3),(1,2,3,0,-5,-6)]).numpy()
sum6_2b_2 = ncon([occH, W], [(-3,-4,-5,-6,0,1,2,3),(0,-1,-2,1,2,3)]).numpy()
sum6_2b_3 = ncon([eta3B, sum6_2b_1], [(0,-1,-2,1,2,3), (1,2,3,0,-3,-4)]).numpy()
sum6_2b_4 = ncon([eta3B, sum6_2b_2], [(1,2,3,0,-3,-4), (0,-1,-2,1,2,3)]).numpy()
sum6_2b = sum6_2b_3 - sum6_2b_4
sum6_2b_1 = ncon([eta3B, occH, W], [(4,-1,-2,5,6,7),
(4,5,6,7,0,1,2,3),
(1,2,3,0,-3,-4)]).numpy()
sum6_2b_2 = ncon([eta3B, occH, W], [(5,6,7,4,-3,-4),
(4,5,6,7,0,1,2,3),
(0,-1,-2,1,2,3)]).numpy()
sum6_2b = sum6_2b_1 - sum6_2b_2
# sum6_2b_1 = ncon([occH, W], [(-1,-2,-3,-4,0,1,2,3),(1,2,3,0,-5,-6)]).numpy()
# sum6_2b_2 = ncon([occH, W], [(-3,-4,-5,-6,0,1,2,3),(0,-1,-2,1,2,3)]).numpy()
# sum6_2b_3 = ncon([eta3B, sum6_2b_1], [(0,-1,-2,1,2,3), (1,2,3,0,-3,-4)]).numpy()
# sum6_2b_4 = ncon([eta3B, sum6_2b_2], [(1,2,3,0,-3,-4), (0,-1,-2,1,2,3)]).numpy()
# sum6_2b = sum6_2b_3 - sum6_2b_4
#seventh term
sum7_2b_1 = ncon([occI, W], [(-1,-2,-3,-4,0,1,2,3), (2,3,-5,0,1,-6)]).numpy()
sum7_2b_2 = ncon([eta3B, sum7_2b_1], [(0,1,-1,2,3,-4),(2,3,-2,0,1,-3)]).numpy()
sum7_2b = sum7_2b_2 - np.transpose(sum7_2b_2,[1,0,2,3]) - \
np.transpose(sum7_2b_2,[0,1,3,2]) + \
np.transpose(sum7_2b_2,[1,0,3,2])
sum7_2b_1 = ncon([eta3B, occI, W], [(4,5,-1,6,7,-4),
(4,5,6,7,0,1,2,3),
(2,3,-2,0,1,-3)]).numpy()
sum7_2b = sum7_2b_1 - np.transpose(sum7_2b_1,[1,0,2,3]) - \
np.transpose(sum7_2b_1,[0,1,3,2]) + \
np.transpose(sum7_2b_1,[1,0,3,2])
# sum7_2b_1 = ncon([occI, W], [(-1,-2,-3,-4,0,1,2,3), (2,3,-5,0,1,-6)]).numpy()
# sum7_2b_2 = ncon([eta3B, sum7_2b_1], [(0,1,-1,2,3,-4),(2,3,-2,0,1,-3)]).numpy()
# sum7_2b = sum7_2b_2 - np.transpose(sum7_2b_2,[1,0,2,3]) - \
# np.transpose(sum7_2b_2,[0,1,3,2]) + \
# np.transpose(sum7_2b_2,[1,0,3,2])
dG += sum4_2b + 0.5*sum5_2b + (1/6)*sum6_2b + 0.25*sum7_2b
......
......@@ -106,10 +106,10 @@ class WegnerGenerator(Generator):
eta2B) -- two-body generator"""
partition = self.decouple_OD()
fd = partition[0]
fod = partition[1]
Gd = partition[2]
God = partition[3]
fd = partition[0].astype(np.float32)
fod = partition[1].astype(np.float32)
Gd = partition[2].astype(np.float32)
God = partition[3].astype(np.float32)
holes = self._holes
particles = self._particles
......@@ -314,15 +314,22 @@ class WegnerGenerator3B(WegnerGenerator):
sum4_1b = sum4_1b_3 - sum4_1b_4
# fifth term
sum5_1b_1 = ncon([occF, Wd.astype(np.float32)],
[(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)]).numpy()
sum5_1b_2 = ncon([occF, Wod.astype(np.float32)],
[(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)]).numpy()
sum5_1b_3 = ncon([sum5_1b_1, Wod.astype(np.float32)],
[(0,1,-1,2,3,4), (2,3,4,0,1,-2)]).numpy()
sum5_1b_4 = ncon([sum5_1b_2, Wd.astype(np.float32)],
[(0,1,-1,2,3,4), (2,3,4,0,1,-2)]).numpy()
sum5_1b = sum5_1b_3 - sum5_1b_4
sum5_1b_1 = ncon([Wd, occF, Wod], [(5,6,-1,7,8,9),
(5,6,7,8,9,0,1,2,3,4),
(2,3,4,0,1,-2)]).numpy()
sum5_1b_2 = ncon([Wod, occF, Wd], [(5,6,-1,7,8,9),
(5,6,7,8,9,0,1,2,3,4),
(2,3,4,0,1,-2)]).numpy()
sum5_1b = sum5_1b_1 - sum5_1b_2
# sum5_1b_1 = ncon([occF, Wd.astype(np.float32)],
# [(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)]).numpy()
# sum5_1b_2 = ncon([occF, Wod.astype(np.float32)],
# [(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)]).numpy()
# sum5_1b_3 = ncon([sum5_1b_1, Wod.astype(np.float32)],
# [(0,1,-1,2,3,4), (2,3,4,0,1,-2)]).numpy()
# sum5_1b_4 = ncon([sum5_1b_2, Wd.astype(np.float32)],
# [(0,1,-1,2,3,4), (2,3,4,0,1,-2)]).numpy()
# sum5_1b = sum5_1b_3 - sum5_1b_4
eta1B += 0.25*sum4_1b + (1/12)*sum5_1b
......@@ -335,28 +342,53 @@ class WegnerGenerator3B(WegnerGenerator):
sum4_2b = sum4_2b_3 - sum4_2b_4
#fifth term
sum5_2b_1 = ncon([occG, God], [(-1,-2,-4,0,1,2), (1,2,0,-3)]).numpy()
sum5_2b_2 = ncon([occG, Gd], [(-1,-2,-4,0,1,2), (1,2,0,-3)]).numpy()
sum5_2b_3 = ncon([Wd, sum5_2b_1], [(0,-1,-2,1,2,-4), (1,2,0,-3)]).numpy()
sum5_2b_4 = ncon([Wod, sum5_2b_2], [(0,-1,-2,1,2,-4), (1,2,0,-3)]).numpy()
sum5_2b_5 = sum5_2b_3 - sum5_2b_4
sum5_2b = sum5_2b_5 - np.transpose(sum5_2b_5, [3,2,0,1]) - \
np.transpose(sum5_2b_5, [0,1,3,2]) + \
np.transpose(sum5_2b_5, [2,3,0,1])
sum5_2b_1 = ncon([Wd, occG, God], [(3,-1,-2,4,5,-4),
(3,4,5,0,1,2),
(1,2,0,-3)]).numpy()
sum5_2b_2 = ncon([Wod, occG, Gd], [(3,-1,-2,4,5,-4),
(3,4,5,0,1,2),
(1,2,0,-3)]).numpy()
sum5_2b = sum5_2b_2 - np.transpose(sum5_2b_2, [3,2,0,1]) - \
np.transpose(sum5_2b_2, [0,1,3,2]) + \
np.transpose(sum5_2b_2, [2,3,0,1])
# sum5_2b_1 = ncon([occG, God], [(-1,-2,-4,0,1,2), (1,2,0,-3)]).numpy()
# sum5_2b_2 = ncon([occG, Gd], [(-1,-2,-4,0,1,2), (1,2,0,-3)]).numpy()
# sum5_2b_3 = ncon([Wd, sum5_2b_1], [(0,-1,-2,1,2,-4), (1,2,0,-3)]).numpy()
# sum5_2b_4 = ncon([Wod, sum5_2b_2], [(0,-1,-2,1,2,-4), (1,2,0,-3)]).numpy()
# sum5_2b_5 = sum5_2b_3 - sum5_2b_4
# sum5_2b = sum5_2b_5 - np.transpose(sum5_2b_5, [3,2,0,1]) - \
# np.transpose(sum5_2b_5, [0,1,3,2]) + \
# np.transpose(sum5_2b_5, [2,3,0,1])
#sixth term
sum6_2b_1 = ncon([occH, Wod], [(-1,-2,-3,-4,0,1,2,3),(1,2,3,0,-5,-6)]).numpy()
sum6_2b_2 = ncon([occH, Wod], [(-3,-4,-5,-6,0,1,2,3),(0,-1,-2,1,2,3)]).numpy()
sum6_2b_3 = ncon([Wd, sum6_2b_1], [(0,-1,-2,1,2,3), (1,2,3,0,-3,-4)]).numpy()
sum6_2b_4 = ncon([Wd, sum6_2b_2], [(1,2,3,0,-3,-4), (0,-1,-2,1,2,3)]).numpy()
sum6_2b = sum6_2b_3 - sum6_2b_4
sum6_2b_1 = ncon([Wd, occH, Wod], [(4,-1,-2,5,6,7),
(4,5,6,7,0,1,2,3),
(1,2,3,0,-3,-4)]).numpy()
sum6_2b_2 = ncon([Wd, occH, Wod], [(5,6,7,4,-3,-4),
(4,5,6,7,0,1,2,3),
(0,-1,-2,1,2,3)]).numpy()
sum6_2b = sum6_2b_1 - sum6_2b_2
# sum6_2b_1 = ncon([occH, Wod], [(-1,-2,-3,-4,0,1,2,3),(1,2,3,0,-5,-6)]).numpy()
# sum6_2b_2 = ncon([occH, Wod], [(-3,-4,-5,-6,0,1,2,3),(0,-1,-2,1,2,3)]).numpy()
# sum6_2b_3 = ncon([Wd, sum6_2b_1], [(0,-1,-2,1,2,3), (1,2,3,0,-3,-4)]).numpy()
# sum6_2b_4 = ncon([Wd, sum6_2b_2], [(1,2,3,0,-3,-4), (0,-1,-2,1,2,3)]).numpy()
# sum6_2b = sum6_2b_3 - sum6_2b_4
#seventh term
sum7_2b_1 = ncon([occI, Wod], [(-1,-2,-3,-4,0,1,2,3), (2,3,-5,0,1,-6)]).numpy()
sum7_2b_2 = ncon([Wd, sum7_2b_1], [(0,1,-1,2,3,-4),(2,3,-2,0,1,-3)]).numpy()
sum7_2b = sum7_2b_2 - np.transpose(sum7_2b_2,[1,0,2,3]) - \
np.transpose(sum7_2b_2,[0,1,3,2]) + \
np.transpose(sum7_2b_2,[1,0,3,2])
sum7_2b_1 = ncon([Wd, occI, Wod], [(4,5,-1,6,7,-4),
(4,5,6,7,0,1,2,3),
(2,3,-2,0,1,-3)]).numpy()
sum7_2b = sum7_2b_1 - np.transpose(sum7_2b_1,[1,0,2,3]) - \
np.transpose(sum7_2b_1,[0,1,3,2]) + \
np.transpose(sum7_2b_1,[1,0,3,2])
# sum7_2b_1 = ncon([occI, Wod], [(-1,-2,-3,-4,0,1,2,3), (2,3,-5,0,1,-6)]).numpy()
# sum7_2b_2 = ncon([Wd, sum7_2b_1], [(0,1,-1,2,3,-4),(2,3,-2,0,1,-3)]).numpy()
# sum7_2b = sum7_2b_2 - np.transpose(sum7_2b_2,[1,0,2,3]) - \
# np.transpose(sum7_2b_2,[0,1,3,2]) + \
# np.transpose(sum7_2b_2,[1,0,3,2])
eta2B += sum4_2b + 0.5*sum5_2b + (1/6)*sum6_2b + 0.25*sum7_2b
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment