0
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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;
}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.