Hello,
It seems that the DMRG in Julia gives inconsistent result for weak interactions. The sample code is below.
ITensors.space(::SiteType"HardCore") = 2
ITensors.state(::StateName"0", ::SiteType"HardCore") = [1.0, 0.0]
ITensors.state(::StateName"1", ::SiteType"HardCore") = [0.0, 1.0]
function ITensors.op!(Op::ITensor,
::OpName"N",
::SiteType"HardCore",
s::Index)
Op[s'=>2,s=>2] = 1
end
function ITensors.op!(Op::ITensor,
::OpName"Adag",
::SiteType"HardCore",
s::Index)
Op[s'=>1,s=>2] = 1
end
function ITensors.op!(Op::ITensor,
::OpName"A",
::SiteType"HardCore",
s::Index)
Op[s'=>2,s=>1] = 1
end
t = 1.0
V1 = 1.0e-3
V2 = 2e-5
N = 20
sites = siteinds("HardCore",N)
ampo = OpSum()
# NN
for j = 1:N-1
ampo += -t,"Adag",j,"A",j+1
ampo += -t,"A",j,"Adag",j+1
ampo += V1,"N",j,"N",j+1
end
# NNN
for j = 1:N-2
ampo += V2,"N",j,"N",j+2
end
H = MPO(ampo,sites)
psi0 = productMPS(sites, n -> isodd(n) ? "0" : "1")
println(inner(psi0,H,psi0))
# Create an initial random matrix product state
psi0 = randomMPS(sites)
sweeps = Sweeps(100)
setmaxdim!(sweeps, 20,20,40,40,60,60,80)
setcutoff!(sweeps, 1E-11)
#noise!(sweeps, 1E-4,1E-6,1E-10,0)
N_observer = DMRGObserver(["N"],sites,energy_tol=1E-12)
energy, psi = dmrg(H,psi0, sweeps, observer=N_observer)
println("Final energy = $energy")
The result from the code is -12.378711861320. But exact diagonalization gives -12.378713430075452. There is a 1.6E-6 discrepancy. Although it is small, the DMRG in C++ version gives -12.378713429701, only 3.7e-10 discrepancy with ED. The sample code in C++ is below.
#include "itensor/all.h"
using namespace std;
using namespace itensor;
int
main(int argc, char* argv[])
{
auto Ns = 20;
auto t = 1.0;
auto V1 = 1E-3;
auto V2 = 2E-5;
auto sites = Boson(Ns,{"ConserveQNs=",true});
auto ampo = AutoMPO(sites);
for(auto j : range1(Ns-1))
{
ampo += -t,"Adag",j,"A",j+1;
ampo += -t,"A",j,"Adag",j+1;
ampo += V1,"N",j,"N",j+1;
}
for(auto j : range1(Ns-2))
{
ampo += V2,"N",j,"N",j+2;
}
auto H = toMPO(ampo);
//Perform 20 sweeps of DMRG
auto sweeps = Sweeps(20);
//Specify max number of states kept each sweep
sweeps.cutoff() = 1E-11;
sweeps.maxdim() = 20, 20, 40, 40, 60;
auto state = InitState(sites);
for(auto j : range1(Ns))
{
state.set(j, j%2 == 1 ? "Emp" : "1");
}
auto psi0 = MPS(state);
auto [energy0, psi0f] = dmrg(H, psi0, sweeps, {"Quiet=", true, "EnergyErrgoal=",1E-12,"EntropyErrgoal=",1E-12});
cout << energy0 << endl;
return 0;
}
What causes this difference in Julia? Thanks.
Jin