# Non-Hermitian operator time evolution

edited

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