I am trying to simulate a specific type of evolution for Hexagonal Boron-Nitrate which I have defined as follows:
Python: Select all
class HBNHamiltonian(CouplingModel, MPOModel):
def __init__(self, model_params, phi):
self.options = model_params
model_params["phi"] = phi
p = model_params["p"]
t0 = p.t0
ub = p.ub
un = p.un
site = SpinHalfFermionSite()
lat = Honeycomb(p.nx, p.ny, site)
CouplingModel.__init__(self, lat)
# we assume a pulse polarized in the x direction
# nearest neighbor in the i direction
s1 = -t0 * np.exp(-1j * phi)
self.add_coupling(s1, 1, "Cdu", 0, "Cu", [1, 0], plus_hc=True)
self.add_coupling(s1, 1, "Cdd", 0, "Cd", [1, 0], plus_hc=True)
# nearest neighbor in the cos(pi/3)i + sin(pi/3)j direction
s2 = -t0 * np.exp(-1j * phi * .5)
self.add_coupling(s2, 0, "Cdu", 1, "Cu", [0, 0], plus_hc=True)
self.add_coupling(s2, 0, "Cdd", 1, "Cd", [0, 0], plus_hc=True)
# nearest neighbor in the cos(2pi/3)i + sin(2pi/3)j direction
s3 = -t0 * np.exp(1j * phi * .5)
self.add_coupling(s3, 1, "Cdu", 0, "Cu", [0, 1], plus_hc=True)
self.add_coupling(s3, 1, "Cdd", 0, "Cd", [0, 1], plus_hc=True)
# add onsite interaction
self.add_onsite(ub, 0, "NuNd")
self.add_onsite(un, 1, "NuNd")
MPOModel.__init__(self, lat, self.calc_H_MPO())
Python: Select all
# maximum bond dimension, used for both DMRG and TEBD, multiple of 200
maxdim = 1000
# number of unit cells in the x and y direction
nx = 2
ny = 2
# number of sites (unit cell has 2 sites)
N = nx * ny * 2
# the number of time steps
nsteps = 4000
# strength of pulse (MV/cm)
iF0 = 10.
maxerr = 1e-12 # used for DMRG
# hopping parameter, in units eV
it = 2.7
# lattice spacing, in angstroms
ia = 2.5
# pulse parameters
iomega0 = 32.9 # driving (angular) frequency, in THz
cycles = 10
pbc = False # periodic boundary conditions
ub = 3.3 # onsite potential for boron in eV
un = -1.7 # onsite potential for nitrogen in eV
# used just to scale parameters
p = Parameters(N, 0, it, ia, cycles, iomega0, iF0, pbc, nx=nx, ny=ny, ub=ub, un=un)
gps = dict(p=p)
model = HBNHamiltonian(gps, 0.)
sites = model.lat.mps_sites()
state = ["up", "down"] * (N // 2)
psi0_i = MPS.from_product_state(sites, state)
# run DMRG to find the ground state
chi_list = {0:20, 1:40, 2:100, 4:200}
chi = 200
iter = 4
while chi < maxdim:
iter += 2
chi += 200
chi_list[iter] = chi
dmrg_dict = {"chi_list":chi_list, "max_E_err":maxerr, "max_sweeps":(maxdim / 100) + 4, "mixer":True, "combine":False}
dmrg_params = Config(dmrg_dict, "DMRG-maxerr{}".format(maxerr))
dmrg = DMRG(psi0_i, model, dmrg_params)
E, psi0 = dmrg.run()
psi = psi0
ti = 0
tf = 2 * np.pi * cycles / p.field
times, delta = np.linspace(ti, tf, num=nsteps, endpoint=True, retstep=True)
# set up additional operators and expectation values
currentop = HBNCurrent(gps, 0.)
phis = [0.]
energies = [model.H_MPO.expectation_value(psi)]
currents = [currentop.H_MPO.expectation_value(psi)]
# evolve 1 step at a time with WPO1 method
for i, t in enumerate(times[1:]):
# apply U to psi
wpo_dict = {"dt":delta, "order":1, "start_time":t, "start_trunc_err":TruncationError(eps=maxerr), "trunc_params":{"svd_min":maxerr, "chi_max":maxdim}, "N_steps":1, "F0":iF0, "compression_method":"SVD"}
wpo_params = Config(wpo_dict, "HBN")
evolver = ExpMPOEvolution(psi, model, wpo_params)
psi = evolver.run()
psi.canonical_form_finite(renormalize=True, cutoff=maxerr)
# calculate next phi and operators
phi = phi_tl(p, t)
gps = dict(p=p)
model = HBNHamiltonian(gps, phi)
currentop = HBNCurrent(gps, phi)
# obtain expectation values
energies.append(model.H_MPO.expectation_value(psi))
currents.append(currentop.H_MPO.expectation_value(psi))
I tried just doing a simple version of time evolution where I just create the operator using make_U and apply it to the state, but this did not give the desired results. I tested this code with a slightly modified Fermi-Hubbard model for which I have obtained the correct evolution using TEBD:
Python: Select all
maxerr = 1e-12
maxdim = 1000
N = 10
cycles = 10
it = .52
iU = 0. * it
a = 4
iomega0 = 32.9
iF0 = 10.
pbc = False
nsteps = 4000
p = Parameters(N, iU, it, a, cycles, iomega0, iF0, pbc)
model = FHHamiltonian(p, 0)
sites = model.lat.mps_sites()
state = ["up", "down"] * (N // 2)
psi0_i = MPS.from_product_state(sites, state)
# the max bond dimension
chi_list = {0:20, 1:40, 2:100, 4:200}
chi = 200
iter = 4
while chi < maxdim:
iter += 2
chi += 200
chi_list[iter] = chi
dmrg_dict = {"chi_list":chi_list, "max_E_err":maxerr, "max_sweeps":(maxdim / 100) + 4, "mixer":True, "combine":False}
dmrg_params = Config(dmrg_dict, "DMRG-maxerr{}".format(maxerr))
dmrg = DMRG(psi0_i, model, dmrg_params)
E, psi0 = dmrg.run()
psi = psi0
ti = 0
tf = 2 * np.pi * cycles / p.field
times, delta = np.linspace(ti, tf, num=nsteps, endpoint=True, retstep=True)
dt = -1j * delta
start_time = time.time()
ham = model
current = FHCurrent(p, 0.)
energies = [ham.H_MPO.expectation_value(psi)]
currents = [current.H_MPO.expectation_value(psi)]
phis = [0.]
Udt = None
for i, t in enumerate(times[1:]):
# make and apply the time evolution operator
del Udt
Udt = ham.H_MPO.make_U_II(dt)
Udt.apply_naively(psi)
psi.canonical_form_finite(renormalize=True, cutoff=maxerr)
print(psi.norm)
# make the new operators at this time
del ham
del current
phi = phi_tl(p, t)
ham = FHHamiltonian(p, phi)
current = FHCurrent(p, phi)
# calculate expectation values
phis.append(phi)
energies.append(ham.H_MPO.expectation_value(psi))
currents.append(current.H_MPO.expectation_value(psi))
Thanks,
Jacob