Commit bcc35662 authored by Davison's avatar Davison
Browse files

small changes here and there; updated docstrings

parent 9b54dec2
This diff is collapsed.
This diff is collapsed.
from tensornetwork import *
import numpy as np
from hamiltonian import PairingHamiltonian
from occupation_tensors import OccupationTensors
from generator import WegnerGenerator
from hamiltonian import *
from occupation_tensors import *
from generator import *
class Flow_IMSRG2(object):
class Flow(object):
"""Parent class for organization purposes. Ideally, all Flow
classes should inherit from this class. In this way, AssertionErrors
can be handled in a general way."""
def flow():
print("Function that iterates flow equation once")
class Flow_IMSRG2(Flow):
"""Calculates the flow equations for the IMSRG(2)."""
def __init__(self, h, occ_t):
assert isinstance(h, PairingHamiltonian), "Arg 0 must be PairingHamiltonian object"
"""Class constructor. Instantiates Flow_IMSRG2 object.
Arguments:
h -- Hamiltonian object
occ_t -- OccupationTensors object"""
assert isinstance(h, Hamiltonian), "Arg 0 must be PairingHamiltonian object"
assert isinstance(occ_t, OccupationTensors), "Arg 1 must be OccupationTensors object"
# self.f = h.f
......@@ -16,10 +32,10 @@ class Flow_IMSRG2(object):
self._holes = h.holes
self._particles = h.particles
# self._occA = occ_t.occA
# self._occB = occ_t.occB
# self._occC = occ_t.occC
# self._occD = occ_t.occD
self._occA = occ_t.occA
self._occB = occ_t.occB
self._occC = occ_t.occC
self._occD = occ_t.occD
# @property
# def f(self):
......@@ -37,26 +53,37 @@ class Flow_IMSRG2(object):
# def G(self, G):
# self._G = G
def flow(self, gen, occ_t):
assert isinstance(gen, WegnerGenerator), "Arg 0 must be WegnerGenerator object"
def flow(self, gen):
"""Iterates the IMSRG2 flow equations once.
Arugments:
gen -- Generator object; generator produces the flow
Returns:
(dE, -- zero-body tensor
df, -- one-body tensor
dG) -- two-body tensor"""
assert isinstance(gen, Generator), "Arg 0 must be WegnerGenerator object"
# f = self.f
# G = self.G
f = gen.f
G = gen.G
eta1B = gen.calc_eta1B()
eta2B = gen.calc_eta2B()
eta1B, eta2B = gen.calc_eta()
occA = occ_t.occA
occB = occ_t.occB
occC = occ_t.occC
occD = occ_t.occD
# occA = occ_t.occA
# occB = occ_t.occB
# occC = occ_t.occC
# occD = occ_t.occD
# occA = self._occA
# occB = self._occB
# occC = self._occC
# occD = self._occD
occA = self._occA
occB = self._occB
occC = self._occC
occD = self._occD
# - Calculate dE/ds
# first term
......
from tensornetwork import *
import numpy as np
from hamiltonian import PairingHamiltonian
from occupation_tensors import OccupationTensors
from hamiltonian import *
from occupation_tensors import *
class WegnerGenerator(object):
class Generator(object):
"""Parent class for organization purposes. Ideally, all Generator
classes should inherit from this class. In this way, AssertionErrors
can be handled in a general way."""
def calc_eta():
print("Function that calculates the generator")
class WegnerGenerator(Generator):
"""Calculate Wegner's generator for a normal ordered Hamiltonian."""
def __init__(self, h, occ_t):
"""Class constructor. Instantiate WegnerGenerator object.
Arguments:
h -- Hamiltonian object (must be normal-ordered)
occ_t -- OccupationTensor object"""
assert isinstance(h, PairingHamiltonian), "Arg 0 must be Hamiltonian object"
assert isinstance(h, Hamiltonian), "Arg 0 must be Hamiltonian object"
assert isinstance(occ_t, OccupationTensors), "Arg 1 must be OccupationTensors object"
self.f = h.f
......@@ -21,28 +35,42 @@ class WegnerGenerator(object):
self._occC = occ_t.occC
self._occD = occ_t.occD
# self._eta1B = self.eta1B
# self._eta2B = self.eta2B
@property
def f(self):
"""Returns:
f -- one-body tensor elements (initialized by Hamiltonian object)"""
return self._f
@property
def G(self):
"""Returns:
f -- two-body tensor elements (initialized by Hamiltonian object)"""
return self._G
@f.setter
def f(self, f):
"""Sets the one-body tensor."""
self._f = f
@G.setter
def G(self, G):
"""Sets the two-body tensor."""
self._G = G
def __decouple_OD(self):
"""Decouple the off-/diagonal elements from each other in
the one- and two-body tensors. This procedure is outlined in
An Advanced Course in Computation Nuclear Physics, Ch.10.
Returns:
(fd, -- diagonal part of f
fod, -- off-diagonal part of f
Gd, -- diagonal part of G
God) -- off-diagonal part of G"""
f = self.f
G = self.G
holes = self._holes
......@@ -61,13 +89,23 @@ class WegnerGenerator(object):
return (fd, fod, Gd, God)
def calc_eta1B(self):
def calc_eta(self):
"""Calculate the generator. The terms are defined in An
Advanced Course in Computation Nuclear Physics, Ch.10.
See also dx.doi.org/10.1016/j.physrep.2015.12.007
Returns:
(eta1B, -- one-body generator
eta2B) -- two-body generator"""
fd, fod, Gd, God = self.__decouple_OD()
holes = self._holes
particles = self._particles
occA = self._occA
occB = self._occB
occC = self._occC
# - Calculate 1B generator
......@@ -90,17 +128,6 @@ class WegnerGenerator(object):
eta1B = sum1_1b + sum2_1b + 0.5*sum3_1b
return eta1B
def calc_eta2B(self):
fd, fod, Gd, God = self.__decouple_OD()
holes = self._holes
particles = self._particles
occA = self._occA
occB = self._occB
# - Calculate 2B generator
# first term (P_ij piece)
sum1_2b_1 = ncon([fd, God], [(-1, 0), (0, -2, -3, -4)]).numpy()
......@@ -135,7 +162,84 @@ class WegnerGenerator(object):
eta2B = sum1_2b + 0.5*sum2_2b + sum3_2b
return eta2B
return (eta1B, eta2B)
# def calc_eta1B(self):
# fd, fod, Gd, God = self.__decouple_OD()
# holes = self._holes
# particles = self._particles
# occA = self._occA
# occC = self._occC
# # - Calculate 1B generator
# # first term
# sum1_1b_1 = ncon([fd, fod], [(-1, 0), (0, -2)]).numpy()
# sum1_1b_2 = np.transpose(sum1_1b_1)
# sum1_1b = sum1_1b_1 - sum1_1b_2
# # second term
# sum2_1b_1 = ncon([fd, God], [(0, 1), (1, -1, 0, -2)]).numpy()
# sum2_1b_2 = ncon([fod, Gd], [(0, 1), (1, -1, 0, -2)]).numpy()
# sum2_1b_3 = sum2_1b_1 - sum2_1b_2
# sum2_1b = ncon([occA, sum2_1b_3],[(-1, -2, 0, 1), (0,1)]).numpy()
# # third term
# sum3_1b_1 = ncon([occC, God], [(-1, -2, -3, 0, 1, 2), (0, 1, 2, -4)]).numpy()
# sum3_1b_2 = ncon([Gd, sum3_1b_1], [(2, -1, 0, 1), (0, 1, 2, -2)]).numpy()
# sum3_1b_3 = np.transpose(sum3_1b_2)
# sum3_1b = sum3_1b_2 - sum3_1b_3
# eta1B = sum1_1b + sum2_1b + 0.5*sum3_1b
# return eta1B
# def calc_eta2B(self):
# fd, fod, Gd, God = self.__decouple_OD()
# holes = self._holes
# particles = self._particles
# occA = self._occA
# occB = self._occB
# # - Calculate 2B generator
# # first term (P_ij piece)
# sum1_2b_1 = ncon([fd, God], [(-1, 0), (0, -2, -3, -4)]).numpy()
# sum1_2b_2 = ncon([fod, Gd], [(-1, 0), (0, -2, -3, -4)]).numpy()
# sum1_2b_3 = sum1_2b_1 - sum1_2b_2
# sum1_2b_4 = np.transpose(sum1_2b_3, [1, 0, 2, 3])
# sum1_2b_5 = sum1_2b_3 - sum1_2b_4
# # first term (P_kl piece)
# sum1_2b_6 = ncon([fd, God], [(0, -3), (-1, -2, 0, -4)]).numpy()
# sum1_2b_7 = ncon([fod, Gd], [(0, -3), (-1, -2, 0, -4)]).numpy()
# sum1_2b_8 = sum1_2b_6 - sum1_2b_7
# sum1_2b_9 = np.transpose(sum1_2b_8, [0, 1, 3, 2])
# sum1_2b_10 = sum1_2b_8 - sum1_2b_9
# sum1_2b = sum1_2b_5 - sum1_2b_10
# # second term
# sum2_2b_1 = ncon([occB, God], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
# sum2_2b_2 = ncon([occB, Gd], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
# sum2_2b_3 = ncon([Gd, sum2_2b_1], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
# sum2_2b_4 = ncon([God, sum2_2b_2], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
# sum2_2b = sum2_2b_3 - sum2_2b_4
# # third term
# sum3_2b_1 = ncon([Gd, God], [(0, -1, 1, -3), (1, -2, 0, -4)]).numpy()
# sum3_2b_2 = np.transpose(sum3_2b_1, [1, 0, 2, 3])
# sum3_2b_3 = np.transpose(sum3_2b_1, [0, 1, 3, 2])
# sum3_2b_4 = np.transpose(sum3_2b_1, [1, 0, 3, 2])
# sum3_2b_5 = sum3_2b_1 - sum3_2b_2 - sum3_2b_3 + sum3_2b_4
# sum3_2b = ncon([occA, sum3_2b_5], [(0, 1, -1, -2), (0, 1, -3, -4)]).numpy()
# eta2B = sum1_2b + 0.5*sum2_2b + sum3_2b
# return eta2B
# @property
# def eta1B(self):
......
import numpy as np
from tensornetwork import *
class PairingHamiltonian(object):
class Hamiltonian(object):
"""Parent class for organization purposes. Ideally, all Hamiltonian
classes should inherit from this class. In this way, AssertionErrors
can be handled in a general way."""
def __construct():
print("Function that constructs the Hamiltonian")
def __normal_order():
print("Function that normal-orders the Hamiltonian")
class PairingHamiltonian2B(Hamiltonian):
"""Generate the two-body pairing Hamiltonian. Inherits from Hamiltonian."""
def __init__(self, n_hole_states, n_particle_states, d=1.0, g=0.5, pb=0.0):
"""Class constructor. Instantiate PairingHamiltonian2B object.
Arguments:
n_hole_states -- number of holes states in the single particle basis
n_particle_states -- number of particles states in the single particle basis
Keyword arguments:
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)"""
def __init__(self, n_hole_states, n_particle_states, d=1, g=0.5, pb=0.0):
self._d = d
self._g = g
self._pb = pb
......@@ -19,58 +44,105 @@ class PairingHamiltonian(object):
@property
def d(self):
"""Returns:
d -- energy level spacing."""
return self._d
@property
def g(self):
"""Returns:
g -- pairing strength."""
return self._g
@property
def pb(self):
"""Returns:
pb -- pair-breaking strength."""
return self._pb
@property
def reference(self):
"""Returns:
reference -- reference state (ground state)."""
return self._reference
@property
def holes(self):
"""Returns:
holes -- indices of hole states in single particle basis."""
return self._holes
@property
def particles(self):
"""Returns:
particles -- indices of particle states in single particle basis."""
return self._particles
@property
def sp_basis(self):
"""Returns:
sp_basis -- indices of full single particle basis."""
return self._sp_basis
@property
def n_sp_states(self):
"""Returns:
n_sp_states -- size of single-particle basis."""
return self._n_sp_states
@property
def H1B(self):
"""Returns:
H1B -- one-body (rank 2) tensor defined by __construct()."""
return self._H1B
@property
def H2B(self):
"""Returns:
H2B -- two-body (rank 4) tensor defined by __construct()."""
return self._H2B
@property
def E(self):
"""Returns:
E -- zero-body (rank 0) tensor defined by __normal_order()."""
return self._E
@property
def f(self):
"""Returns:
f -- one-body (rank 2) tensor defined by __normal_order()."""
return self._f
@property
def G(self):
"""Returns:
G -- two-body (rank 4) tensor defined by __normal_order()."""
return self._G
def __delta2B(self, p,q,r,s):
"""Determines if a two-body tensor elements should be zero,
positive, or negative. This behavior is dicated by the pairing
term in pairing Hamiltonian.
Arguments:
p,q,r,s -- indices in single-particle basis"""
pp = np.floor_divide(p,2)
qp = np.floor_divide(q,2)
rp = np.floor_divide(r,2)
......@@ -93,6 +165,14 @@ class PairingHamiltonian(object):
return 0
def __deltaPB(self, p,q,r,s):
"""Determines if a two-body tensor elements should be zero,
positive, or negative. This behavior is dicated by the pair-
breaking term in pairing Hamiltonian.
Arguments:
p,q,r,s -- indices in single particle basis"""
pp = np.floor_divide(p,2)
qp = np.floor_divide(q,2)
rp = np.floor_divide(r,2)
......@@ -114,6 +194,15 @@ class PairingHamiltonian(object):
return 0
def __construct(self):
"""Constructs the one- and two-body pieces of the pairing
Hamiltonian.
Returns:
(H1B, -- one-body tensor elements (defined by one-body operator)
H2B) -- two-body tensor elements (defined by two-body operator)"""
bas1B = self.sp_basis # get the single particle basis
# one body part of Hamiltonian is floor-division of basis index
......@@ -133,6 +222,14 @@ class PairingHamiltonian(object):
return (H1B, H2B)
def __normal_order(self):
"""Normal-orders the pairing Hamiltonian.
Returns:
(E, -- zero-body piece
f, -- one-body piece
G) -- two-body piece"""
bas1B = self.sp_basis # get the single particle basis
H1B_t = self.H1B # get the 1B tensor
H2B_t = self.H2B # get the 2B tensor
......@@ -172,4 +269,9 @@ class PairingHamiltonian(object):
return (E, f, G)
class PairingHamiltonian3B(PairingHamiltonian2B):
def __delta3B():
pass
pass
from hamiltonian import PairingHamiltonian
from occupation_tensors import OccupationTensors
from generator import WegnerGenerator
from flow import Flow_IMSRG2
# Main program for IM-SRG.
# Author: Jacob Davison
# Date: 07/10/2019
# import packages, libraries, and modules
from hamiltonian import *
from occupation_tensors import *
from generator import *
from flow import *
from scipy.integrate import odeint, ode
import numpy as np
import time
ha = PairingHamiltonian(4,4)
ha = PairingHamiltonian2B(4,4)
ot = OccupationTensors(ha.sp_basis, ha.reference)
wg = WegnerGenerator(ha, ot)
fl = Flow_IMSRG2(ha, ot)
def derivative(t, y, hamiltonian, occ_tensors, generator, flow):
assert isinstance(hamiltonian, PairingHamiltonian), "Arg 2 must be Hamiltonian object"
"""Defines the derivative to pass into ode object.
Arguments:
(required by scipy.integrate.ode)
t -- points at which to solve for y
y -- in this case, 1D array that contains E, f, G
(additional parameters)
hamiltonian -- Hamiltonian object
occ_tensors -- OccupationTensors object
generator -- Generator object
flow -- Flow object
Returns:
dy -- next step in flow"""
assert isinstance(hamiltonian, Hamiltonian), "Arg 2 must be Hamiltonian object"
assert isinstance(occ_tensors, OccupationTensors), "Arg 3 must be OccupationTensors object"
assert isinstance(generator, WegnerGenerator), "Arg 4 must be WegnerGenerator object"
assert isinstance(flow, Flow_IMSRG2), "Arg 5 must be a Flow_IMSRG2 object"
assert isinstance(generator, Generator), "Arg 4 must be Generator object"
assert isinstance(flow, Flow), "Arg 5 must be a Flow object"
E, f, G = ravel(y, hamiltonian.holes, hamiltonian.particles)
E, f, G = ravel(y, hamiltonian.n_sp_states)
generator.f = f
generator.G = G
dE, df, dG = flow.flow(generator, occ_tensors)
dE, df, dG = flow.flow(generator)
dy = unravel(dE, df, dG)
return dy
def unravel(E, f, G):
"""Transforms E, f, and G into a 1D array. Facilitates
compatability with scipy.integrate.ode.
Arguments:
E, f, G -- normal-ordered pieces of Hamiltonian
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)
return np.concatenate([unravel_E, unravel_f, unravel_G], axis=0)
def ravel(y, holes, particles):
def ravel(y, bas_len):
"""Transforms 1D array into E, f, and G. Facilitates
compatability with scipy.integrate.ode.
Arugments:
y -- 1D data array (output from unravel)
bas_len -- length of single particle basis
Returns:
E, f, G -- normal-ordered pieces of Hamiltonian"""
bas_len = len(np.append(holes,particles))
# 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))
......@@ -48,7 +93,10 @@ def ravel(y, holes, particles):
def main():
"""Main method uses scipy.integrate.ode to solve the IMSRG flow
equations."""
start = time.time()
print("""Pairing model IM-SRG flow:
d = {:2d}
g = {:2.4f}
......@@ -77,7 +125,7 @@ def main():
while solver.successful() and solver.t < sfinal:
ys = solver.integrate(sfinal, step=True)
Es, fs, Gs = ravel(ys, ha.holes, ha.particles)
Es, fs, Gs = ravel(ys, ha.n_sp_states)
s_vals.append(solver.t)
E_vals.append(Es)
......@@ -93,7 +141,11 @@ def main():
if len(E_vals) > 20 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