Skip to content
GitLab
Commit c65ad832 authored by Jacob August Davison's avatar Jacob August Davison
Browse files

many changes; optimized for IM-SRG(2)

parent 99abe768
Hide whitespace changes
Inline Side-by-side
TODO.md
View file @ c65ad832
......@@ -20,10 +20,15 @@
* [ ] scan reference state configurations to find ground state for values of pb and g
* [x] using TN to implement the IMSRG(3)
* [ ] data for committee meeting
* [ ] get scaling data for various IMSRG(2) implementations
* [ ] benchmark IMSRG(2) results against full CI and Heiko's python code (check energy convergence discrepancy)
* [X] data for committee meeting
* [X] get scaling data for various IMSRG(2) implementations
* [X] benchmark IMSRG(2) results against full CI and Heiko's python code (check energy convergence discrepancy)
* [ ] test IMSRG(3)
* [ ] TT decompose IMSRG(3)-relevant occupation tensors (should speed things up)
* [ ] compare energy divergence to IMSRG(2)
# [ ] implement batch processing
* [ ] modulate the flow equations (need more control over memory allocation)
* [ ] might move away from NCON framework in favor of TensorNetwork syntax
* [ ] need function to isolate batches, run data through flow equations, and bring them back together
benchmarking_wd/imsrg_pairing.py
View file @ c65ad832
......@@ -22,6 +22,7 @@ from scipy.integrate import odeint, ode
from sys import argv
import time
#-----------------------------------------------------------------------------------
# basis and index functions
#-----------------------------------------------------------------------------------
......@@ -899,7 +900,7 @@ def main(n_holes, g=0.5):
H1B, H2B = pairing_hamiltonian(delta, g, user_data)
E, f, Gamma = normal_order(H1B, H2B, user_data)
print(E)
#print(E)
# reshape Hamiltonian into a linear array (initial ODE vector)
y0 = np.append([E], np.append(reshape(f, -1), reshape(Gamma, -1)))
......@@ -907,7 +908,7 @@ def main(n_holes, g=0.5):
t = 1
dy = derivative_wrapper(t, y0, user_data)
dE, df, dG = get_operator_from_y(dy, dim1B, dim1B*dim1B)
print(dE)
#print(dE)
# integrate flow equations
solver = ode(derivative_wrapper,jac=None)
......@@ -952,7 +953,11 @@ def main(n_holes, g=0.5):
# make executable
#------------------------------------------------------------------------------
if __name__ == "__main__":
main(4)
for n in range(2,14,2):
ti = time.time()
main(n)
tf = time.time()
print('Total time: {:2.5f}'.format(tf-ti))
import pytest
......
main.py
View file @ c65ad832
......@@ -125,7 +125,7 @@ def ravel(y, bas_len):
return(ravel_E, ravel_f, ravel_G)
# @profile
def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0, verbose=1):
"""Main method uses scipy.integrate.ode to solve the IMSRG(2) flow
equations."""
......@@ -145,20 +145,22 @@ def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
initf = time.time() # finish instantiation timer
print("Initialized objects in {:2.4f} seconds\n".format(initf-initi))
print("""Pairing model IM-SRG(2) flow:
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,
len(ha.holes), len(ha.particles),
d=ref) )
print("Flowing...")
if verbose:
print("Initialized objects in {:2.4f} seconds\n".format(initf-initi))
if verbose:
print("""Pairing model IM-SRG(2) flow:
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,
len(ha.holes), len(ha.particles),
d=ref) )
if verbose:
print("Flowing...")
flowi = time.time()
# --- Solve the IM-SRG flow
......@@ -186,81 +188,52 @@ def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
iters += 1
# 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 iters %10 == 0 and verbose: 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]))
if verbose: 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]))
if verbose: 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]))
if verbose: print("---- Energy diverged at iter {:>06d} with energy {:3.8f}\n".format(iters,E_vals[-1]))
break
if iters % 1000 == 0:
if iters % 1000 == 0 and verbose:
print('Iteration {:>06d}'.format(iters))
flowf = time.time()
end = time.time()
time_str = "{:2.5f}\n".format(end-start)
print("IM-SRG(2) converged in {:2.5f} seconds".format(flowf-flowi))
time_str = "{:2.5f}".format(end-start)
if verbose: print("IM-SRG(2) converged in {:2.5f} seconds".format(flowf-flowi))
del ha, ot, wg, fl, solver, y0, sfinal, ds
return (convergence, iters, d, g, pb, n_holes+n_particles, s_vals, E_vals, time_str)
if __name__ == '__main__':
#test_refs('logs_refs\\')
# print(ci_matrix.exact_diagonalization(1.0, 0.5,0.0))
#test = main(4,4)
# occt = OccupationTensors(np.array([0,1,2,3,4,5,6,7]), np.array([1,1,1,1,0,0,0,0]))
# bas1B = [0,1,2,3,4,5,6,7]
# occE = occt.occE
# for a in bas1B:
# for b in bas1B:
# for c in bas1B:
# for d in bas1B:
# for e in bas1B:
# for f in bas1B:
# val = occE[a,b,c,d,e,f]
# if val != 0:
# with open('occE_nonzero.txt', 'a') as fi:
# fi.write('%s %s %s %s %s %s -- %s\n' % (a,b,c,d,e,f,val))
test_exact('plots_exact_2b/', main)
#main(4,4)
# bas1B = np.array([0,1,2,3,4,5,6,7])
# ref = np.array([1,1,1,1,0,0,0,0],dtype=np.float16)
# occt = OccupationTensors(bas1B,ref)
# occF = occt.occF
# test1 = np.einsum('i,j,k,l,m->ijklm',ref,ref,(1-ref),(1-ref),(1-ref))
# test1e = np.einsum('ijklm,nopqr->ijklmnopqr',test1,test1)
# test2 = np.einsum('i,j,k,l,m->ijklm',(1-ref),(1-ref),ref,ref,ref)
# test2e = np.einsum('ijklm,nopqr->ijklmnopqr',test2,test2)
# test = test1e + test2e
# # test = np.einsum('ijklm,nopqr->ijklmnopqr',test3,test3)
# print(occF.dtype, test.dtype)
# print(occF.shape, test.shape)
# print(np.array_equal(occF, test))
#test_exact('plots_exact_2b/', main)
main(4,4)
#holes = int(sys.argv[1])
# print('convergence,iters,d,g,pb,num states,GS energy,total time')
# for num_pairs in range(1,21):
# #print(n)55
# #n = int(num_sp/2)
# #holes = num_pairs*2
# particles = num_pairs*2
# convergence,iters,d,g,pb,sp_states,s_vals,E_vals,time_str = main(holes, particles, verbose=0)
# print('{},{},{},{},{},{},{},{}'.format(convergence,iters,d,g,pb,sp_states,E_vals[-1],time_str))
# for a in bas1B:
# for b in bas1B:
# for c in bas1B:
# for d in bas1B:
# for e in bas1B:
# for f in bas1B:
# val1 = test[a,b,c,d,e,f]
# val2 = occC[a,b,c,d,e,f]
# if val1 != val2:
# print('conflict found for ',a,b,c,d,e,f'-',val1,val2)
# with open('occE_nonzero_test.txt', 'a') as fi:
# fi.write('%s %s %s %s %s -- %s\n' % (a,b,c,d,e,val))
# del convergence, iters, d, g, pb, sp_states, s_vals, E_vals, time_str
oop_imsrg/generator.py
View file @ c65ad832
......@@ -42,6 +42,8 @@ class WegnerGenerator(Generator):
self._occC = occ_t.occC
self._occD = occ_t.occD
# self._occRef1 =
@property
def f(self):
"""Returns:
......@@ -85,6 +87,7 @@ class WegnerGenerator(Generator):
particles = self._particles
# - Decouple off-diagonal 1B and 2B pieces
# fod1 = tn.ncon([])
fod = np.zeros(f.shape, dtype=np.float32)
fod[np.ix_(particles, holes)] += f[np.ix_(particles, holes)]
fod[np.ix_(holes, particles)] += f[np.ix_(holes, particles)]
......
oop_imsrg/hamiltonian.py
View file @ c65ad832
......@@ -29,6 +29,7 @@ class PairingHamiltonian2B(Hamiltonian):
Keyword arguments:
ref -- the reference state. must match dimensions imposed by arugments (default: [1,1,1,1,0,0,0,0])
p -
d -- the energy level spacing (default: 1.0)
g -- the pairing strength (default: 0.5)
pb -- strength of the pair-breaking term (operates in double particle basis) (default: 0.0)"""
......
oop_imsrg/occupation_tensors.py
View file @ c65ad832
......@@ -29,16 +29,22 @@ class OccupationTensors(object):
self._occB4 = self.__get_occB(flag=1)
self._occC = self.__get_occC()
self._occD = self.__get_occD(flag=1)
self._occE = self.__get_occE()
self._occF = self.__get_occF()
self._occG = self.__get_occG()
self._occH = self.__get_occH()
self._occI = self.__get_occI()
self._occJ = self.__get_occJ()
# self._occE = self.__get_occE()
# self._occF = self.__get_occF()
# self._occG = self.__get_occG()
# self._occH = self.__get_occH()
# self._occI = self.__get_occI()
# self._occJ = self.__get_occJ()
if not os.path.exists("occ_storage/"):
os.mkdir("occ_storage/")
# @property
# def occRef1(self):
# """Returns:
# occRef1 -- represents n_a(1-n_b)."""
@property
......@@ -130,7 +136,19 @@ class OccupationTensors(object):
# ---- BUILD OCCUPATION TENSORS ---
# def __get_occRef1(self):
# """Builds the occupation tensor occRef1, necessary for occupation number
# representation of the Hamiltonian.
# Returns:
# occRef1 -- n_a(1-n-b)
# """
#@jit#(nopython=True)
def __get_occA(self, flag=0):
"""Builds the occupation tensor occA.
......@@ -271,7 +289,7 @@ class OccupationTensors(object):
Ga = tn.Node(np.append(1-ref[:,np.newaxis], np.ones((n,1)),axis=1).astype(int))
Gb = tn.Node(np.transpose(np.append(np.ones((n,1)), -1*ref[:,np.newaxis],axis=1).astype(int)))
Gab = tn.ncon([Ga,Gb], [(-1,1),(1,-2)])
final = tn.outer_product(Gab, tn.Node(np.ones((8,8))))
final = tn.outer_product(Gab, tn.Node(np.ones((n,n))))
occB = final
......@@ -396,10 +414,15 @@ class OccupationTensors(object):
# (1-ref[c])*(1-ref[d])
Ga = tn.Node(np.array([ [1,0],[1,0],[1,0],[1,0],[0,0],[0,0],[0,0],[0,0] ]))
Gb = tn.Node(np.transpose(np.array([ [1,0],[1,0],[1,0],[1,0],[0,0],[0,0],[0,0],[0,0] ])))
Gc = tn.Node(np.array([ [0,0],[0,0],[0,0],[0,0],[1,0],[1,0],[1,0],[1,0] ]))
Gd = tn.Node(np.transpose(np.array([ [0,0],[0,0],[0,0],[0,0],[1,0],[1,0],[1,0],[1,0] ])))
# Ga = tn.Node(np.array([ [1,0],[1,0],[1,0],[1,0],[0,0],[0,0],[0,0],[0,0] ]))
# Gb = tn.Node(np.transpose(np.array([ [1,0],[1,0],[1,0],[1,0],[0,0],[0,0],[0,0],[0,0] ])))
# Gc = tn.Node(np.array([ [0,0],[0,0],[0,0],[0,0],[1,0],[1,0],[1,0],[1,0] ]))
# Gd = tn.Node(np.transpose(np.array([ [0,0],[0,0],[0,0],[0,0],[1,0],[1,0],[1,0],[1,0] ])))
Ga = tn.Node(np.transpose(np.append(ref[np.newaxis,:], np.zeros((1,n)), axis=0).astype(int)))
Gb = tn.Node(np.transpose(Ga.tensor))
Gc = tn.Node(np.transpose(np.append(ref[::-1][np.newaxis,:],np.zeros((1,n)), axis=0).astype(int)))
Gd = tn.Node(np.transpose(Gc.tensor))
Gabcd = tn.ncon([Ga,Gb,Gc,Gd], [(-1,1),(1,-2),(-3,2),(2,-4)])
......@@ -435,10 +458,15 @@ class OccupationTensors(object):
# occD[a,b,c,d] = ref[a]*ref[b]*\
# (1-ref[c])*(1-ref[d])
Ga = tn.Node(np.array([ [1,0],[1,0],[1,0],[1,0],[0,0],[0,0],[0,0],[0,0] ]))
Gb = tn.Node(np.transpose(np.array([ [1,0],[1,0],[1,0],[1,0],[0,0],[0,0],[0,0],[0,0] ])))
Gc = tn.Node(np.array([ [0,0],[0,0],[0,0],[0,0],[1,0],[1,0],[1,0],[1,0] ]))
Gd = tn.Node(np.transpose(np.array([ [0,0],[0,0],[0,0],[0,0],[1,0],[1,0],[1,0],[1,0] ])))
# Ga = tn.Node(np.array([ [1,0],[1,0],[1,0],[1,0],[0,0],[0,0],[0,0],[0,0] ]))
# Gb = tn.Node(np.transpose(np.array([ [1,0],[1,0],[1,0],[1,0],[0,0],[0,0],[0,0],[0,0] ])))
# Gc = tn.Node(np.array([ [0,0],[0,0],[0,0],[0,0],[1,0],[1,0],[1,0],[1,0] ]))
# Gd = tn.Node(np.transpose(np.array([ [0,0],[0,0],[0,0],[0,0],[1,0],[1,0],[1,0],[1,0] ])))
Ga = tn.Node(np.transpose(np.append(ref[np.newaxis,:], np.zeros((1,n)), axis=0).astype(int)))
Gb = tn.Node(np.transpose(Ga.tensor))
Gc = tn.Node(np.transpose(np.append(ref[::-1][np.newaxis,:],np.zeros((1,n)), axis=0).astype(int)))
Gd = tn.Node(np.transpose(Gc.tensor))
Gabcd = tn.ncon([Ga,Gb,Gc,Gd], [(-1,1),(1,-2),(-3,2),(2,-4)])
......
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