使用 itensor 计算激发态(使用 DMRG)

计算科学 矩阵 C++ 向量 张量
2021-11-30 15:37:32

我正在尝试计算一些哈密顿量的前几个激发态(我正在使用itensor它的DMRG算法)。为此,我正在尝试修改官方网页上给出的示例:

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

我不想只计算基态和第一个激发态,我还想计算第二个激发态。我的代码如下所示:

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

但我收到错误:

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)

在文档页面上:

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

我读过,std::vector wfs必须是0-indexed但我不明白,如果它必须包含多个先前计算的状态,而当前计算的状态必须与之正交,它怎么能是 0 索引?

1个回答

对于任何对此问题感兴趣的人,我找到了以下解决方案:

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