Difficulty with simulating a specific model with VUMPS (that works perfectly fine with regular DMRG)

ITensor Julia Questions
1

Hi,

Hope this is the right forum for this. I am currently using the ITensorInfiniteMPS implementation of VUMPS (GitHub - ITensor/ITensorInfiniteMPS.jl: A package for working with infinite matrix product states (MPS) with ITensor.) and I am having trouble simulating a particular model with VUMPS. The model in question is

\hat{H} = \sum_{i}d_i^\dag c_{i+1} d_{i+2}c^\dag_{i+3} + h.c. \sim \sum_{i}S_{i+1}^+S^-_{i+2} S^-_{i+3}S^+_{i+4} + h.c.

where c, d are two species of fermion per cell and where I have also given an equivalent (via Jordan-Wigner) spin/hardcore boson model. I’ve elected to implement the model using hardcore boson/spin-1/2 sites. My implementation is coded below:

function unit_cell_terms(::Model"DipoleHopping"; t)
  opsum = OpSum()
  
  dipoleHopNN     = t,"S+",2,"S-",3,"S-",4,"S+",5
  dipoleHopNN_HC  = t,"S-",2,"S+",3,"S+",4,"S-",5
  
  opsum += dipoleHopNN
  opsum += dipoleHopNN_HC

  return opsum
end

Though it seems like this model has only one fermion/spin-1/2 per unit cell, I have other terms in the model that enforce the two spin-1/2 in the primitive unit cell. I’ve neglected to include these terms for ease of presentation and because I’d like to understand this simpler model first.

I’m able to simulate this model perfectly fine using regular open-chain DMRG with a bond dimension of 50 (for N=52 spin 1/2’s and a cutoff of 1e-10). (I can also benchmark with an exact solution of the model, as it maps onto an XX chain with half the DOF when projected to a special conserved subspace.) However, when I simulate with VUMPS (a rather naïve implementation of it, as I’ve just copied the Hubbard model example and made only minimal modifications), I get the following ‘errors’:

The two-site subspace expansion produced a zero-norm expansion at (1, 2). This is likely due to the long-range nature of the QN conserving Hamiltonian.

(and similarly for other pairs of sites), and the algorithm only ever uses bond dimensions from 1-4. The results from this run are obviously nonsense (e.g. single cut entropies give me alternating log(2) and 0). Below, I’m posting a minimal section of my code that contains the VUMPS parameters.

using ITensors, ITensorMPS
using ITensorInfiniteMPS

base_path = joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src")
src_files = ["vumps_subspace_expansion.jl", "entropy.jl"]
for f in src_files
  include(joinpath(base_path, f))
end


##############################################################################

# VUMPS parameters

#

maxdim = 100 # Maximum bond dimension

cutoff = 1e-10 # Singular value cutoff when increasing the bond dimension

max_vumps_iters = 200 # Maximum number of iterations of the VUMPS algorithm at each bond dimension

vumps_tol = 1e-5

outer_iters = 5 # Number of times to increase the bond dimension

localham_type = MPO # or ITensor

conserve_qns = true

eager = true

model_params = (t=1)

##############################################################################

N = 8 # Unit cell size

initstate(n) = isodd(floor(Int,n/2)+1) ? "↑" : "↓"

s = infsiteinds("S=1/2", N; initstate, conserve_qns)

ψ = InfMPS(s, initstate)

model = Model"DipoleHopping"()

@show model, model_params

# Form the Hamiltonian

H = InfiniteSum{localham_type}(model, s; model_params...)

outputlevel = 1

vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters, outputlevel, eager)

subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)

# For now, to increase the bond dimension you must alternate

# between steps of VUMPS and subspace expansion (which outputs

# a new state that is equal to the original state but with

# a larger bond dimension)

println("\nRun VUMPS on initial product state, unit cell size $N")

ψ = vumps_subspace_expansion(H, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs)

S_infinite = [entropy(ψ, b) for b in 1:N]
@show S_infinite

nothing

I understand that my usage of VUMPS here could be very naïve, hence why I’m running into this issue. I would appreciate any suggestion on how I can make the results of VUMPS line up with the finite DMRG.

N.B.: To work with my custom model, I made a DipoleHopping.jl file with that function and put include("models/DipoleHopping.jl") in the ITensorInfiniteMPS.jl file. This is a rather complicated step that can make reproducing my code quite unwieldy, but it’s the only thing that’s worked for me.

Thanks very much, and please let me know if additional information is needed!

Best,

Amogh

2

If you include the package name, i.e. function ITensorInfiniteMPS.unit_cell_terms(::Model"DipoleHopping"; t) then you shouldnt need do that.

One easy thing to first check is that your Hamiltonian matches the open-chain DMRG Hamiltonian, by using your unit cell version to create a finite one (see the examples ).
Ruling out that there isn’t an issue or a trivial product state lurking, it also could be that the initial state isn’t good and you’ll need to provide something different (others have found success initializing from a small infinite DMRG calculation, or initializing from another model for example).

Lastly, this would take some work, but you can also try to get working the new global subspace expansion from ITensorTDVP.jl (see the news here: GitHub - ITensor/ITensorTDVP.jl: Time dependent variational principle (TDVP) of MPS based on ITensors.jl. )

3

Dear Ryan,

Thank you for the quick reply! Indeed my unit cell hamiltonian creates the correct MPO for open-chain DMRG (cross-checked with finite chain MPO that I put in “by hand”).

Definitely it would be good to try a more complicated ansatz than these product states. Is there any way I can import MPS’s from open chain DMRG computations, or infinite MPS’s from prior VUMPS calculations (e.g. the output of something easy like the Fermi-Hubbard model), for use as VUMPS ansätze? I’ve been saving/loading MPS’s using HDF5 in open-chain DMRG but I’m not sure how this functionality transfers over to VUMPS/InfiniteMPS.

4

There is limited functionality here unfortunately (there is a bit of a question of how to take a finite MPS and make it into a uniform state). Matt has added something to attempt this

However last time I tried this I don’t believe it had quantum numbers working (PRs welcome!)

The same save/load functionality exists here as well. Here’s an example on saving loading

using HDF5
using ITensorInfiniteMPS
using ITensorMPS

# do stuff to get a state
psi = ...
h5open("my_wavefunction.h5","w") do fo
  write(fo,"psi",psi)
  # example metadata, needs vector for hdf5 conversion
  write(fo, "model_params", [model_params]) 
end

#######

psi, model_params = h5open("my_wavefunction.h5","r") do fi
  read(fi,"psi",InfiniteCanonicalMPS), read(fi,"model_params")[1]
end
s = siteinds(psi) 

model = #...
# load H using matching indices from old calculation
H = InfiniteSum{localham_type}(model, s; model_params...)
5

The subspace expansion used in the VUMPS code in ITensorInfiniteMPS.jl is a 2-site expansion, and may just fail for certain Hamiltonians (there is a nice explanation of that issue in [2005.06104] Time Dependent Variational Principle with Ancillary Krylov Subspace). This isn’t a limitation of the VUMPS algorithm per se, it is just one subspace expansion approach (specifically, it is the algorithm outlined in Appendix B of [1701.07035] Variational optimization algorithms for uniform matrix product states), is relatively cheap, and which works for many cases.

Something else to try, if the suggestions by @ryanlevy are hard to get working, would be to implement infinite TEBD and do imaginary time evolution (see [quant-ph/0310089] Efficient simulation of one-dimensional quantum many-body systems, [0711.3960] The iTEBD algorithm beyond unitary evolution) for some number of steps before running VUMPS, which will push the system towards the ground state and also entangle it and hopefully circumvent the subspace expansion issues you are seeing, which are likely caused by starting out on a product state. Another option for starting out with a state beyond a product state is to take a product state in the symmetry sector you are interested in and apply a random unitary circuit with infinite TEBD, that kind of state works well as a starting state for finite DMRG so likely would work well for VUMPS as well. Otherwise, you can use infinite DMRG ([0804.2509] Infinite size density matrix renormalization group, revisited) at the beginning of the calculation, which is related to VUMPS but is closer to finite DMRG in that it starts by running DMRG on a small finite system and grows it towards the thermodynamic limit, so if finite DMRG works it may work better for starting out the calculation.

Unfortunately we don’t have infinite TEBD or infinite DMRG implemented in ITensorInfiniteMPS.jl so you will have to implement those yourself.

6

Thanks very much to Ryan and Matt for the help! Trying a more intelligently constructed ansatz avoided the issue in my original question. It’s likely that simply not having a trivial product state as an ansatz fixed the issue.