Skip to content

GitLab

Commit 4074bb3b authored by Davison, Jacob's avatar Davison, Jacob
Browse files

looking for flowing coincidence

parent 61d5449b
Hide whitespace changes
Inline Side-by-side
main.py
View file @ 4074bb3b
......@@ -55,6 +55,8 @@ def get_vacuum_coeffs(E, f, G, basis, holes):
return (H0B, H1B, H2B)
# @profile
def derivative(t, y, inputs):
"""Defines the derivative to pass into ode object.
......@@ -142,7 +144,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.
......@@ -182,7 +184,7 @@ 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)
......@@ -253,7 +255,28 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1, flow_dat
dens_weights,v = np.linalg.eigh(fci_hme)
if verbose:
print("iter, \t s, \t E, \t SS0B, \t SS1B, \t SS2B, \t ||eta1B||, \t ||eta2B||, \t commute1bd, \t commute1bod, \t commute2bd,\t commute2bod")
print_columns = ['iter',
's',
'E',
'CI_gs',
'0bN_SS',
'<SS>',
'0b_SS',
'1bc_SS',
'2bc_SS',
'||eta1b||',
'||eta2b||',
'commute1bd',
'commute1bod',
'commute2bd',
'commute2bod']
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:
......@@ -308,22 +331,44 @@ def main(n_holes, n_particles, ref=[], d=1.0, g=0.5, pb=0.0, verbose=1, flow_dat
# 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))
print("{:>6d}, \t {: .8f}, \t {: .8f}, \t {: .8f}, \t {: .8f}, \t {: .8f}, \t {: .8f}, \t {: .8f}, \t {: .8f}, \t {: .8f}, \t {: .8f}, \t {: .8f}".format(iters,
solver.t,
Es,
SS0B,
contract_1b,
contract_2b,
norm_eta1B,
norm_eta2B,
0.0,
0.0,
0.0,
0.0))
# commute1bd,
# commute1bod,
# commute2bd,
# commute2bod))
data_columns = [iters,
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]
# print(column_string.format(*data_columns))
print("{:>6d}, {: .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)
......@@ -432,22 +477,29 @@ if __name__ == '__main__':
# plt.legend(['E(s)/E(s=0)', 'SS(s)/SS(s=0)'])
# plt.savefig('flow_conservation.png')
g = 2.0
pb = 0.01
ensemble = re.ReferenceEnsemble(8, g, pb, 'white', '', 'vac_coeffs')
opt_x = ensemble.optimize_reference(4)
ref1 = ensemble.refs.T.dot(opt_x)
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[:,0]
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])
#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)
main(4,4, g=g, ref=ref1, pb=pb, generator='white')
#ref = 0.8*basis[0,:] + 0.2*basis[1,:]
#ref = basis.T.dot(v0*v0)
ref = basis[0,:]
#ref = pickle.load(open('reference_g2.00_pb0.01_4-4.p', 'rb'))
main(4,4, g=g, ref=ref, pb=pb, generator='white')
print('FCI ENERGY = {: .8f}'.format(w[0]))
data = pickle.load(open('expect_flow.p', 'rb'))
fig = plt.figure(figsize=(8,4))
......
oop_imsrg/flow.py
View file @ 4074bb3b
......@@ -233,7 +233,7 @@ 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._occC = occ_t.occC
self._occD = occ_t.occD
......
oop_imsrg/generator.py
View file @ 4074bb3b
......@@ -287,9 +287,11 @@ class WegnerGenerator3B(WegnerGenerator):
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
......@@ -298,6 +300,18 @@ class WegnerGenerator3B(WegnerGenerator):
self._occI = occ_t.occI
self._occJ = occ_t.occJ
ref = h.reference
n = len(self._holes)+len(self._particles)
Ga = tn.Node(np.transpose(np.append(ref[np.newaxis,:], np.zeros((1,n)), axis=0).astype(float)))
Gb = tn.Node(np.append(ref[::-1][np.newaxis,:],np.zeros((1,n)), axis=0).astype(float))
self._occRef1 = tn.ncon([Ga,Gb], [(-1,1),(1,-2)]) # n_a(1-n_b)
self._occRef2 = tn.ncon([tn.Node(np.transpose(Gb.tensor)), tn.Node(np.transpose(Ga.tensor))], [(-1,1),(1,-2)]) # (1-n_a)n_b
self._eta1B = np.zeros_like(self.f)
self._eta2B = np.zeros_like(self.G)
@property
def W(self):
"""Returns:
......@@ -371,10 +385,12 @@ class WegnerGenerator3B(WegnerGenerator):
holes = self._holes
particles = self._particles
occA4 = self._occA4
occA = self._occA
occA2 = self._occA2
occB = self._occB
#occA2 = self._occA2
occB4 = self._occB4
occC = self._occC
occC6 = self._occC6
occD = self._occD
occF = self._occF
occG = self._occG
......@@ -385,8 +401,10 @@ class WegnerGenerator3B(WegnerGenerator):
# Calculate 1B generator
# fourth term
sum4_1b_1 = np.matmul(np.transpose(occD,[2,3,0,1]), God)
sum4_1b_2 = np.matmul(np.transpose(occD,[2,3,0,1]), Gd)
occD_transpose = np.transpose(occD.tensor, [2,3,0,1])
sum4_1b_1 = np.multiply(occD_transpose, God)
sum4_1b_2 = np.multiply(occD_transpose, Gd)
sum4_1b_3 = tn.ncon([Wd, sum4_1b_1], [(1,2,-1,3,4,-2),(3,4,1,2)])#.numpy()
sum4_1b_4 = tn.ncon([Wod, sum4_1b_2], [(1,2,-1,3,4,-2),(3,4,1,2)])#.numpy()
sum4_1b = sum4_1b_3 - sum4_1b_4
......@@ -414,8 +432,8 @@ class WegnerGenerator3B(WegnerGenerator):
# Calculate 2B generator
# fourth term
sum4_2b_1 = np.matmul(-1.0*np.transpose(occA2), fod)
sum4_2b_2 = np.matmul(-1.0*np.transpose(occA2), fd)
sum4_2b_1 = np.matmul(-1.0*np.transpose(occA.tensor), fod)
sum4_2b_2 = np.matmul(-1.0*np.transpose(occA.tensor), fd)
sum4_2b_3 = tn.ncon([Wd, sum4_2b_1], [(1,-1,-2,2,-3,-4), (2,1)])#.numpy()
sum4_2b_4 = tn.ncon([Wod, sum4_2b_2], [(1,-1,-2,2,-3,-4), (2,1)])#.numpy()
sum4_2b = sum4_2b_3 - sum4_2b_4
......@@ -500,22 +518,22 @@ class WegnerGenerator3B(WegnerGenerator):
sum1_3b = sum1_3b_4 + sum1_3b_8 + sum1_3b_13
#fourth term
sum4_3b_1 = tn.ncon([Gd, occB, Wod], [(-1,-2,3,4),(3,4,1,2),(1,2,-3,-4,-5,-6)])#.numpy()
sum4_3b_2 = tn.ncon([God, occB, Wd], [(-1,-2,3,4),(3,4,1,2),(1,2,-3,-4,-5,-6)])#.numpy()
sum4_3b_1 = tn.ncon([Gd, occB4, Wod], [(-1,-2,3,4),(3,4,1,2),(1,2,-3,-4,-5,-6)])#.numpy()
sum4_3b_2 = tn.ncon([God, occB4, Wd], [(-1,-2,3,4),(3,4,1,2),(1,2,-3,-4,-5,-6)])#.numpy()
sum4_3b_3 = sum4_3b_1 - sum4_3b_2
sum4_3b = sum4_3b_3 - np.transpose(sum4_3b_3, [1,0,2,3,4,5]) - \
np.transpose(sum4_3b_3, [2,1,0,3,4,5])
#fifth term
sum5_3b_1 = tn.ncon([Gd, occB, Wod], [(3,4,-4,-5),(3,4,1,2),(-1,-2,-3,1,2,-6)])#.numpy()
sum5_3b_2 = tn.ncon([God, occB, Wd], [(3,4,-4,-5),(3,4,1,2),(-1,-2,-3,1,2,-6)])#.numpy()
sum5_3b_1 = tn.ncon([Gd, occB4, Wod], [(3,4,-4,-5),(3,4,1,2),(-1,-2,-3,1,2,-6)])#.numpy()
sum5_3b_2 = tn.ncon([God, occB4, Wd], [(3,4,-4,-5),(3,4,1,2),(-1,-2,-3,1,2,-6)])#.numpy()
sum5_3b_3 = sum5_3b_1 - sum5_3b_2
sum5_3b = sum5_3b_3 - np.transpose(sum5_3b_3, [0,1,2,5,4,3]) - \
np.transpose(sum5_3b_3, [0,1,2,3,5,4])
#sixth term
sum6_3b_1 = tn.ncon([Gd, occA, Wod], [(4,-1,3,-4),(3,4,1,2),(1,-2,-3,2,-5,-6)])#.numpy()
sum6_3b_2 = tn.ncon([God, occA, Wd], [(4,-1,3,-4),(3,4,1,2),(1,-2,-3,2,-5,-6)])#.numpy()
sum6_3b_1 = tn.ncon([Gd, occA4, Wod], [(4,-1,3,-4),(3,4,1,2),(1,-2,-3,2,-5,-6)])#.numpy()
sum6_3b_2 = tn.ncon([God, occA4, Wd], [(4,-1,3,-4),(3,4,1,2),(1,-2,-3,2,-5,-6)])#.numpy()
sum6_3b_3 = sum6_3b_1 - sum6_3b_2
sum6_3b_4 = sum6_3b_3 - np.transpose(sum6_3b_3, [0,1,2,4,3,5]) - \
np.transpose(sum6_3b_3, [0,1,2,5,4,3])
......@@ -528,8 +546,9 @@ class WegnerGenerator3B(WegnerGenerator):
sum7_3b = sum7_3b_1 - sum7_3b_2
#eighth term
sum8_3b_1 = tn.ncon([Wd, occC, Wod], [(4,5,-3,6,-5,-6), (4,5,6,1,2,3), (3,-1,-2,1,2,-4)])#.numpy()
sum8_3b_2 = tn.ncon([Wd, occC, Wod], [(6,-2,-3,4,5,-6), (4,5,6,1,2,3), (-1,1,2,-4,-5,3)])#.numpy()
sum8_3b_1 = tn.ncon([Wd, occC6, Wod], [(4,5,-3,6,-5,-6), (4,5,6,1,2,3), (3,-1,-2,1,2,-4)])#.numpy()
sum8_3b_2 = tn.ncon([Wd, occC6, Wod], [(6,-2,-3,4,5,-6), (4,5,6,1,2,3), (-1,1,2,-4,-5,3)])#.numpy()
sum8_3b_3 = sum8_3b_1 - sum8_3b_2
sum8_3b_4 = sum8_3b_3 - np.transpose(sum8_3b_3, [0,1,2,4,3,5]) - \
np.transpose(sum8_3b_3, [0,1,2,5,4,3])
......
oop_imsrg/occupation_tensors.py
View file @ 4074bb3b
......@@ -30,13 +30,14 @@ class OccupationTensors(object):
self._occB = self.__get_occB()
self._occB4 = self.__get_occB(flag=1)
self._occC = self.__get_occC()
self._occC6 = self.__get_occC(flag=1)
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/")
......@@ -86,6 +87,14 @@ class OccupationTensors(object):
occC -- represents n_a*n_b*(1-n_c) + (1-n_a)*(1-n_b)*n_c"""
return self._occC
@property
def occC6(self):
"""Returns:
occC -- represents n_a*n_b*(1-n_c) + (1-n_a)*(1-n_b)*n_c rank 6"""
return self._occC6
@property
def occD(self):
"""Returns:
......@@ -312,74 +321,74 @@ class OccupationTensors(object):
ref = self._reference
n = len(bas1B)
if flag == 0: # default
# occC = np.zeros((n,n,n,n,n,n),dtype=np.float32)
# for a in bas1B:
# for b in bas1B:
# for c in bas1B:
# occC[a,b,c,a,b,c] = ref[a]*ref[b]*(1-ref[c]) + \
# (1-ref[a])*(1-ref[b])*ref[c]
# print(sys.getsizeof(occC)/10**6)
# TENSOR TRAIN DECOMPOSITION of rank 6 tensor with elements
# given by n_a*n_b + (1 - n_b - n_a)*n_c.
#Ga = tn.Node(np.array([ [1,1,1,1], [1,1,1,1], [1,1,1,1], [1,1,1,1], [0,1,1,0], [0,1,1,0], [0,1,1,0], [0,1,1,0] ]))
# Gb = tn.Node(np.array([ [[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]],
# [[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]],
# [[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]],
# [[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]],
# [[0,0,0,0], [0,1,0,0],[0,0,0,0],[0,0,0,-1]],
# [[0,0,0,0], [0,1,0,0],[0,0,0,0],[0,0,0,-1]],
# [[0,0,0,0], [0,1,0,0],[0,0,0,0],[0,0,0,-1]],
# [[0,0,0,0], [0,1,0,0],[0,0,0,0],[0,0,0,-1]] ]))
#Gc = tn.Node(np.transpose(np.array([ [1,1,1,1], [1,1,1,1], [1,1,1,1], [1,1,1,1], [1,0,0,0], [1,0,0,0], [1,0,0,0], [1,0,0,0] ])))
Ga1 = np.append(ref[:,np.newaxis], np.ones((n,2)),axis=1).astype(self.DATA_TYPE)
Ga2 = np.append(Ga1, ref[:,np.newaxis],axis=1).astype(self.DATA_TYPE)
Ga = tn.Node(Ga2)
Gb1= np.append(ref[np.newaxis,np.newaxis,:],np.zeros((1,1,n)),axis=1).astype(self.DATA_TYPE)
Gb2= np.append(Gb1, np.append(np.zeros((1,1,n)), np.ones((1,1,n)),axis=1), axis=0).astype(self.DATA_TYPE)
Gb3 = np.array([[1,0],[0,-1]])
Gb = tn.Node(np.kron(Gb3, np.transpose(Gb2)))
Gc1 = np.append(np.ones((n,1)), np.repeat(ref[:,np.newaxis],3,axis=1), axis=1).astype(self.DATA_TYPE)
Gc = tn.Node(np.transpose(Gc1))
Gabc = tn.ncon([Ga, Gb, Gc], [(-1, 1), (-2, 1, 2), (2, -3)])
# occC = np.zeros((n,n,n,n,n,n),dtype=np.float32)
# for a in bas1B:
# for b in bas1B:
# for c in bas1B:
# occC[a,b,c,a,b,c] = ref[a]*ref[b]*(1-ref[c]) + \
# (1-ref[a])*(1-ref[b])*ref[c]
# final = tn.outer_product(Gabc, tn.Node(np.ones((n,n,n)))).tensor
# print(sys.getsizeof(occC)/10**6)
# @numba.jit(nopython=True)#, parallel=True)
# def enforce_delta(n, tensor):
# TENSOR TRAIN DECOMPOSITION of rank 6 tensor with elements
# given by n_a*n_b + (1 - n_b - n_a)*n_c.
#Ga = tn.Node(np.array([ [1,1,1,1], [1,1,1,1], [1,1,1,1], [1,1,1,1], [0,1,1,0], [0,1,1,0], [0,1,1,0], [0,1,1,0] ]))
# Gb = tn.Node(np.array([ [[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]],
# [[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]],
# [[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]],
# [[1,0,0,0], [0,1,0,0],[0,0,-1,0],[0,0,0,-1]],
# [[0,0,0,0], [0,1,0,0],[0,0,0,0],[0,0,0,-1]],
# [[0,0,0,0], [0,1,0,0],[0,0,0,0],[0,0,0,-1]],
# [[0,0,0,0], [0,1,0,0],[0,0,0,0],[0,0,0,-1]],
# [[0,0,0,0], [0,1,0,0],[0,0,0,0],[0,0,0,-1]] ]))
#Gc = tn.Node(np.transpose(np.array([ [1,1,1,1], [1,1,1,1], [1,1,1,1], [1,1,1,1], [1,0,0,0], [1,0,0,0], [1,0,0,0], [1,0,0,0] ])))
Ga1 = np.append(ref[:,np.newaxis], np.ones((n,2)),axis=1).astype(self.DATA_TYPE)
Ga2 = np.append(Ga1, ref[:,np.newaxis],axis=1).astype(self.DATA_TYPE)
Ga = tn.Node(Ga2)
Gb1= np.append(ref[np.newaxis,np.newaxis,:],np.zeros((1,1,n)),axis=1).astype(self.DATA_TYPE)
Gb2= np.append(Gb1, np.append(np.zeros((1,1,n)), np.ones((1,1,n)),axis=1), axis=0).astype(self.DATA_TYPE)
Gb3 = np.array([[1,0],[0,-1]])
Gb = tn.Node(np.kron(Gb3, np.transpose(Gb2)))
Gc1 = np.append(np.ones((n,1)), np.repeat(ref[:,np.newaxis],3,axis=1), axis=1).astype(self.DATA_TYPE)
Gc = tn.Node(np.transpose(Gc1))
# bas1B = range(n)
Gabc = tn.ncon([Ga, Gb, Gc], [(-1, 1), (-2, 1, 2), (2, -3)])
# for a in bas1B:
# for b in bas1B:
# for c in bas1B:
# for d in bas1B:
# for e in bas1B:
# for f in bas1B:
# if not(a == d and b == e and c == f):
# tensor[a,b,c,d,e,f] = 0
# return tensor
# final = tn.outer_product(Gabc, tn.Node(np.ones((n,n,n)))).tensor
# @numba.jit(nopython=True)#, parallel=True)
# def enforce_delta(n, tensor):
# bas1B = range(n)
# for a in bas1B:
# for b in bas1B:
# for c in bas1B:
# for d in bas1B:
# for e in bas1B:
# for f in bas1B:
# if not(a == d and b == e and c == f):
# tensor[a,b,c,d,e,f] = 0
# return tensor
if flag == 0:
occC = tn.outer_product(Gabc, tn.Node(np.ones(n))) #tn.Node(enforce_delta(n, final))
if flag == 1:
occC = np.zeros((n,n,n),dtype=np.float32)
occC = tn.outer_product(Gabc, tn.Node(np.ones((n,n,n)))) #tn.Node(enforce_delta(n, final))
# occC = np.zeros((n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
for c in bas1B:
occC[a,b,c] = ref[a]*ref[b]*(1-ref[c]) + \
(1-ref[a])*(1-ref[b])*ref[c]
# for a in bas1B:
# for b in bas1B:
# for c in bas1B:
# occC[a,b,c] = ref[a]*ref[b]*(1-ref[c]) + \
# (1-ref[a])*(1-ref[b])*ref[c]
return occC
def __get_occD(self, flag=0):
......
Markdown is supported
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