0 votes
asked by (260 points)
edited by

Dear all,

I would like to apply TDVP to solve the ODE as $$\frac{d}{dt}\rho(t) = -i \Delta ( \hat{a} +\hat{a}^{\dagger} )\rho(t) + i \Delta \rho(t) ( 2\hat{a} +\hat{a}^{\dagger} ),$$ where @@\rho@@ is @@2\times2@@ density matrix. The reuslts are similar but not correct (Here in red line). I will be appreciated if someone help me to point out the problem.

#include "itensor/all.h"
#include "tdvp.h" 
#include "basisextension.h"
#include "itensor/util/print_macro.h"

using namespace itensor;

int main()
{
//  parameters
int nsite = 2;
int ndvr = 1;
auto delta = 1.0;
auto eps = 1.0;
auto dt = 0.05;
auto nsteps = 800;

auto sites = Boson(nsite,{"ConserveQNs",false,"MaxOcc=",ndvr});
auto ampo = AutoMPO(sites);

ampo +=  -1_i * delta,     "A",   1;
ampo +=  -1_i * delta,  "Adag",   1;
// asymmetric
ampo +=   2_i * delta,     "A",   2;
ampo +=   1_i * delta,  "Adag",   2;

auto H = toMPO(ampo);
// initialization
auto state = InitState(sites);
state.set(1,"Emp");
state.set(2,"Emp");

auto psi = MPS(state);

auto sweeps = Sweeps(1);

sweeps.maxdim() = 20;
sweeps.niter() = 20;

// TDVP propagator
std::ofstream of("rho1.dat");
for(int n = 1; n <= nsteps; ++n)
    {
    if(n < 200)
        {
        // Global subspace expansion
        std::vector<Real> epsilonK = {1E-10, 1E-8};
        addBasis(psi,H,epsilonK,{"Cutoff",1E-8,
                                  "Method","DensityMatrix",
                                  "KrylovOrd",3,
                                  "DoNormalize",true,
                                  "Quiet",true});
        }

    auto [energy,psi1] = tdvp(psi,H,dt,sweeps,{"DoNormalize",true,
                                    "Quiet",true,
                                    "NumCenter",1});

    // calculate the element corresponding to $\rho(0,0)$
    auto el = std::vector<int>{{1,1}};

    auto rho1 = ITensor(1.);
    for(auto j : range1(2))
        {
        rho1 *= (psi1(j)*setElt(sites(j)=el[j-1]));
        }

    auto t1 = dt * n;

    auto rtmp1 = elt(realPart(rho1));

    printf(of,"%d %d \n", t1, rtmp1);


    }

of.close();
  }

Inserting following cod into tdvp.h solverMingru Yang code to return $\psi$.

std::tuple<Real,MPS> inline
tdvp(MPS & psi,
     MPO const& H,
     Cplx t,
     const Sweeps& sweeps,
     const Args& args = Args::global())
    {
    LocalMPO PH(H,args);
    auto energy = TDVPWorker(psi,PH,t,sweeps,args);
    return std::tuple<Real,MPS>(energy,psi);
    }

By the way, we can apply the inner function "applyMPO" instead of TDVP to solve this equation, however, which gives the same results as TDVP in above.

    #include "itensor/all.h"
    #include "tdvp.h"
    #include "basisextension.h"
    #include "itensor/util/print_macro.h"

    using namespace itensor;

    int main()
    {
    int nsite = 2;
    int ndvr = 1;
    auto delta = 1.0;
    auto eps = 1.0;
    auto dt = 0.01;
    auto nsteps = 4000;

    auto sites = Boson(nsite,{"ConserveQNs",false,"MaxOcc=",ndvr});
    auto ampo = AutoMPO(sites);

    ampo +=  -1_i * delta,     "A",   1;
    ampo +=  -1_i * delta,  "Adag",   1;

    ampo +=   2_i * delta,     "A",   2;
    ampo +=   1_i * delta,  "Adag",   2;

    auto H = toMPO(ampo);
    auto expH = toExpH(ampo,dt);

// Set the initial state
//    auto state = InitState(sites);
//    state.set(1,"Emp");
//    state.set(2,"Emp");
//    auto psi = MPS(state);

    auto psi = MPS(sites);
    auto s1 = sites(1);
    auto s2 = sites(2);
    auto rho = ITensor(s1,s2);
    rho.set(s1(1),s2(1),1.0);
    ITensor D;
    psi.ref(1) = ITensor(s1);
    psi.ref(2) = ITensor(s2);
    svd(rho,psi.ref(1),D,psi.ref(2));
    psi.ref(2) *= D;

// new
//    PrintData(psi);
// time propagation

    auto args = Args("Method=","DensityMatrix","Cutoff=",1E-9,"MaxDim=",500);

    std::ofstream of("rho1.dat");
    for(int n = 1; n <= nsteps; ++n)
        {
        psi = applyMPO(expH,psi,args);
        psi.noPrime().normalize();


        auto el = std::vector<int>{{1,1}};

        auto rho1 = ITensor(1.);
        for(auto j : range1(2))
            {
            rho1 *= (psi(j)*setElt(sites(j)=el[j-1]));
            }

        auto t1 = dt * n;

        auto rtmp1 = elt(realPart(rho1));

        printf(of,"%d %d \n", t1, rtmp1);


        }

    of.close();
    }
commented by (260 points)
Thanks to Dr. Mingru Yang.  She pointed out my error that the MPS shouldn't be normalized.

1 Answer

0 votes
answered by (62.7k points)

(Please see discussion above.)

Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.

Categories

...