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