Commit d53e39a6 authored by Jacob Davison's avatar Jacob Davison
Browse files

precision from floats to doubles

parent 74df6285
......@@ -438,24 +438,7 @@ if __name__ == '__main__':
#test_refs('logs_refs\\')
# test = main(4,4)
import os
print(os.environ)
# # test_exact('plots_exact\\')
# # print(ci_matrix.exact_diagonalization(1.0, 0.5,0.0))
# tracemalloc.start()
# for i in range(5):
# test = main(4,4)
# print(test[-1])
# snapshot = tracemalloc.take_snapshot()
# top_stats = snapshot.statistics('lineno')
# total = sum(stat.size for stat in top_stats)
# print("Total allocated size: %.1f KiB" % (total / 1024))
# snapshot = tracemalloc.take_snapshot()
# top_stats = snapshot.statistics('lineno')
# total = sum(stat.size for stat in top_stats)
# print("Final allocated size: %.1f KiB" % (total / 1024))
......@@ -73,10 +73,10 @@ def unravel(E, f, G, W):
Returns:
concatenation of tensors peeled into 1D arrays"""
unravel_E = np.reshape(E, -1)
unravel_f = np.reshape(f, -1)
unravel_G = np.reshape(G, -1)
unravel_W = np.reshape(W, -1)
unravel_E = np.reshape(E, -1).astype(np.float32)
unravel_f = np.reshape(f, -1).astype(np.float32)
unravel_G = np.reshape(G, -1).astype(np.float32)
unravel_W = np.reshape(W, -1).astype(np.float32)
return np.concatenate([unravel_E, unravel_f, unravel_G, unravel_W], axis=0)
......@@ -96,14 +96,14 @@ def ravel(y, bas_len):
# bas_len = len(np.append(holes,particles))
ravel_E = np.reshape(y[0], ())
ravel_f = np.reshape(y[1:bas_len**2+1], (bas_len, bas_len))
ravel_E = np.reshape(y[0], ()).astype(np.float32)
ravel_f = np.reshape(y[1:bas_len**2+1], (bas_len, bas_len)).astype(np.float32)
ravel_G = np.reshape(y[bas_len**2+1:bas_len**2+1+bas_len**4],
(bas_len, bas_len, bas_len, bas_len))
(bas_len, bas_len, bas_len, bas_len)).astype(np.float32)
ravel_W = np.reshape(y[bas_len**2+1+bas_len**4:bas_len**2+1+bas_len**4+bas_len**6],
(bas_len,bas_len,bas_len,bas_len,bas_len,bas_len))
(bas_len,bas_len,bas_len,bas_len,bas_len,bas_len)).astype(np.float32)
return(ravel_E, ravel_f, ravel_G, ravel_W)
return (ravel_E, ravel_f, ravel_G, ravel_W)
def main3b(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0):
"""Main method uses scipy.integrate.ode to solve the IMSRG3 flow
......
......@@ -233,7 +233,7 @@ class Flow_IMSRG3(Flow_IMSRG2):
eta1B = partition[0]
eta2B = partition[1]
eta3B = partition[2]
# print(eta3B.shape)
occA = self._occA
occA2 = self._occA2
occB = self._occB
......@@ -260,13 +260,13 @@ class Flow_IMSRG3(Flow_IMSRG2):
sum4_1b = sum4_1b_3 - sum4_1b_4
# fifth term
sum5_1b_1 = ncon([occF, eta3B],
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],
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],
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],
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
......
......@@ -30,6 +30,7 @@ class WegnerGenerator(Generator):
self.f = h.f
self.G = h.G
self._holes = h.holes
self._particles = h.particles
......@@ -81,12 +82,12 @@ class WegnerGenerator(Generator):
particles = self._particles
# - Decouple off-diagonal 1B and 2B pieces
fod = np.zeros(f.shape)
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)]
fd = f - fod
God = np.zeros(G.shape)
God = np.zeros(G.shape, dtype=np.float32)
God[np.ix_(particles, particles, holes, holes)] += G[np.ix_(particles, particles, holes, holes)]
God[np.ix_(holes, holes, particles, particles)] += G[np.ix_(holes, holes, particles, particles)]
Gd = G - God
......@@ -117,6 +118,7 @@ class WegnerGenerator(Generator):
occB = self._occB
occC = self._occC
# - Calculate 1B generator
# first term
sum1_1b_1 = ncon([fd, fod], [(-1, 0), (0, -2)]).numpy()
......@@ -201,7 +203,7 @@ class WegnerGenerator3B(WegnerGenerator):
self.f = h.f
self.G = h.G
self.W = np.zeros(h.n_sp_states*np.ones(6,dtype=np.int64))
self.W = np.zeros(h.n_sp_states*np.ones(6,dtype=np.int32),dtype=np.float32)
self._holes = h.holes
self._particles = h.particles
......@@ -252,13 +254,13 @@ class WegnerGenerator3B(WegnerGenerator):
fod = partition[1]
Gd = partition[2]
God = partition[3]
W = self.W
holes = self._holes
particles = self._particles
Wod = np.zeros(W.shape)
Wod = np.zeros(W.shape, dtype=np.float32)
Wod[np.ix_(particles, particles, particles, holes, holes, holes)] += W[np.ix_(particles, particles, particles, holes, holes, holes)]
Wod[np.ix_(holes, holes, holes, particles, particles, particles)] += W[np.ix_(holes, holes, holes, particles, particles, particles)]
Wd = W - Wod
......@@ -285,7 +287,7 @@ class WegnerGenerator3B(WegnerGenerator):
God = partition[3]
Wd = partition[4]
Wod = partition[5]
eta1B, eta2B = super().calc_eta()
holes = self._holes
......@@ -302,6 +304,7 @@ class WegnerGenerator3B(WegnerGenerator):
occI = self._occI
occJ = self._occJ
# Calculate 1B generator
# fourth term
sum4_1b_1 = np.matmul(np.transpose(occD,[2,3,0,1]), God)
......@@ -311,13 +314,13 @@ class WegnerGenerator3B(WegnerGenerator):
sum4_1b = sum4_1b_3 - sum4_1b_4
# fifth term
sum5_1b_1 = ncon([occF, Wd],
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],
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],
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],
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
......@@ -375,8 +378,8 @@ class WegnerGenerator3B(WegnerGenerator):
np.transpose(sum1_3b_7, [0,1,2,5,4,3])
#terms with P(ij/k)P(l/mn) -- line 3
sum1_3b_9 = ncon([Gd, God], [(-1,-2,-4,0),(0,-3,-5,-6)])
sum1_3b_10 = ncon([God, Gd], [(-1,-2,-4,0),(0,-3,-5,-6)])
sum1_3b_9 = ncon([Gd, God], [(-1,-2,-4,0),(0,-3,-5,-6)]).numpy()
sum1_3b_10 = ncon([God, Gd], [(-1,-2,-4,0),(0,-3,-5,-6)]).numpy()
sum1_3b_11 = sum1_3b_9 - sum1_3b_10
sum1_3b_12 = sum1_3b_11 - np.transpose(sum1_3b_11, [0,1,2,4,3,5]) - \
np.transpose(sum1_3b_11, [0,1,2,5,4,3])
......@@ -424,5 +427,4 @@ class WegnerGenerator3B(WegnerGenerator):
eta3B = sum1_3b + 0.5*sum4_3b + (-0.5)*sum5_3b + (-1.0)*sum6_3b + (1/6)*sum7_3b + 0.5*sum8_3b
return (eta1B, eta2B, eta3B)
......@@ -40,9 +40,9 @@ class PairingHamiltonian2B(Hamiltonian):
else:
self._reference = ref
self._holes = np.arange(n_hole_states, dtype=np.int64)
self._holes = np.arange(n_hole_states, dtype=np.int32)
self._n_sp_states = n_hole_states + n_particle_states
self._particles = np.arange(n_hole_states,self.n_sp_states, dtype=np.int64)
self._particles = np.arange(n_hole_states,self.n_sp_states, dtype=np.int32)
self._sp_basis = np.append(self.holes, self.particles)
self._H1B, self._H2B = self.construct()
......@@ -217,7 +217,7 @@ class PairingHamiltonian2B(Hamiltonian):
# two body part of Hamiltonian constructed from four indices
# with non-zero elements defined by pairing term
H2B = np.zeros(np.ones(4, dtype=np.int64)*self.n_sp_states)
H2B = np.zeros(np.ones(4, dtype=np.int32)*self.n_sp_states,dtype=np.float32)
for p in bas1B:
for q in bas1B:
for r in bas1B:
......@@ -260,7 +260,7 @@ class PairingHamiltonian2B(Hamiltonian):
tb_contract = 0.5*net.contract(flatten).tensor.numpy()
E = ob_contract + tb_contract
E = E.astype(np.float32)
# - Calculate 1B piece
ob_node1b = net.add_node(H1B_t)
......@@ -270,7 +270,7 @@ class PairingHamiltonian2B(Hamiltonian):
tb_contract = net.contract(tb_ihjh)
f = ob_node1b.tensor.numpy() + tb_contract.tensor.numpy()
f = f.astype(np.float32)
G = H2B_t
del net
......
......@@ -130,14 +130,14 @@ class OccupationTensors(object):
n = len(bas1B)
if flag == 0: # default
occA = np.zeros((n,n,n,n))
occA = np.zeros((n,n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
occA[a,b,a,b] = ref[a] - ref[b]
if flag == 1:
occA = np.zeros((n,n))
occA = np.zeros((n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
......@@ -161,14 +161,14 @@ class OccupationTensors(object):
n = len(bas1B)
if flag == 0: # default
occB = np.zeros((n,n,n,n))
occB = np.zeros((n,n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
occB[a,b,a,b] = 1 - ref[a] - ref[b]
if flag == 1:
occB = np.zeros((n,n))
occB = np.zeros((n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
......@@ -192,7 +192,7 @@ class OccupationTensors(object):
n = len(bas1B)
if flag == 0: # default
occC = np.zeros((n,n,n,n,n,n))
occC = np.zeros((n,n,n,n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
......@@ -202,7 +202,7 @@ class OccupationTensors(object):
if flag == 1:
occC = np.zeros((n,n,n))
occC = np.zeros((n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
......@@ -227,7 +227,7 @@ class OccupationTensors(object):
n = len(bas1B)
if flag == 0: # default
occD = np.zeros((n,n,n,n,n,n,n,n))
occD = np.zeros((n,n,n,n,n,n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
......@@ -237,7 +237,7 @@ class OccupationTensors(object):
(1-ref[c])*(1-ref[d])
if flag == 1:
occD = np.zeros((n,n,n,n))
occD = np.zeros((n,n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
......@@ -260,7 +260,7 @@ class OccupationTensors(object):
ref = self._reference
n = len(bas1B)
occE = np.zeros((n,n,n,n,n,n))
occE = np.zeros((n,n,n,n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
......@@ -308,7 +308,7 @@ class OccupationTensors(object):
ref = self._reference
n = len(bas1B)
occG = np.zeros((n,n,n,n,n,n))
occG = np.zeros((n,n,n,n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
......@@ -328,7 +328,7 @@ class OccupationTensors(object):
ref = self._reference
n = len(bas1B)
occH = np.zeros((n,n,n,n,n,n,n,n))
occH = np.zeros((n,n,n,n,n,n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
......@@ -350,7 +350,7 @@ class OccupationTensors(object):
ref = self._reference
n = len(bas1B)
occI = np.zeros((n,n,n,n,n,n,n,n))
occI = np.zeros((n,n,n,n,n,n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
......@@ -371,7 +371,7 @@ class OccupationTensors(object):
ref = self._reference
n = len(bas1B)
occJ = np.zeros((n,n,n,n,n,n))
occJ = np.zeros((n,n,n,n,n,n),dtype=np.float32)
for a in bas1B:
for b in bas1B:
......
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