Computing excited states using itensor (with DMRG)

Question

I am trying to compute first few excited states of some Hamiltonian (I am using itensor and its DMRG algorithm). To do so, I am trying to modify the example given on the official webpage:

http://itensor.org/docs.cgi?vers=cppv3&page=formulas/excited_dmrg

Instead of computing only the ground state and the first excited one, I want to compute also the second excited one. My code looks like this:

#include "itensor/all.h"

using namespace itensor;

int 
main()
    {
    int N = 100;

    //
    // Initialize the site degrees of freedom.
    //
    auto sites = SpinHalf(N,{"ConserveQNs=",false}); //make a chain of N spin 1/2's

    //Transverse field
    Real h = 4.0;

    //
    // Use the AutoMPO feature to create the 
    // transverse field Ising model
    //
    // Factors of 4 and 2 are to rescale
    // spin operators into Pauli matrices
    //
    auto ampo = AutoMPO(sites);
    for(int j = 1; j < N; ++j)
        {
        ampo += -4,"Sz",j,"Sz",j+1;
        }
    for(int j = 1; j <= N; ++j)
        {
        ampo += -2*h,"Sx",j;
        }
    auto H = toMPO(ampo);

    //
    // Set the parameters controlling the accuracy of the DMRG
    // calculation for each DMRG sweep. 
    //
    auto sweeps = Sweeps(30);
    sweeps.maxdim() = 10,20,100,100,200;
    sweeps.cutoff() = 1E-10;
    sweeps.niter() = 2;
    sweeps.noise() = 1E-7,1E-8,0.0;
    println(sweeps);

    //
    // Begin the DMRG calculation
    // for the ground state
    //
    auto [en0,psi0] = dmrg(H,randomMPS(sites),sweeps,{"Quiet=",true});

    println("\n----------------------\n");

    //
    // Make a vector of previous wavefunctions;
    // code will penalize future wavefunctions
    // for having any overlap with these
    //
    auto wfs = std::vector<MPS>(2); // <-- First change
    wfs.at(0) = psi0;

    //
    // Here the Weight option sets the energy penalty for
    // psi1 having any overlap with psi0
    //
    auto [en1,psi1] = dmrg(H,wfs,randomMPS(sites),sweeps,{"Quiet=",true,"Weight=",20.0});
    wfs.at(1) = psi1;  // <-- Second change

    auto [en2,psi2] = dmrg(H,wfs,randomMPS(sites),sweeps,{"Quiet=",true,"Weight=",20.0}); // <-- Third change

    //
    // Print the final energies reported by DMRG
    //
    printfln("\nGround State Energy = %.10f",en0);
    printfln("\nFirst excited State Energy = %.10f",en1);
    printfln("\nSecond excited State Energy = %.10f",en2); // <-- Fourth change

    //
    // The expected gap of the transverse field Ising
    // model is given by Eg = 2*|h-1|
    //
    // (The DMRG gap will have finite-size corrections.)
    //
    printfln("\nDMRG energy gap = %.10f",en1-en0);
    printfln("\nTheoretical gap = %.10f",2*std::fabs(h-1));

    //
    // The overlap <psi0|psi1> should be very close to zero
    //
    printfln("\nOverlap <psi0|psi1> = %.2E",inner(psi0,psi1));

    return 0;
    }

But I am getting error:

terminate called after throwing an instance of 'std::out_of_range'
  what():  vector::_M_range_check: __n (which is 2) >= this->size() (which is 1)

On the documentation page:

http://itensor.org/docs.cgi?page=classes/dmrg

I've read, that the std::vector wfs has to be 0-indexed. But I don't understand, how can it be 0-indexed, if it has to include more than one previously computed states, to which the currently computed one has to be orthogonal?

Answer 1

For anyone interested in this problem, I found the following solution:

#include "itensor/all.h"

using namespace itensor;

int 
main()
    {
    int N = 100;

    //
    // Initialize the site degrees of freedom.
    //
    auto sites = SpinHalf(N,{"ConserveQNs=",false}); //make a chain of N spin 1/2's

    //Transverse field
    Real h = 4.0;

    //
    // Use the AutoMPO feature to create the 
    // transverse field Ising model
    //
    // Factors of 4 and 2 are to rescale
    // spin operators into Pauli matrices
    //
    auto ampo = AutoMPO(sites);
    for(int j = 1; j < N; ++j)
        {
        ampo += -4,"Sz",j,"Sz",j+1;
        }
    for(int j = 1; j <= N; ++j)
        {
        ampo += -2*h,"Sx",j;
        }
    auto H = toMPO(ampo);

    //
    // Set the parameters controlling the accuracy of the DMRG
    // calculation for each DMRG sweep. 
    //
    auto sweeps = Sweeps(30);
    sweeps.maxdim() = 10,20,100,100,200;
    sweeps.cutoff() = 1E-10;
    sweeps.niter() = 2;
    sweeps.noise() = 1E-7,1E-8,0.0;
    println(sweeps);

    //
    // Begin the DMRG calculation
    // for the ground state
    //
    auto [en0,psi0] = dmrg(H,randomMPS(sites),sweeps,{"Quiet=",true});

    println("\n----------------------\n");

    //
    // Make a vector of previous wavefunctions;
    // code will penalize future wavefunctions
    // for having any overlap with these
    //
    auto wfs = std::vector<MPS>(1);
    wfs.push_back(psi0);

    //
    // Here the Weight option sets the energy penalty for
    // psi1 having any overlap with psi0
    //
    auto [en1,psi1] = dmrg(H,wfs,randomMPS(sites),sweeps,{"Quiet=",true,"Weight=",20.0});
    wfs.push_back(psi1);

    auto [en2,psi2] = dmrg(H,wfs,randomMPS(sites),sweeps,{"Quiet=",true,"Weight=",20.0});

    //
    // Print the final energies reported by DMRG
    //
    printfln("\nGround State Energy = %.10f",en0);
    printfln("\nFirst excited State Energy = %.10f",en1);
    printfln("\nSecond excited State Energy = %.10f",en2);

    return 0;
}