Commit 092e815c authored by Jacob August Davison's avatar Jacob August Davison
Browse files

you can now specify the generator in main()

parent 1c60dc7c
......@@ -280,7 +280,10 @@ def eta_white(f, Gamma, user_data):
denom = f[a,a] - f[i,i] + Gamma[idx2B[(a,i)], idx2B[(a,i)]]
val = f[a,i]/denom
eta1B[a, i] = val
eta1B[i, a] = -val
eta1B[i, a] = -val
# if denom < 1:
# print("one body {}{}, ".format(a,i), denom)
#print(denom)
# two-body part of the generator
eta2B = np.zeros_like(Gamma)
......@@ -303,6 +306,10 @@ def eta_white(f, Gamma, user_data):
eta2B[idx2B[(a,b)],idx2B[(i,j)]] = val
eta2B[idx2B[(i,j)],idx2B[(a,b)]] = -val
#print(denom)
if denom < 1:
print("two body {}{}{}{}, ".format(a,b,i,j), denom)
#print(Gamma[idx2B[(a,b)], idx2B[(i,j)]], a,b,i,j)
return eta1B, eta2B
......@@ -322,6 +329,8 @@ def eta_white_mp(f, Gamma, user_data):
val = f[a,i]/denom
eta1B[a, i] = val
eta1B[i, a] = -val
if denom < 1:
print("one body {}{}, ".format(a,i), denom)
# two-body part of the generator
eta2B = np.zeros_like(Gamma)
......@@ -339,6 +348,10 @@ def eta_white_mp(f, Gamma, user_data):
eta2B[idx2B[(a,b)],idx2B[(i,j)]] = val
eta2B[idx2B[(i,j)],idx2B[(a,b)]] = -val
if denom < 1:
print("two body {}{}{}{}, ".format(a,b,i,j), denom)
return eta1B, eta2B
def eta_white_atan(f, Gamma, user_data):
......@@ -830,8 +843,44 @@ def calc_mbpt3(f, Gamma, user_data):
#------------------------------------------------------------------------------
# Main program
#------------------------------------------------------------------------------
def dump_flow_data(count, s_val, E, f, Gamma, user_data):
bas1B = user_data["bas1B"]
idx2B = user_data["idx2B"]
holes = user_data['holes']
dim1B = len(bas1B)
Gamma_rank4 = np.zeros((dim1B, dim1B, dim1B, dim1B))
basis = bas1B
for p in basis:
for q in basis:
for r in basis:
for s in basis:
Gamma_rank4[p,q,r,s] = Gamma[idx2B[(p,q)], idx2B[(r,s)]]
# Compute 2B vacuum coeffs
H2B = Gamma_rank4
# Compute 1B vacuum coeffs
H1B = f - np.trace(Gamma_rank4[np.ix_(basis,holes,basis,holes)], axis1=1,axis2=3)
# Compute 0B vacuum
Gamma_trace = 0
for i in holes:
for j in holes:
Gamma_trace += Gamma_rank4[i,j,i,j]
f_trace = 0
for i in holes:
f_trace += f[i,i]
H0B = E - f_trace + 0.5*Gamma_trace
print(s_val)
pickle.dump((s_val, H0B, H1B, H2B), open('vac_coeffs_reference_c{}.p'.format(count), 'wb'))
#@profile
def main(n_holes, g=0.5):
def main(n_holes, g=2):
# grab delta and g from the command line
# delta = float(argv[1])
# g = float(argv[2])
......@@ -847,7 +896,7 @@ def main(n_holes, g=0.5):
# 1st state
holes = np.arange(particles)
particles = np.arange(particles,dim1B)
print(holes, particles)
# 2nd state
# holes = [0,1,4,5]
# particles = [2,3,6,7]
......@@ -893,7 +942,7 @@ def main(n_holes, g=0.5):
"dE": 0.0, # and main routine
"calc_eta": eta_wegner, # specify the generator (function object)
"calc_eta": eta_white, # specify the generator (function object)
"calc_rhs": flow_imsrg2 # specify the right-hand side and truncation
}
......@@ -937,10 +986,12 @@ def main(n_holes, g=0.5):
norm_fod = calc_fod_norm(f, user_data)
norm_Gammaod = calc_Gammaod_norm(Gamma, user_data)
s_val = solver.t
if count % 10 == 0:
print("%8.5f %14.8f %14.8f %14.8f %14.8f %14.8f %14.8f %14.8f %14.8f"%(
solver.t, E , DE2, DE3, E+DE2+DE3, user_data["dE"], user_data["eta_norm"], norm_fod, norm_Gammaod))
print(s_val)
dump_flow_data(count, s_val, E, f, Gamma, user_data)
count += 1
if abs(DE2/E) < 10e-8: break
......
......@@ -45,6 +45,7 @@ def get_vacuum_coeffs(E, f, G, basis, holes):
Gnode[0] ^ Gnode[2]
Gnode[1] ^ Gnode[3]
result = Gnode @ Gnode
H0B = E - np.trace(f[np.ix_(holes,holes)]) + 0.5*result.tensor
return (H0B, H1B, H2B)
......@@ -141,7 +142,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):
def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1, flow_data_log=0, generator='wegner'):
"""Main method uses scipy.integrate.ode to solve the IMSRG(2) flow
equations."""
......@@ -156,7 +157,12 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1):
ha = PairingHamiltonian2B(n_holes, n_particles, ref=ref, d=d, g=g, pb=pb)
ot = OccupationTensors(ha.sp_basis, ha.reference)
wg = WegnerGenerator(ha, ot)
generator_dict = {'wegner':WegnerGenerator(ha, ot),
'white':WhiteGenerator(ha),
'white_mp':WhiteGeneratorMP(ha)}
wg = generator_dict[generator] #WegnerGenerator(ha, ot)
fl = Flow_IMSRG2(ha, ot)
initf = time.time() # finish instantiation timer
......@@ -166,13 +172,14 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1):
if verbose:
print("""Pairing model IM-SRG(2) flow:
Generator = {}
d = {:2.4f}
g = {:2.4f}
pb = {:2.4f}
SP basis size = {:2d}
n_holes = {:2d}
n_particles = {:2d}
ref = {d}""".format(ha.d, ha.g, ha.pb, ha.n_sp_states,
ref = {d}""".format(generator, ha.d, ha.g, ha.pb, ha.n_sp_states,
len(ha.holes), len(ha.particles),
d=ref) )
if verbose:
......@@ -185,12 +192,12 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1):
pickle.dump( coeffs, open( "./vac_coeffs_unevolved.p", "wb" ) )
solver = ode(derivative,jac=None)
solver.set_integrator('vode', method='bdf', order=5, nsteps=500)
solver.set_integrator('vode', method='bdf', order=5, nsteps=1000)
solver.set_f_params(ha, ot, wg, fl)
solver.set_initial_value(y0, 0.)
sfinal = 100
ds = 0.1
sfinal = 50
ds = 1
s_vals = []
E_vals = []
......@@ -198,15 +205,22 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1):
convergence = 0
while solver.successful() and solver.t < sfinal:
#print('solver success,', solver.successful())
ys = solver.integrate(sfinal, step=True)
Es, fs, Gs = ravel(ys, ha.n_sp_states)
s_vals.append(solver.t)
E_vals.append(Es)
iters += 1
# if iters == 176:
# break
if iters %10 == 0 and verbose: print("iter: {:>6d} \t scale param: {:0.4f} \t E = {:0.8f}".format(iters, solver.t, Es))
if flow_data_log:# and iters %10 == 0:
H0B, H1B, H2B = get_vacuum_coeffs(Es, fs, Gs, ha.sp_basis, ha.holes)
fname = 'vac_coeffs_flow_c{}.p'.format(iters)
pickle.dump((solver.t, H0B, H1B, H2B), open(fname, 'wb'))
# if iters %20 == 0 and verbose:
# coeffs = get_vacuum_coeffs(Es, fs, Gs, ha.sp_basis, ha.holes)
......@@ -218,9 +232,9 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1):
convergence = 1
break
if len(E_vals) > 100 and abs(E_vals[-1] - E_vals[-2]) > 1:
if verbose: 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:
# if verbose: print("---- Energy diverged at iter {:>06d} with energy {:3.8f}\n".format(iters,E_vals[-1]))
# break
if iters > 20000:
if verbose: print("---- Energy diverged at iter {:>06d} with energy {:3.8f}\n".format(iters,E_vals[-1]))
......@@ -229,7 +243,8 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1):
if iters % 1000 == 0 and verbose:
print('Iteration {:>06d}'.format(iters))
iters += 1
#print(solver.successful())
flowf = time.time()
end = time.time()
time_str = "{:2.5f}".format(end-start)
......@@ -294,13 +309,14 @@ 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]]
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.985*np.asarray(refs[0]) + (1.0-0.985)/4.0*(np.asarray(refs[1]) + np.asarray(refs[2]) + np.asarray(refs[3]) + np.asarray(refs[4]))
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])
test_exact('plots_exact_2b/', main)
main(4,4,g=2, ref=ref, flow_data_log=0, generator='white')
# H1B_true, H2B_true = pickle.load(open('comparison.p','rb'))
# H1B, H2B = pickle.load(open('vac_coeffs_unevolved.p', 'rb'))
......
......@@ -71,16 +71,16 @@ class Flow_IMSRG2(Flow):
df, -- one-body tensor
dG) -- two-body tensor"""
assert isinstance(gen, Generator), "Arg 0 must be WegnerGenerator object"
assert isinstance(gen, Generator), "Arg 0 must be Generator object"
# f = self.f
# G = self.G
f = gen.f.astype(np.float32)
G = gen.G.astype(np.float32)
partition = gen.calc_eta()
eta1B = partition[0]#.astype(np.float32)
eta2B = partition[1]#.astype(np.float32)
eta1B, eta2B = gen.calc_eta()
# eta1B = partition[0]#.astype(np.float32)
# eta2B = partition[1]#.astype(np.float32)
# occA = occ_t.occA
# occB = occ_t.occB
......@@ -208,7 +208,7 @@ class Flow_IMSRG2(Flow):
# sum3_2b_5 = sum3_2b_1 - sum3_2b_2 - sum3_2b_3 + sum3_2b_4
# sum3_2b = tn.ncon([occA, sum3_2b_5], [(0, 1, -1, -2), (0, 1, -3, -4)])#.numpy()
dG = sum1_2b + 0.5*sum2_2b + sum3_2b
dG = sum1_2b + 0.5*sum2_2b - sum3_2b
return (dE, df, dG)
......
......@@ -184,11 +184,11 @@ def test_exact(plots_dir, main):
# plt.title('Correlation energy with pb = {:2.4f}'.format(pb))
# plt.savefig(plots_dir+'pb{:2.4f}.png'.format(pb))
if not os.path.exists('data_pickles/'):
os.mkdir('data_pickles/')
if not os.path.exists('data_pickles1/'):
os.mkdir('data_pickles1/')
fulldata = [g_vals, E_exacts, E_corrs, E_pys]
with open('data_pickles/fulldata.pickle', 'wb') as f:
with open('data_pickles1/fulldata.pickle', 'wb') as f:
pickle.dump(fulldata, f, pickle.HIGHEST_PROTOCOL)
#plt.close()
......
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