+1 vote
asked by (130 points)
edited by

I have been trying to recreate the dynamics of a driven Fermi-Hubbard Hamiltonian in Julia:

\hat{H}=-t_0 \sum_{j,\sigma}\left( e^{-i\Phi(t)}\hat{c}^{\dag}_{j,\sigma}\hat{c}{j+1,\sigma} + h.c. \right)+\sum_{j} \hat{n}_{j,\uparrow}\hat{n}{j,\downarrow}
(don't know why this isn't formatting this is the Hamiltonian given by equation 12 in https://journals.aps.org/pra/abstract/10.1103/PhysRevA.101.053408 )

where @@ \sigma @@ is spin, but my first step is attempting this with @@\Phi(t)@@ set to 0, which results in an undriven, time independent model. I have successfully performed DMRG to obtain the ground state (which I know because I have compared energies to exact calculations in the QuSpin package), but I have been unsuccessful in using TEBD (and RK2 and RK4) to evolve the ground state using the time independent Hamiltonian. I know applying the time evolution operator as an exponential of the TI Hamiltonian to the ground state should preserve the ground state, but instead the energy of the state resulting from the application of the time evolution operator increases throughout the evolution. Though I have not tested dependence rigorously, the energy increase seems to be mostly independent of the time step even though it should be @@\mathcal{O}(\delta t)@@ according to https://www.sciencedirect.com/science/article/pii/S0003491619302532. I have implemented a dynamic time step as well, and tested a couple different maximum bond dimensions in DMRG, as well as decreasing the cutoff for both DMRG and the apply function. Any help is greatly appreciated, thanks in advance! My code is below:

Use second order TEBD to evolve psi(t) to psi (t + dt)
function TEBD(psi, dt, time, params, independent, cutoff)

    phi = phi_tl(time + dt / 2, params, independent)

    # odd gates (1,2),(3,4),(5,6),...
    ogates = ITensor[]
    # even gates (2,3),(4,5),(6,7),...
    egates = ITensor[]
    for j=1:params.nsites
        s1 = params.space[j]
        #periodic BC
        if j == params.nsites
            s2 = params.space[1]
            s2 = params.space[j+1]
        # we have to define the two site operator so contributions
        # of each site will be counted twice
        ul = ur = params.U / 2
        # exp(i phi(t)) and exp(-i phi(t))
        eiphi = exp(1.0im * phi)
        eiphiconj = conj(eiphi)
        # create operator (only scaled parameters are passed in so t is always 1)
        hj = -eiphiconj * op("Cdagup",s1) * op("Cup",s2) +
             -eiphiconj * op("Cdagdn",s1) * op("Cdn",s2) +
             -eiphi * op("Cdagup",s2) * op("Cup",s1) +
             -eiphi * op("Cdagdn",s2) * op("Cdn",s1) +
             ul * op("Nupdn", s1) * op("Id",s2) +
             ur * op("Id",s1) * op("Nupdn", s2)
        # odd gate
        if j % 2 == 1
            Gj = exp(-1.0im * dt * hj)
            push!(ogates, Gj)
        # even gate
            Gj = exp(-1.0im * dt / 2 * hj)
            push!(egates, Gj)

    gates = ITensor[]
    append!(gates, egates)
    append!(gates, ogates)
    append!(gates, egates)

    return apply(gates, psi; cutoff=cutoff)


Propogates ground from 0 to tf using method to propogate by an adaptive timestep
    ground - ground state (MPS)
    tf - final time
    method - the function used to propogate a wavefunction over single dt
    dti - initial guess for timestep
    epsilon - total change allowed in a single time step
    independent - boolean that is True if we are evolving with a time independent
              Hamiltonian, and False otherwise
function propogation(ground::MPS, params, tf, method, dti, epsilon, independent, cutoff)
    # error measures (en is for current step, en1 is for previous step)
    en = en1 = 0.0
    time = 0.0
    # copy states
    psi = deepcopy(ground)
    # set dt and previous dt
    dt = pdt = dti

    # initialize vectors for saving times and expectation values
    times = [0]
    energies = [inner(ground, get_ham(0, params, independent), ground)]
    currents = [inner(ground, get_current(0, params, independent), ground)]

    # run for the entire time interval
    while time < tf

        # get the MPS at time + dt
        next_psi = method(psi, dt, time, params, independent, cutoff)

        # calculate difference between current and next psi
        en = difference(psi, next_psi)

        # run propogation while the difference is greater than acceptable error
        while en > epsilon
            # adjust time step
            dt *= epsilon / en
            # get the next MPS and calculate difference
            next_psi = method(psi, dt, time, params, independent, cutoff)
            en = difference(psi, next_psi)

        # incremement time and add it to the array
        time += dt
        times = vcat(times, [time])

        # accept wavefunction
        psi = deepcopy(next_psi)

        # calculate expectations
        energies = vcat(energies, [inner(psi, get_ham(time, params, independent), psi)])
        currents = vcat(currents, [inner(psi, get_current(time, params, independent), psi)])

        en1 = (en1 > 0) ? en1 : en

        # adjust for next time step
        # https://www.sciencedirect.com/science/article/pii/S0377042705001123
        ndt = dt * ((epsilon^2 * dt)/ (en * en1 * pdt))^(1/12)

        # update values for next iteration e_{n-1} -> e_n, dt_{n-1} = dt_n,
        # dt_n -> dt_{n+1}
        en1 = en
        pdt = dt
        dt = ndt
    return times, energies, currents
commented by (70.1k points)
Thanks for the question. I have two questions back:
1. how quickly does the energy change in time? I don't think there's any guarantee that the TEBD method will preserve the energy of the ground state due to various technical issues about the approximations it makes, but it should preserve it moderately well

2. what is the reason for the factor dt/2 in some of your gates but only dt without the division by 2 in others? I think the correct Trotter formula should always have the /2 here
commented by (130 points)
1. My evolution should run from 0 to about 241, but half the initial energy is lost within .5. My dynamic time step converges on about .005. I compared to the Heisenberg model example (https://itensor.github.io/ITensors.jl/stable/tutorials/MPSTimeEvolution.html) where the only change I made was to start evolution from the ground state as given by DMRG. This model loses about less than a hundredth of a percent of its initial energy over the entire evolution with a timestep of .1.

2. I am using the formula for TEBD2 given by https://www.sciencedirect.com/science/article/pii/S0003491619302532. I have also tried using the same gates used in the example referenced above which made no difference.
commented by (70.1k points)
Are you normalizing your state after each call to apply?
commented by (130 points)
I just tried that, it did not make a difference.
commented by (70.1k points)
Hi, so I would like to help with this, but to do so it would be really useful to have a fully working, minimal code that reproduces the issue. For example, from what I can see, the code above requires a ground state MPS as input, but I don't have the code you used to compute that ground state. Also the code calls a function phi_tl in it that isn't defined above.

Overall, it's hard to know a priori why you are finding that the energy isn't conserved. It could be a very subtle bug in the code somewhere, or perhaps a bug in an ITensor function. I think it's a case where the only way to find out is to run the code on the simplest case that reproduces the issue and then come up with various checks or test cases that draw out the bug.

Could you please provide the minimal working code? Thanks!

Please log in or register to answer this question.

Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.