Real Time Evolution of a Honeycomb Lattice

How do I use this algorithm? What does that parameter do?
Post Reply
jmasur
Posts: 3
Joined: 05 Apr 2022, 22:43

Real Time Evolution of a Honeycomb Lattice

Post by jmasur »

Hello,

I am trying to simulate a specific type of evolution for Hexagonal Boron-Nitrate which I have defined as follows:

Code: 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())
I have tried two ways to evolve this thing, first is evolving by the WPO1 method implemented in ExpMPOEvolution one step at a time. I first find the ground state using DMRG (this works), and then I try to evolve by my time dependent Hamiltonian (I do not use the time dependent version of WPO1 because I need to calculate expectation values after each step), this causes the norm to overflow as it increases exponentially with each time step. The minimum example for this is:

Code: 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))
where phi_tl is just\( \Phi(t) = \frac{aF_0}{\omega_0} \sin^2 \left(\omega_0 t / 20 \right) \sin (\omega_0 t) \).

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:

Code: 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))
In summary, using WPO leads to an overflow as the norm gets way too big, and simply creating the time evolution operator and applying it does not give the correct evolution. Any advice on how to fix this would be greatly appreciated.

Thanks,
Jacob
Post Reply