# Computing excited states using itensor (with DMRG)

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?

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;
}