Hi miles
I find one example.The hamiltonian is H=He+Hph+He-ph
$$
H_{e}=-t\sum_{i,\sigma}(c^{\dag}_{i\sigma}+h.c.)+U\sum_{i}n_{i\uparrow}n_{i\downarrow} \\
Hph=\Omega b^{dag}b \\
He-ph=-\frac{g}{\sqrt(N)}(b^{\dag}+b)\sum_{i,\sigma}(-1)^{i}n_{i\sigma}
$$
And I write a code. When I execute the code, the energy of every iteration is not always descending. It will be fluctuating, as -10,-11,-10,-9.... I try different parameters. However the result is same.
#include "itensor/all.h"
using namespace itensor;
using Holstein = MixedSiteSet<ElectronSite,BosonSite>;
int
main()
{
auto N = 20;
auto U = 0.5;
auto J = 1.0;
auto g = 1.0;
auto sigma = 1.0;
auto sites = Holstein(N,{"ConserveNf=",true,
"ConserveNb=",false,
"MaxOcc=",50});
auto ampo = AutoMPO(sites);
for(int j=1;j <= N-2; j+=2)
{
ampo += -J,"Cdagup",j+2,"Cup",j;
ampo += -J,"Cdagup",j,"Cup",j+2;
ampo += -J,"Cdagdn",j+2,"Cdn",j;
ampo += -J,"Cdagdn",j,"Cdn",j+2;
}
ampo += -J,"Cdagup",1,"Cup",19;
ampo += -J,"Cdagup",19,"Cup",1;
ampo += -J,"Cdagdn",1,"Cdn",19;
ampo += -J,"Cdagdn",19,"Cdn",1;
for(int j=1;j < N; j += 2)
{
ampo += U,"Nup",j,"Ndn",j;
}
for(int j=1;j<N;j +=4)
{
ampo += g,"Adag",2,"Nup",j;
ampo += g,"Adag",2,"Ndn",j;
}
for(int j=3;j<N;j +=4)
{
ampo += -g,"Adag",2,"Nup",j;
ampo += -g,"Adag",2,"Ndn",j;
}
ampo += sigma,"Adag",2,"A",2;
auto H = toMPO(ampo);
auto sweeps = Sweeps(20);
//very important to use noise for this model
sweeps.noise() = 1E-6,1E-6,1E-8,1E-12;
sweeps.maxdim() = 10,20,100,100,200;
sweeps.cutoff() = 1E-10;
auto state = InitState(sites);
for(auto j : range1(N))
{
if(j%4==1) state.set(j,"Up");
else if(j%4==3) state.set(j,"Dn");
else state.set(j,"0");
}
auto psi0 = MPS(state);
auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet=",true});
printfln("Ground State Energy = %.12f",energy);