How to run DMRG over 3+ sites in C++ ?

Question

I'm attempting to adapt Sujay Kazi's work in
https://www.itensor.org/support/2212/modify-dmrg-optimize-over-different-number-sites-time-julia?show=2219#a2219
(On github: https://github.com/sujay-kazi/ITensors.jl/blob/master/src/mps/dmrg.jl)
using the dmrg.cc file in the tutorials
https://github.com/ITensor/ITensor/blob/v3/tutorial/06_DMRG/dmrg.cc ,
but I'm getting an error that the effective Hamiltonian and effective state have different dimensions.

I set up the Hamiltonian, called H using the autoMPO, then have an initial state guess psi0, and which to find the ground state and energy (psig, energyg). Then also set up the sweeps for the DMRG in the usual way. The problem I have is in the Davidson step where Heff and psig have different sizes. Heff is too small, but how do I grow Heff to encompass the length needed for the multi-site DMRG?

auto Heff = LocalMPO(H);

int dmrgStepSize = 3;
int finalLatticeSiteIndex = PhysicalLatticeLength;
float energyg; 
auto psig = psi0
    //Loop over sweeps
    for(auto sweepIteration : range1(sweeps.nsweep()))
        {
        //Loop over bonds
        for(int bondIndex = 1, subsweepIndex = 1; subsweepIndex <= dmrgStepSize; sweepnext(bondIndex,subsweepIndex,finalLatticeSiteIndex))
            {
            printfln("Sweep=%d, HS=%d, Bond=(%d,%d)",sweepIteration,subsweepIndex,bondIndex,bondIndex+dmrgStepSize);

            //Grow effective H
            psig.position(bondIndex);
            Heff.position(bondIndex,psig);

            // Form state for effective eigenvalue problem for this set of bonds 
            auto phig = psig(bondIndex);
            for(int stepsFromEdge : range1(dmrgStepSize-1) ){
                phig *= psig(bondIndex+stepsFromEdge) ; 
            }
            std::cout << " Makes it to davidson algorithm. \n" ;
            //Solve effective eigenvalue problem
            energyg = davidson(Heff,phig);

            //Update accuracy parameters
            //to pass to svd
            auto args = Args("Cutoff",sweeps.cutoff(sweepIteration),
                             "MaxDim",sweeps.maxdim(sweepIteration),
                             "MinDim",sweeps.mindim(sweepIteration));

            //Add code:
            //
            // SVD phi and store results 
            // 
            auto spec = psig.svdBond(bondIndex,phig,(subsweepIndex==1?Fromleft:Fromright), Heff, args) ; 

            } // for loop over bondIndex

        } // for loop over sweepIteration

Answer 1

Hi Jared,
The short answer to your question is that you will need to modify the LocalMPO code if you want to have it leave three sites exposed, or unprojected. Right now, that code does support a variable number of sites, but the only values supported are 0,1, and 2. One reason for this is that the most efficient pattern of contraction in the .product method of LocalMPO depends on the number of sites, so it can be important to optimize the code for each case. Another reason of course is that 3 site DMRG is just a less standard technique.

So the main steps would be to go through that code and find every place that the variable nc_ is used, and add cases for nc_==3. Then when the LocalMPO is constructed, make sure to set Heff.numCenter(3)

Best regards,
Miles

Comments

Hi Miles,

Thanks for the insight. I'll try to make those adjustments over the weekend.

Jared

---

Hi Miles,

For sure the things that have if( nc_ == 2) else if (nc_ == 1) etc in them need  to be extended, but I'm confused at the level of overhauling needed in the localmpo code. For instance,if I'm just inputting an MPO and a NumCenter argument, then I think this is the primary snippet of code to focus on for calling localMPO (line 275 of https://github.com/ITensor/ITensor/blob/v3/itensor/mps/localmpo.h)

inline LocalMPO::
LocalMPO(const MPO& H,
         const Args& args)
    : Op_(&H),
      PH_(H.length()+2),
      LHlim_(0),
      RHlim_(H.length()+1),
      nc_(2),
      Psi_(0)
    {
    if(args.defined("NumCenter"))
        numCenter(args.getInt("NumCenter"));
    }
Do RHlim and PH_ need to depend on the NumCenter argument? What I'm imagining is that we add fictitious sites past the actual physical system and set the bonds to those sites as 0 and so the +1 or +2 in the arguments of RHlim_ and PH_ account for the extra sites in the usual 2-site DMRG?
 

If I adapt something like
       if(nc_ == 2)
            lop_.update(Op_->A(b), Op_->A(b+1), L(), R());
        else if(nc_ == 1)
            lop_.update(Op_->A(b), L(), R());
        else if(nc_ == 0)
            lop_.update(L(),R());

to add an additional operator as M1 (middle 1)

        if(nc_ == 3)
            lop_.update(Op_->A(b), Op_->A(b+1), Op_->A(b+2), L(), M1(), R());
        else if(nc_ == 2)
            lop_.update(Op_->A(b), Op_->A(b+1), L(), R());
        else if(nc_ == 1)
            lop_.update(Op_->A(b), L(), R());
        else if(nc_ == 0)
            lop_.update(L(),R());
        end


then I get an error that M1 is undefined, which makes sense since L and R are defined in statements like

    ITensor const&
    L() const { return PH_[LHlim_]; }
    // Replace left edge tensor at current bond
    void
    L(ITensor const& nL) { PH_[LHlim_] = nL; }
    // Replace left edge tensor bordering site j
    // (so that nL includes sites < j)
    void
    L(int j, ITensor const& nL);

But for M1, would something like this work

    ITensor const&
    M1() const { return PH_[LHlim_ + 1]; }
    // Replace left edge tensor at current bond
    void
    M1(ITensor const& nM1) { PH_[LHlim_+1] = nM1; }
    // Replace left edge tensor bordering site j
    // (so that nM1 includes sites < j)
    void
    M1(int j, ITensor const& nM1);

?

PH_ is some vector of ITensors, but I'm not really sure what nL/nR/nM1 really means in this code. They're defined as ITensor objects (object used loosely here - I'm not sure if it's exactly a CC+ object in the jargon) and a constant, so I suppose it shouldn't be changed mid-call to L, R, or M1. Searching the localmpo.h file linked to above, I don't see any calls to L() or R() that don't use empty parentheses for me to try to figure it out.

I attempted (admittedly, not very hard) to look up tutorials on classes in C++, and they define things like a rectangle or other cute but not very useful classes that always take the same inputs.


As an aside, I'm not sure that this is the best path to my goal.
My main purpose of this is to find more accurate MPS approximations of a pair of nearly degenerate states.

My *hope* is that enlarging the number of sites that the DMRG optimizes over would help improve the ability to distinguish between the two nearly degenerate states. Using a randomMPS fails (and a random initial ansatz performs worse on exact diagonalization code in Fortran compared to other initial ansatz that may as well be random with respect to the ground state), and adding noise doesn't seem to help (I've tried variants of maximum number of bonds, noise size, etc).

- I haven't attempted using pinning fields at the opposite ends of the lattice to get the state close, then rerunning without the pinning fields to optimize.
- Another thing I haven't tried is taking a MPS, and (somehow) splitting it into two to try to find better approximations to the states then rerunning the DMRG using those as the starting states.

Jared

---

Hi Jared,
Indeed adapting the code to 3 sites is a good bit of work. Instead of extending the system with extra or fictitious sites, the better thing to do would be to change the sweeping logic so that the sweep stops 3 sites from the rightmost end instead of 2 sites or 1 site.

But you're right there are many other details that have to be changed. Most of them can be found, however, by looking for where the nc_ field is used in the code, and comparing the difference between the 1 and 2 site cases, then figuring out what the 3 site case would be of that same part of the code.

But I agree it may not be needed for your problem. The 3 site approach would mainly help only in cases where the local parts of the Hamiltonian are sort of "underspecified" and the physics is inherently longer ranged. An example would be if the Hamiltonian had *only* 3rd neighbor terms and no nearest-neighbor terms. Or if one was treating spinful fermions (electrons) by mapping to spinless fermion sites where odd sites represent spin-up electrons and even sites spin-down electrons. So again there the 2 site approach can get stuck because it's effectively a 1-site DMRG in disguise.

So you should consider whether your problem is really likely of that type, or more to do with inherently difficult physics. If the gap between your states really is tiny, then it may be better to just do a ton of sweeps, or come up with ways to prepare careful / physically motivated initial states or as you mentioned modify the Hamiltonian such as at the boundaries. Definitely study very small systems too at first. Finally, you may want to try excited-state DMRG where you find one ground state then put an "orthogonality penalty" for the second state you find, and maybe you'll get both states that way (this approach works well for the S=1 Heisenberg chain which is an example of an SPT phase with nearly-zero-energy edge states).

I think you know how to access this excited state DMRG code in ITensor, but if not please ask about it. We have a tutorial in the Julia docs about it -

Best,
Miles

P.S. thanks for being so active on the forum!

---

Hi Miles,

Thanks for the detailed reply. I think it's mostly a tough physics problem and multi-site DMRG is probably not the best time investment for this nearest-neighbor model.

My custom site set to keep interactions to nearest-neighbor didn't fix it automatically with the random MPS, but using a non-random trial state helped immensely. That non-random ansatz for everything except the phase transition/crossover regime, where the excited states can sometimes have a lower energy than the already calculated states. If it counts for anything, I think the documentation for excited states in the C++ code was easy enough to use (though, I don't recall at this point where I first learned it).

I was hoping to find a solution that would be a little more robust to the trial MPS to start the DMRG. I think because I only ever have 2-3 states crossing over one another, I should be able to make something that works and produces small residuals even in the phase transition regime (that is smeared out from being on such a small lattice).

For fermions in C++, I still have a tutorial in the works to demonstrate some of the features, but it's not priority #1 for me with other deadlines and work. Once I have things working sufficiently well, I will make that and share for folks here.

Happy to be on the forums even thought it's tough for me to give advice - I'm still learning the features and how to optimize for the models I'm actually interested in.

Cheers,
Jared

---

Hi Jared, yes sounds like a tough case. I'd recommend studying very small systems and working up from there. Also I'd recommend studying "nearby" Hamiltonians to see if you can get those under control (perturbing the edges, adding a field, strengthening certain terms over others, etc.) then if you can you can use the ground state of that Hamiltonian as an initial state for yours.