Initial state calculation on 2D DMRG with Gaussian MPS

Hi there,

Recently I ran some DMRG calculations of a basic 2D Hubbard model with NN hopping and some onsite interaction, namely

H = \sum_{\langle i, j\rangle} c^{\uparrow}_i c^{\uparrow\dagger}_j + c^{\downarrow}_i c^{\downarrow\dagger}_j + h.c. +U\sum_i \hat{n_i }^{\uparrow \downarrow}

In practice, I’ve found that the DMRG algorithm, starting from a random guess of a product state with QN conservation, results in extremely slow convergences ( >300 sweeps and still not close to convergence). So I looked to the Gaussian MPS algorithm implemented in iTensor, and with my code


   # we assume a Lx * Ly, 2D square lattice, Ltotal = Lx * Ly
    # single particle H
    H = zeros((Ltotal, Ltotal))

    #hopping
    for i = 1:Ltotal
      
     # util function to get the nearest neighbor from certain order and geometry
      nns = get_nn(i, L)

      for nn in nns
        r = distance(i, nn)
        hop = hopping(decay, r)
        
        H[i, nn] = -t * hop
        H[nn, i] = -t * hop
      end 

    end 
    
   

    @show ishermitian(H)
    _, u = eigen(H)

  
    # Get the Slater determinant
    # Create an mps for the free fermion ground state

      # N = [Nup, Ndn] is the configuration
     
      @show N

      ϕup = u[:, 1:N[1]]
      ϕdn = u[:, 1:N[2]]

      ψ0 = slater_determinant_to_mps(sites, ϕup, ϕdn)


  return ψ0

Since up and down electrons have no interaction in the free electron case, I use the same H for both.
However, the initial state guess generated this way seem to be a lot worse than the random guess, and I’ve seen bizarre behavior of the DMRG solver, such as energy actually going up during sweeps.

I wonder if there’s something wrong with the way I generate these initial states, or the Gaussian MPS is only strictly applicable to 1D states.

Thanks in advance

The method to make Gaussian MPS does use a kind of 1D structure of the non-interacting problem you give it. Actually it must, since MPS itself is a 1D structure. So it’s true that Gaussian MPS is mainly applicable to quasi-1D systems, but I would not say it’s “strictly” applicable to 1D. It may work fine for quasi-1D i.e. small 2D systems if one tunes the parameters well.

The key parameters to try to adjust are these:

  • eigval_cutoff (default is 1E-8)
  • maxblocksize (default is the number of sites of the system)

Try setting the maxblocksize to Ly to start with, then see if a somewhat smaller or larger value of maxblocksize gives better results.

Also are you computing the energy of your Gaussian MPS after making it (using inner(psi',H,psi) where H is your MPO? You should see if it at least gives an ok energy before using it as an initial state. If the energy is very high, you could use this as a diagnostic to help you to tune the inputs to the slater_determinant_to_mps function to hopefully get a lower energy.

@miles Could you please comment more on what are the problems that once would incourr when trying to applying the GMPS algorithm to a 2D system? In the examples folder: ITensors.jl/ITensorGaussianMPS/examples at main · ITensor/ITensors.jl · GitHub there are examples with 2d systems.