Excited State DMRG doesnt converge - Julia version

[Manas]

I have tried it and the excited state converges for certain values of parameter set in the Hamiltonian. The code snippets and output file seen is attached as jpeg images. However for certain other values of the parameter within the Hamiltonian the code is always returning the ground state. This is seen from the inner product between ground and excited state which returns ±1 as opposed to 0 . I have maintained large number of sweeps , added noise, increased bond-dimension steadily and even made the weight parameter very high (much higher than the gap) and changed initial guess MPS. All of these can be seen in the code. Do you have any idea why am I seeing this behavior ? What can be the possible solutions to fix it ?

Any help on the above problems would be extremely beneficial to me. Thanks in advance

Manas

[miles]

Could you please cut and past the code and outputs instead of posting images please? The images are hard to see. Thank you –

[Manas]

Hi Miles,
Sure, here are the code snippets

function Ham_constructor(sites, coeff, bnd_str, on_site =false, D_vec=nothing)

    i=1
    target = i+1

    # Input operator terms which define a Hamiltonian

    ampo = OpSum()

    while i <= (length(sites))

        if (bnd_str == "pbc" && i == length(sites))
           target = 1
        elseif (bnd_str == "obc" && i == length(sites))
           break
        end        
                    
        ampo += (0.5*coeff[1], "Sz * Sz", i) 
        ampo += (0.5*coeff[1], "Sy * Sy", i) 
        ampo += (0.5*coeff[1], "Sx * Sx", i)

        ampo += (coeff[2], "Sz", i, "Sz", target)

        ampo += (coeff[3], "Sx", i, "Sx", target)
        ampo += (coeff[3], "Sy", i, "Sy", target)


        ampo += (coeff[4], "Sz * Sz", i, "Sz * Sz", target)

        ampo += (coeff[5], "Sx * Sx", i, "Sx * Sx", target) 
        ampo += (coeff[5], "Sx * Sy", i, "Sx * Sy", target)
        ampo += (coeff[5], "Sy * Sx", i, "Sy * Sx", target)
        ampo += (coeff[5], "Sy * Sy", i, "Sy * Sy", target)

        ampo += (coeff[6], "Sz * Sx", i, "Sz * Sx", target)         
        ampo += (coeff[6], "Sz * Sy", i, "Sz * Sy", target)         
        ampo += (coeff[6], "Sx * Sz", i, "Sx * Sz", target)         
        ampo += (coeff[6], "Sy * Sz", i, "Sy * Sz", target)         

        ampo += (coeff[7], "Sx * Sx", i, "Sy * Sy", target)
        ampo += (-coeff[7], "Sx * Sy", i, "Sy * Sx", target)
        ampo += (-coeff[7], "Sy * Sx", i, "Sx * Sy", target)
        ampo += (coeff[7], "Sy * Sy", i, "Sx * Sx", target)

        if (on_site == true)
            println(D_vec[i])
            ampo += (D_vec[i], "Sz * Sz", i) 
            ampo += (-D_vec[i]*(1/3), "Sz * Sz", i) 
            ampo += (-D_vec[i]*(1/3), "Sy * Sy", i) 
            ampo += (-D_vec[i]*(1/3), "Sx * Sx", i)
        end

        i+=1
        target = i+1

    end
    
  # Convert these terms to an MPO tensor network
    H = MPO(ampo, sites)
    return H
end



function DMRG_driver(sites, H, state_reqd, init_state=nothing, gr_state=nothing, weights=nothing)

    let
        # Create an initial random matrix product state
        #Random.seed!(1234)
        #psi0 = randomMPS(sites, 60)

        state = [isodd(i) ? "Up" : "Dn" for i=1:length(sites)]
        #state = ["Z0" for i=1:length(sites)]
        psi0 = productMPS(sites,state)
        @show flux(psi0)

        # Plan to do 40 DMRG sweeps:
        sweeps = Sweeps(60)

        # Set maximum MPS bond dimensions for each sweep
        #setmaxdim!(sweeps, 10, 30, 50, 80, 100, 150, 200,300)
        setmaxdim!(sweeps, 10, 10, 10, 20, 20, 30, 30, 50, 80, 80, 100, 100, 150, 180, 200, 250, 300)

        # Set maximum truncation error allowed when adapting bond dimensions
        setcutoff!(sweeps, 1e-5)
        setnoise!(sweeps, 1e-5, 1e-5, 1e-7, 1e-7, 1e-9, 1e-9, 1e-10, 1e-11, 1e-12, 1e-12)
        @show sweeps


        # Run the DMRG algorithm, returning energy and optimized MPS
        if state_reqd == "gr"
             energy, psi = dmrg(H, psi0, sweeps)
        elseif state_reqd == "1st_exc"
             energy, psi = dmrg(H, [gr_state], init_state, sweeps; weights)

        end

        return energy, psi
    end

end


function main(spin_type, no_spins, state_reqd, bnd_str)

    # Create N spin-one degrees of freedom
    #sites = siteinds("S=$spin_type", no_spins, conserve_qns=true)
    sites = siteinds("S=$spin_type", no_spins)

    # Alternatively can make spin-half sites instead
    #sites = siteinds("S=1/2",no_spins)

    D_tuple=IterTools.product([-8,0,8], [-8,0,8])

    energy_arr = zeros(length(D_tuple))
    var_arr = zeros(length(D_tuple))
    ovlp_gr = zeros(length(D_tuple))


    for (index, D_tup) in enumerate(D_tuple)

        @show D_tup
        #print("one begins", index)
        D1 = D_tup[1]
        D2 = D_tup[2]
        D_vec = [D1 + ((-1)^i)*D2 for i in 0:length(sites)-1]   #local anisotropy as parameter for Hamiltonian
        #println("D array", D_vec)

        H = Ham_constructor(sites,[-5.58, 9.53, -8.97, 1.27, 6.59, -3.18, 5.04], bnd_str, true, D_vec)

        if state_reqd=="gr"        
            energy, psi = DMRG_driver(sites, H, state_reqd)

        elseif state_reqd == "1st_exc"
            energy_gr, psi_gr = DMRG_driver(sites, H, "gr")
            psi_ex_init = psi_gr
            weights = 300
            energy, psi = DMRG_driver(sites, H, state_reqd, psi_ex_init, psi_gr, weights)
            @show inner(psi, psi_gr)
        end
        

        # Compute energy variance
        #************************

        H_sq = inner(H,psi,H,psi)
        H_avg = inner(psi', H , psi)
        var = H_sq-H_avg^2
        @printf("Var of final state = %.12f\n", real(var))
        @printf("Avg energy in final state = %.12f\n", real(H_avg))
        #@printf("Per exciton energy in final state computed from <H> = %.12f\n", real(H_avg)/length(sites))
        @printf("Per exciton energy in final state computed directly = %.12f\n", real(energy)/length(sites))

        energy_arr[index] = real(H_avg)/length(sites)
        var_arr[index] = real(var)
        ovlp_gr[index] = real(inner(psi, psi_gr))

       **** REST OF THE CODE FOLLOWS ****

       end

end

The outputs are as follows:

For no_spins = 10, S=1 spin chains

Parameter set (D1 = 0, D2 =-8)

#**************************************

EXCITED STATE NOT CONVERGED; RETURNS GROUND STATE

.
.
.
After sweep 59 energy=-36.72260723778795 maxlinkdim=4 maxerr=3.20E-07 time=0.029
After sweep 60 energy=-36.722607237787955 maxlinkdim=4 maxerr=3.20E-07 time=0.029

inner(psi, psi_gr) = 0.9999999999906718 + 6.730727086790012e-16im #NOT CONVERGED EXCITED STATE AS OVERLAP WITH GR STATE ABSOLUTE

Var of final state = 0.000825051831
Avg energy in final state = -37.722607237769
Per exciton energy in final state computed directly = -3.672260723779

Parameter set (D1 = -8, D2 =-8)

#**************************************

EXCITED STATE CONVERGED

.
.
.
After sweep 59 energy=-56.3847361997004 maxlinkdim=128 maxerr=9.86E-12 time=2.008
After sweep 60 energy=-56.38473619970033 maxlinkdim=128 maxerr=9.86E-12 time=2.021

inner(psi, psi_gr) = 4.98411099630891e-6 - 1.21952253850921e-10im # WELL CONVERGED EXCITED STATE AS INNER WITH GROUND STATE IS NEARLY ZERO

Var of final state = 0.000000374540
Avg energy in final state = -56.384736199725
Per exciton energy in final state computed directly = -5.638473619970

[miles]

Thanks but also could you please explain, in a minimal way, what you were trying to do, what output you expected, and what output you got instead? The code is rather complicated for me to read and e.g. I’m not sure what D1 and D2 are etc.

[Manas]

Hi Miles,
Sure, the points are answered as below

1) What am I trying to do ?

A: I am trying to simulate the ground and 1st excited state of the following Hamiltonian

H_{act}= H + \sum_j D_j (Sj^z)^2

where H = \sum_j H_{bond}(S_j^z, S_{j+1}^z) with the latter defined as

c_i are coefficients defined as [-5.58, 9.53, -8.97, 1.27, 6.59, -3.18, 5.04] in serial order from c_0 to c_6. Also S_1^{x,y,z} are 3 \times 3 S=1 angular momentum operators.
Note that D_j = \{D_1 + (-1)^i D_2 \:\: \forall i \in N \} with N being the length of the chain. Thus (D_1, D_2) are parameters that define the anisotropies in the Hamiltonian matrix.

2) What output I am expecting and what did I get?

A: For each parameter pair (D_1, D_2) I was expecting to have the converged 1st excited state for following (D_1, D_2) as

(D_1, D_2) = [(-8, -8), (-8, 0), -(-8,8), (0,-8), (0, 0), (0, 8), (8, -8), (8, 0), (8, 8)]

While the ground state converges fine, the first excited state calculations returns the right result only for certain (D_1, D_2) (like in example output for D_1 = -8, D_2 = -8) but does’nt converge at all in other cases. I have personally verified through exact diagonalization that the first excited state is non-degenerate and is gapped especially for low magnitudes of D_1, D_2 . This is true even for small chain length like N = 4, 6, 10 etc where the gap is even higher. I have also obtained the full wavefunction of the 1st excited state from exact diagonalization in N=6 for (D_1=D_2=0) , however the above code fails to produce the excited state result. I was hoping to get the same for N=4,6,8,16 with DMRG under periodic boundary conditions. However the ground state is returned each time as can be seen from

@show inner(psi, psi_gr)

which yields ± 1 instead of a very low nearly zero number (orthogonality)

3) What did I try?

A: I have tried the above code with high number of sweeps (60), slowly increasing bond dimension steadily from 10 all the way to even 400. I have incorporated noise for each sweep from 1e-5 to 1e-12 and set a cutoff at 1e-5. I have changed the initial guess from the ground state wavefunction retrieved (which is correct) to a randomMPS, a product state etc but to no avail. Each time certain (D_1, D_2) excited states are correctly returned but others fail.

[miles]

I see - thanks for the clear explanation. The best thing you can do here is not to make a loop over many different parameter values but just do calculations for them one at a time. Then you can focus on a case where the result does not come out correctly and dig in to investigating it in more detail. For example, you might need a different weight for some Hamiltonian parameters versus other Hamiltonian parameters.

Two things that stands out to me are:

• the weight you are using seems very high and could be causing problems
• for the case that didn’t work, I see that the bond dimension of the excited state is stuck at 4, which seems unlikely (unless the ground state also has similarly low entanglement and maybe it’s just a parameter regime where the entanglement is low?)

I would try studying that case that failed and experimenting with different weight values. It’s good that you are already trying many sweeps and using noise.