First Excited State Not Converging For Larger Systems

ITensor Julia Questions
1

I am trying to compute the first excited state of the Ising model in Julia and the DMRG converges without a problem for the ground state for any system size N but fails to do so for N > 10 for the first excited state. It does converge for the first excited state for up to N = 10. I have not tested N between 10 and 20 but then for N = 20 for example the ground state converges as expected and then the first excited state energy blows up immediately and is much more negative than the ground state.

I am having the same problem for the Schwinger model and there the situation is even worse, maybe because the Hamiltonian is more complicated, but there the first excited state algorithm fails to converge for N > 2 and in fact there is such a big numerical instability that I get errors of the form: ERROR: LoadError: operator does not appear to be hermitian: 17.330705523490913 vs 2.772575076898464e9.

I am checking my ansatz for the first excited state and has very small overlap with the ground state and the same is true for the result of the first excited state with the ground state even if the energy of the first excited state blew up beyond the energy of the ground state.

I give below a full working code for the Ising model which ends up giving a ground state energy of -4.471756355671585 and a first excited energy of -81856.75774379515. I have tried increasing the bond dimension, or adding some noise or increasing the maximum number of sweeps but nothing worked.

Is there anything you can suggest?

using ITensors

function get_Ising_OpSum(N, J, g_z, g_x)
        
    ampo = OpSum()

    for n=1:N-1
        ampo += 4*J,"Sz",n,"Sz",n+1
    end

    for n=1:N
        ampo += 2*g_x,"Sx",n
        ampo += 2*g_z,"Sz",n
    end

    return ampo

end

function get_MPO_from_OpSum(OpSum, sites)

    return MPO(OpSum, sites)

end

N = 20
J = 0.0001
g_z = 0.1
g_x = 0.2
ns = 500
D = 20

sites = siteinds("S=1/2", N)
opsum_ising = get_Ising_OpSum(N, J, g_z, g_x)
H = get_MPO_from_OpSum(opsum_ising, sites)
initial_ansatz_0 = randomMPS(sites, D)
sweeps = Sweeps(ns, maxdim = D)

energy_0, psi_0 = dmrg(H, initial_ansatz_0, sweeps, ishermitian = true, maxdim = D)

# Ms = [psi_0]
# w = energy_0
# initial_ansatz_1 = randomMPS(sites, D)

# energy_1, psi_1 = dmrg(H, Ms, initial_ansatz_1, sweeps, weight = w, ishermitian = true, maxdim = D)

Heff = H + energy_0.*outer(psi_0', psi_0)
initial_ansatz_1 = randomMPS(sites, D)

println("Overlap of ansatz with gs: ", inner(psi_0, initial_ansatz_1))

initial_noise = 1e-2

noise_vector = LinRange(initial_noise, 0.0, ns)

energy_1, psi_1 = dmrg(Heff, initial_ansatz_1, sweeps, ishermitian = true, maxdim = D, noise = noise_vector)

println("Overlap of 1st excited state with gs: ", inner(psi_0, psi_1))

println(energy_0)
println(energy_1)
2

Hi, so the part of what you’re doing that looks a little suspicious to me is where you numerically add outer(psi_0',psi_0) to H. That step is done numerically (at a very high cost) and could make your resulting Hamiltonian not exactly Hermitian by a small amount, which might explain the unphysically large energy you got from DMRG.

Relatedly, you should instead use the “excited state targeting” mode of our dmrg function to do this kind of thing, as it uses what is conceptually the exact same approach of penalizing the overlap with psi_0 but implements it in a much more efficient and numerically stable way. (Not by adding the outer product of psi_0 to H at the beginning but instead by overlapping one copy of psi_0 with the current state inside of the DMRG algorithm and adding the resulting product with psi_0' to the “effective Hamiltonian” used within DMRG. For more details on this you can see this paper.)

3

Hi Miles,

Thanks a lot for the quick reply.

I have also tried

Ms = [psi_0]
w = energy_0
initial_ansatz_1 = randomMPS(sites, D)

energy_1, psi_1 = dmrg(H, Ms, initial_ansatz_1, sweeps, weight = w, ishermitian = true, maxdim = D)

which I believe is what you suggested above but it did not solve the issue. I also tried a different weight like w = 1.0 but the problem persists. These lines are the commented lines in my full working code above. I have written my own code which implements the one site update algorithm and I test it with N that will make iT dmrg E_1 diverge and my code is able to find the correct 1st excited state as compared to exact diagonalization. So from this test I don’t believe the physical parameters are at fault here and my code implements Heff = H - E_0*|psi_0><psi_0|, I am doing something wrong with iT dmrg… The reason I cannot use my code is that it is much slower than iT.

4

I see, so the issue is that you have used a negative value for the weight. It didn’t occur to us that anyone would do that so we’ll probably fix the code to throw an error in the future if someone passes a negative-valued weight parameter.

As you may know, putting a negative weight value would actually reward rather than prohibit a non-zero overlap of the excited state with the ground state.

Please try again with the excited-state-dmrg method but taking the absolute value of the weight or ensuring it is positive. Does it work better now?

1 Like
5

It is working now thanks a lot for the help! So a priori I did not know for a given set of parameters that the gs energy can be negative so in general it is good to take the absolute value. Maybe the library can do that instead of throwing an exception and just print a warning that it received a negative weight and has taken the absolute value? Also when I tried weight = 1.0 it did not work because of course one needs to push the ground state level above the 1st excited so that the 1st excited becomes the ground state for this effective Hamiltonian and the gap was greater than the weight. This is another a priori unknown so the larger the weight the better I guess just to be on the safe side. For example I give as weight the abs(E_0). If E_0 > 0 then its not guaranteed in general that 2E_0 > E_1 so it is again somewhat of a risk.

6

Yes, based on your post I added an error message to the code so it will notify anyone in the future who runs into this situation.

About the magnitude of the weight, I can’t think of a good automatic way for the code to determine it, so I think we will just have to continue to explain in the documentation and in the forum here that the user must choose it to be at least as large as the estimate size of the gap.