DMRG fails for complex Hamiltonian (still Hermitian)

DMRG and Numerical Methods
, , , ,
1

So I have a very simple spin 1 complex Hamiltonian, with hopping terms but with an added complex phase for every hop to make it complex. I will write all my code out, so it’s long but you can copy everything and see it (doesn’t) work. First I initialize two complex hopping operators bp and bm

using ITensors,ITensorMPS, LinearAlgebra,Random
#Initialize complex local operators
dim = 3 #spin 1 model
theta = 0.5*pi #e^i*pi*theta to make it complex
dia = sqrt.(collect(1:dim-1)) / sqrt(2)
nopp = diagm(exp.(im .* (0:dim-1) .* theta))
nopm = diagm(exp.(im .* (0:dim-1) .* -theta))
bp_mat = diagm( 1=>dia )
bm_mat = diagm( -1=>dia)
bp_mat = bp_mat*nopp
bm_mat = nopm*bm_mat
display(bp_mat) #how "bp" look like in matrix
display(bm_mat) #how "bm" look like in matrix

ITensors.op(::OpName"bp",::SiteType"S=1") = bp_mat #bp
ITensors.op(::OpName"bm",::SiteType"S=1") = bm_mat #bm

Overall you can print bp_mat and bm_mat out to see it’s complex. Then I create a hopping Hamiltonian similar to Fermi-Hubbard

N=10 #number of sites
sites = siteinds("S=1",N;conserve_qns=true)
	
os = OpSum() #hopping operator like Hubbard, but with complex phase
for j=1:N-1
os += "bp",j,"bm",j+1
os += "bm",j,"bp",j+1
end
H = MPO(os,sites)

Then I try to solve with conserved quantum number

# Create the initial guess with "Up" and "Dn" 
density = 0.9 # density of "Up" to pick conserved sector
num_up = ceil(Int, density * N)
num_dn = N - num_up
psi_init = vcat(fill("Up", num_up), fill("Dn", num_dn))
psi_init = psi_init[Random.randperm(N)] #shuffle the list
psi0 = MPS(sites,psi_init)

Then I run DMRG

#Parameters of DMRG
nsweeps = 20
maxdim = [400]
cutoff = [1E-15]
weight = 1
no_state = 5
energies = [] #array of energy
states = [] #array of excited states
	
for energy_lvl in 1:no_state
    if energy_lvl == 1
        energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff,outputlevel=0,ishermitian=false)
		states = [psi] #insert ground state
		energies = [energy] #insert ground energy
    else
        energy, psi = dmrg(H, states, psi0; nsweeps, maxdim, cutoff,outputlevel=0,ishermitian=false,weight)
		push!(states, copy(psi)) #insert excited state
		push!(energies, copy(energy)) #insert excited energy
	end
end
display(energies)

However Julia DMRG is failing. As you can see if you run all this code, the energy as well as ground/excited states are all just the one for the real Hamiltonian case (when theta = 0), with some “noise” imaginary compenent that’s at order e-31, but at least it shows DMRG do recognize the complex Hamiltonian. Yet, we can vary theta and see that nothing changes, even though bp and bm obviously do change.

I have turned on ishermitian=false. I tried this first on SciPy eigsh (at N=10 it’s a small matrix so exact diagonalization also work) and it does show both the energy and the states change. What else can I do to improve this? Is there anything I did wrong? I’m new to DMRG, so I’ll be glad if you can help me with this. Thank you! :smiley:

2

The first thing I would suggest trying is making a random MPS as your initial state, and making it complex too.

You can do it like this:

linkdims = 10
psi_init = random_mps(ComplexF64, sites; linkdims)

where I just chose an arbitrary value for linkdims (the dimensions of each link, except the ones near the edges which have to be smaller).

3

Thank you for telling me a new way to generate initial guess! So with what you wrote I took the liberty to edit it as this to add in my code

psi0 = random_mps(ComplexF64, sites, psi_init; linkdims)

However it didn’t work. I actually tried modifying the initial guess before by adding im .* to it but it still didn’t work.

So I assume that it should work, but just doesn’t? I’m thinking maybe something’s wrong with ITensors.op or MPO, but I added ComplexF64 into them to and it still doesn’t work… What seems to be the problem?

4

I’m not sure what is wrong exactly. There may just be a bug in your code somewhere, such as in the definition of your Hamiltonian.

Are you sure that your Hamiltonian is Hermitian, the way you defined it in the code?

5

I think the error is that you also need to make your MPO complex.
I took your code, slightly modified it

using ITensors,ITensorMPS, LinearAlgebra,Random
#Initialize complex local operators
dim = 3 #spin 1 model
theta = 0.5*pi #e^i*pi*theta to make it complex
dia = sqrt.(collect(1:dim-1)) / sqrt(2)
nopp = diagm(exp.(im .* (0:dim-1) .* theta))
nopm = diagm(exp.(im .* (0:dim-1) .* -theta))
bp_mat = diagm( 1=>dia )
bm_mat = diagm( -1=>dia)
bp_mat = bp_mat*nopp
bm_mat = nopm*bm_mat
display(bp_mat) #how "bp" look like in matrix
display(bm_mat) #how "bm" look like in matrix

ITensors.op(::OpName"bp",::SiteType"S=1") = bp_mat #bp
ITensors.op(::OpName"bm",::SiteType"S=1") = bm_mat #bm

N=10 #number of sites
sites = siteinds("S=1",N;conserve_qns=true)
	
os = OpSum() #hopping operator like Hubbard, but with complex phase
for j=1:N-1
os += "bp",j,"bm",j+1
os += "bm",j,"bp",j+1
end
H = MPO(ComplexF64, os,sites)

# Create the initial guess with "Up" and "Dn" 
density = 0.9 # density of "Up" to pick conserved sector
num_up = ceil(Int, density * N)
num_dn = N - num_up
psi_init = vcat(fill("Up", num_up), fill("Dn", num_dn))
psi_init = psi_init[Random.randperm(N)] #shuffle the list
psi0 = MPS(ComplexF64, sites,psi_init)

#Parameters of DMRG
nsweeps = 20
maxdim = [400]
cutoff = [1E-15]
weight = 1
no_state = 5
energies = [] #array of energy
states = [] #array of excited states
	
for energy_lvl in 1:no_state
    if energy_lvl == 1
        energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff,outputlevel=0,ishermitian=false)
		states = [psi] #insert ground state
		energies = [energy] #insert ground energy
    else
        energy, psi = dmrg(H, states, psi0; nsweeps, maxdim, cutoff,outputlevel=0,ishermitian=false,weight)
		push!(states, copy(psi)) #insert excited state
		push!(energies, copy(energy)) #insert excited energy
	end
end
display(energies)

And I see different (though not dramatically different) energies for the states. To verify I also ran 10 states and see the expected trend (larger states have more positive energy)

energies = [ -1.9189859472289947 - 4.8172322733985737e-32im,
 -1.8007465064456778 + 1.111634010578569e-20im,
 -1.6825070656623584 + 3.4038184075122966e-32im,
 -1.6143537075597847 + 4.036843474568825e-26im,
 -1.4961142667764673 + 2.5661449275192192e-25im,
 -1.3749079866111464 - 1.5768148980800724e-27im,
 -1.3097214676743427 - 1.2414013868731206e-31im,
 -1.2566685460342446 - 1.6274348771315184e-30im,
 -1.1017966751289117 - 1.3050816443925843e-30im,
  -1.070286745302481 + 1.5395200310534437e-29im,
]