Avoiding local minima in DMRG

Question

Hi,

I am trying to use DMRG to compute a gap between Sz=0 and Sz=1 states.
My (fermionic) model is by no means intricate: a given 1D tight-binding model at half-filling with Hubbard interactions.
I am sure there is a gap in my particular system, but I am getting 0.
I tested several different things, so let me try to summarize my struggles:

1) Getting stuck in local minima (which depends on the initial state) while conserving symmetries

In principle, I want to converse N and Sz in this type of calculations: first, calculations are faster; second, I can get the gap by simply computing groundstates with the desired QNs. When I do that, I have to initialize my system in a state that has the QNs that I want. If I use Fock states, I can easily get stuck in local minima (as replied by Miles in http://itensor.org/support/105/sweeps-system-size-clarification-choice-intial-state-dmrg). I believe this is what happens in my code, even though I've tried many different initial guesses (all of them Fock states). For that reason, now I want to feed DMRG with an "arbitrary" combination of states (similarly to what is asked in http://www.itensor.org/support/1359/how-to-create-a-particular-initial-state-as-a-mps, http://itensor.org/support/1357/setting-random-initial-product-state, and http://itensor.org/support/490/starting-initial-state-as-product-state-of-singlets?show=490#q490) that still has the right N and Sz. However, I have not been able to create states as in the first and second links while conserving QNs (I get errors saying "Setting Tensor element non-zero would violate its symmetry"). In addition, I don't fully understand the syntax of the solution proposed in the third link to make something similar for electrons.
Let me give a concrete example. Let's say that instead of starting with a Neel-like Fock state, |up,dn,up,dn,...>, I want to start with something like 1/sqrt(2)(|up>+|dn>)1/sqrt(2)(|up>-|dn>)..., which has the same N and Sz. How can I run DMRG, starting from this initial state, and still conserving N and Sz?

2) "ConserveQNs"==false not working properly in version v3 for electrons?

While I could not fix the previous issue, other option that I had in mind was to not conserve symmetries at all (which allows me to start with a randomized wave function via the "randomMPS" function) and then look for the ground state (which has Sz=0) and to excited states (until I get one with Sz=1). Since the Sz=1 can be degenerate, to explore this option I would probably need to add a Zeeman splitting term in the Hamiltonian as well. However, it looks like there is a bug in version v3 when "ConserveQNs"==false is set for electrons (I think this is acknowledged in https://github.com/ITensor/ITensor/issues/311). For that reason, I am installing version v2 and exploring that route.

I believe this should be a quite common problem for people more experienced than me in DMRG calculations (I have just started). So, if someone has any other idea in mind to tackle this difficulty, please let me know.

Let me add that I have also played with the "noise" parameter in DMRG, but I couldn't fix my problems with that.

Best,
Gonçalo Catarina

EDIT: minor modifications to my questions, before any reply, for the sake of clarity.

Comments

Hi Goncalo, thanks. This is a lot of different questions though, so I think I'd have to ask some back to you to know how to really help. E.g. some of the middle part of your question (about making other initial states) may not be necessary if we can think of a simpler way to get your code un-stuck.

But briefly, before getting into many details:

1) what is your Hamiltonian and lattice, and how do you know your code is getting stuck?

2) if you want to do a calculation where you don't conserve particle number nor total Sz, what I'd recommend is that you just turn those off explicitly, like this:

auto sites = Electron(N,{"ConserveNf=",false,"ConserveSz=",false});

What that will do is leave fermion parity as a good quantum number, which will ensure that AutoMPO will still make a correct Hamiltonian (in the sense of having correct fermion signs in the operators).

---

Hi Miles,

Thanks a lot for the reply.

1) I know that my code is getting stuck for two reasons. First, I have other calculations that show a gap between Sz=0 and Sz=1 states. Second, I made a DMRG code to compute the groundstate energy of the SSH Hubbard model, starting from a random Fock state,

//function to generate MPS of random Fock state
auto randFock(auto sites, int Nsites)
{
    auto state = InitState(sites);
       
    for(int i=1; i <= Nsites; i++)
    {
        int randnr = rand()%4; //random integer between 0 and 3
       
        if(randnr == 0)
            state.set(i,"Up");
        else if(randnr == 1)
            state.set(i,"Dn");
        else if(randnr == 2)
            state.set(i,"UpDn");
        else
            state.set(i,"Emp");
    }
   
    auto stateMPS = MPS(state);
   
    return(stateMPS);
}

and conserving only parity,

auto sitesNfSzf = Electron(Nsites,{"ConserveNf",false,"ConserveSz",false});
//conserves parity given by the initial state

and I get different results in every run. I understand that, by conserving parity, if my initial state has an odd number of electrons and the groundstate has an even number of electrons, I will never reach it. But in some cases it actually happens to converge for other states with even number of electrons (more often for smaller systems). This way, I got convinced that I need to feed DMRG with initial states that are not Fock states, but maybe I'm wrong.

2) I've done that, as I said above. However, as far as I understand, parity is <N> mod 2, and should not(??) be related with the fermionic signs in the operators. In version v2, at least, it looks like it is working properly.
In addition, while conserving parity, I can't start with a randomized wave function. For that reason, if the above problem persists, this will likely not solve it.

---

Regarding my code, below I send a minimal working example for my particular problem.
The function DMRGconvE is something that I've made to run DMRG iteratively until I reach convergence in the groundstate energy. All the rest is pretty much standard.
I've also chosen some particular parameters for which calculations run fast and for which I have other type of calculations (configuration interaction method) that give me a gap.
Finally, let me add that I've tried different Sz=0 and Sz=1 initial states (all of them Fock states), and different values of Nuc, and so far I always got a 0 gap.

#include "itensor/all.h"

using namespace std;
using namespace itensor;

//function to generate <Sz>=0 Neel MPS at hf
auto Neelhf(auto sites, int Nsites)
{
    auto state = InitState(sites);
       
    for(int i=1; i <= Nsites; i++)
    {
        if(i%2 == 1)
            state.set(i,"Up");
        else
            state.set(i,"Dn");
    }
   
    auto stateMPS = MPS(state);
   
    return(stateMPS);
}

//function to generate a given <Sz>=1 MPS at hf
auto Sz1test1(auto sites, int Nsites)
{
    auto state = InitState(sites);
   
    for(int i=1; i <= Nsites; i++)
    {
        if(i%2 == 1)
            state.set(i,"Up");
        else
            state.set(i,"Dn");
    }
    state.set(Nsites/2,"Up");
    state.set(Nsites/2 + 1,"Up");
   
    auto stateMPS = MPS(state);
   
    return(stateMPS);
}

//function to generate MPO of TB Hubbard Hamiltonian
auto TBHub(double t1, double t2, double U, int Nuc, int Natuc, auto sites)
{
    auto ampoH = AutoMPO(sites);
   
    for(int i=1; i<Nuc; i++)
    {
        for(int j=1; j<Natuc; j++)
        {   
            ampoH += t1,"Cdagup",Natuc*i,"Cup",Natuc*i-j;
            ampoH += t1,"Cdagup",Natuc*i-j,"Cup",Natuc*i;
            ampoH += t1,"Cdagdn",Natuc*i,"Cdn",Natuc*i-j;
            ampoH += t1,"Cdagdn",Natuc*i-j,"Cdn",Natuc*i;
           
            ampoH += U,"Nupdn",Natuc*i-j;
        }

        ampoH += t2,"Cdagup",Natuc*i,"Cup",Natuc*i+Natuc;
        ampoH += t2,"Cdagup",Natuc*i+Natuc,"Cup",Natuc*i;
        ampoH += t2,"Cdagdn",Natuc*i,"Cdn",Natuc*i+Natuc;
        ampoH += t2,"Cdagdn",Natuc*i+Natuc,"Cdn",Natuc*i;
       
        ampoH += U,"Nupdn",Natuc*i;
    }
    for(int j=1; j<Natuc; j++)
    {   
        ampoH += t1,"Cdagup",Natuc*Nuc,"Cup",Natuc*Nuc-j;
        ampoH += t1,"Cdagup",Natuc*Nuc-j,"Cup",Natuc*Nuc;
        ampoH += t1,"Cdagdn",Natuc*Nuc,"Cdn",Natuc*Nuc-j;
        ampoH += t1,"Cdagdn",Natuc*Nuc-j,"Cdn",Natuc*Nuc;
           
        ampoH += U,"Nupdn",Natuc*Nuc-j;
    }
    ampoH += U,"Nupdn",Natuc*Nuc;
   
    auto H = toMPO(ampoH);
   
    return(H);
}

//routine to implement DMRG until convergence in energy
tuple<Real,MPS> DMRGconvE(auto H, auto psi0, auto sweepsdflt,
    auto sweepsextra, int maxdimi, int maxdimstep, double precE)
{
    //default run
    auto [Edflt,psidflt] = dmrg(H,psi0,sweepsdflt,"Quiet");
   
    //extra runs until convergence set by precE
    double E = 0;
    auto psi = psi0;
    int maxdim = maxdimi;
    while(true)
    {
        auto [Eextra,psiextra] = dmrg(H,psidflt,sweepsextra,"Quiet");
       
        maxdim += maxdimstep;
        sweepsextra.maxdim() = maxdim;
       
        if(fabs(Edflt-Eextra) < precE)
        {
            E = Eextra;
            psi = psiextra;
            break;
        }

        Edflt = Eextra;
        psidflt = psiextra;
    }
   
    return( tuple<Real,MPS>(E,psi) );
}

// main
int main()
{
    /* parameters */
    double t1 = 1; //all energies in units of t1
   
    double t2 = 0.1;
    double U = 1;
    int Nuc = 2;
    double precE = 0.001;
    /* */
   
    /* system */
    int Natuc = 4;
    int Nsites = Natuc*Nuc; //nr of sites
   
    auto sitesNtSzt = Electron(Nsites,{"ConserveNf",true,"ConserveSz",true});
    //conserves <N> and <Sz> given by the intial state (default)
    /* */
   
    /* initial state */
    auto psi0Neelhf = Neelhf(sitesNtSzt,Nsites); //<Sz>=0 Neel MPS at hf
   
    auto psi0Sz1test1 = Sz1test1(sitesNtSzt,Nsites); //a given <Sz>=1 MPS at hf
    /* */
   
    /* Hamiltonian */
    auto HNtSzt = TBHub(t1,t2,U,Nuc,Natuc,sitesNtSzt); //MPO of SSH Hubbard Hamiltonian (for NtSzt)
    /* */

    /* DMRG parameters */
    //default sweeps (corresponding runs should be fast)
    auto sweepsdflt = Sweeps(5);
    sweepsdflt.maxdim() = 20,50,100,100,200;
    sweepsdflt.cutoff() = 1E-5,1E-5,1E-5,1E-8,1E-8;
    sweepsdflt.noise() = 1E-5,1E-6,1E-8,1E-10,1E-12;
   
    //extra sweeps
    auto sweepsextra = Sweeps(1);
    sweepsextra.cutoff() = 1E-8; //fixed at every sweep
    int maxdim;
    int maxdimi = 300;
    int maxdimstep = 100;
    sweepsextra.maxdim() = maxdimi; //updated at every sweep
    /* */
   
    /* run DMRG */
    //(for NtSzt, Neelhf)
    auto [ENtSztNeelhf, psiNtSztNeelhf] =
        DMRGconvE(HNtSzt,psi0Neelhf,sweepsdflt,sweepsextra,maxdimi,maxdimstep,precE);
       
    //(for NtSzt, Sz1test1)
    auto [ENtSztSz1test1, psiNtSztSz1test1] =
        DMRGconvE(HNtSzt,psi0Sz1test1,sweepsdflt,sweepsextra,maxdimi,maxdimstep,precE);
    /* */
    
    /* outputs */
    //gs energy
    cout << endl;
    cout << "#gs energy (NtSzt, Neelhf) = " << ENtSztNeelhf << " [t1]" << endl;
    cout << "#gs energy (NtSzt, Sz1test1) = " << ENtSztSz1test1 << " [t1]" << endl;
    /* */

    return 0;
}

---

Hi, thanks for posting some code.

I'm still missing some important information I think: your code conserves particle number and spin, and uses a randomly generated initial state (i.e. a random total particle number and total Sz), so then naturally it will give different, random results, correct?

---

Hi,

The most relevant code that I want to fix (my working example, second comment), conserves N and Sz, and does **not** use a randomly generated initial state. It starts with some Sz=0 and Sz=1 states.

The comment previous to that one was just to state why I think that, in general, if I start with a Fock state, I can get stuck in local minima. To show that, I took the SSH Hubbard model, initialized it with a randomly generated Fock state and did **not** conserve either N and Sz.

So, to recap:
I am almost sure my code (working example, second comment) is not giving me the correct results.
I believe I am getting stuck in local minima because I initialize with Fock states that probably have a small overlap with the groundstate.
I tried many different Fock states, and I never managed to solve it.
I then used a simpler model (SSH Hubbard) to see that indeed it is quite easy to get stuck in local minima.
Given these checks, I now want to initialize DMRG with a more arbitrary guess. One option is to build MPS that are not Fock states but linear combinations (in the original post I give an example of something that I would like to know how to do). The other option, which is not convenient but that I am also exploring, is to not converse symmetries at all and look for groundstates and excited states.

---

Ok thanks for explaining. The code is fairly complicated, though, so if I start investigating it I might not study the same Hamiltonian as you. Could you please describe a particular Hamiltonian you are finding you can't get to converge? And for that Hamiltonian, could you say what is the result you expect to find? Thanks -

---

Hi,

First of all, thanks a lot for your patience.

The Hamiltonian is inside the function 'TBHub' in my minimal working example code.
It describes a tight-binding model with a unit cell of 4 atoms that has C3 symmetry.
There is intracell hopping (t1) between the middle atom with the other 3, and intercell hopping (t2) between middle atoms.
In addition to that, there is Hubbard repulsion (U).

I expect to get an energy gap between Sz=0 and Sz=1 states at half filling.

Let me try to clarify my code such that it doesn't look so complicated.
In this code, I run DMRG, conserving N and Sz, starting from i) an initial state with Sz=0 at half filling (created by the 'Neelhf' function) and ii) an initial state with Sz=1 at half filling (created by the 'Sz1test1' function).
As I said above, the 'TBHub 'function creates the Hamiltonian.
The 'DMRGconvE' function is just a routine that I have implemented to do more and more sweeps until I reach convergence in energy set by the variable 'precE'. I believe this is the only part of the code that is not so standard, but it is working properly. Basically, I do 5 sweeps with "typical" DMRG parameters, and then I keep doing routines of 1 sweep (using as initial guess the wave function calculated in the previous sweep and increasing maxdim) until I achieve convergence.
In the end, I print the groundstate energy that I obtain starting with i) the initial state with Sz=0 and ii) the initial state with Sz=1. I am getting the same value, which implies a 0 gap.

My guess is that the initial Sz=0 state that I set is leading me into a local minima that is degenerate with the Sz=1 excited state. I tried other initial Sz=0 states (always in the form of Fock states), and never managed to get a different result.