Hi,
I try creating a Hubbard model Hamiltonian in the two different ways using AutoMPO:
In the first case, I used Jordan-Wigner strings:
auto ampoGS = AutoMPO(sites);
ampoGS += -1.0,"Adagup",b,"F",b,"Aup",b+1; //Hopping for spin up
ampoGS += 1.0,"Aup",b,"F",b,"Adagup",b+1;
ampoGS += -1.0,"Adagdn",b,"F",b+1,"Adn",b+1; //Hopping for spin down
ampoGS += 1.0,"Adn",b,"F",b+1,"Adagdn",b+1;
In the second case, I used fermion operators assuming AutoMPO handles this fine.
auto ampo = AutoMPO(sites);
ampo += -1.0,"Cdagup",b,"Cup",b+1; //Hopping for spin up
ampo += -1.0,"Cup",b,"Adagup",b+1;
ampo += -1.0,"Cdagdn",b,"Cdn",b+1; //Hopping for spin down
ampo += -1.0,"Cdn",b,"Cdagdn",b+1;
Next, I found wildly different results performing the DMRG operation on the Hamiltonian constructed using HGS = IQMPO(ampoGS) vs. H = IQMPO(ampo). For free fermions at half-filling, I checked that the DMRG with HGS gives the correct results, while H does not. For the interacting case with Hubbard U, I found that dmrg method on H gives results that vary from run to run?
Ultimately, I want to do time evolution, and for this I need toExpH which flags an error when working with HGS because it contains more than "2 operator" terms. (And I can't trust "H" but maybe I am making a simple error working with AutoMPO?)
Any help would be really appreciated! Thanks!
my complete code is:
auto sites = Hubbard(N);
auto state = InitState(sites);
for(int i = 1; i <= N; ++i)
{
if(i%2 == 1)
state.set(i,"Up");
else
state.set(i,"Dn");
}
auto psioriginal = IQMPS(state);
// Create the Target Hamiltonian and find the Ground State Correlations
auto ampo = AutoMPO(sites);
auto ampoGS = AutoMPO(sites);
auto ampoLocalEnergy = AutoMPO(sites);
//AutoMPO version with C operators doesn't work correctly. Tested for ground state energy at half-filling expecting -1.27 and got something else.
for (int b = 1; b < N; b++)
{
ampoGS += -1.0,"Adagup",b,"F",b,"Aup",b+1; //Hopping for spin up
ampoGS += 1.0,"Aup",b,"F",b,"Adagup",b+1;
ampoGS += -1.0,"Adagdn",b,"F",b+1,"Adn",b+1; //Hopping for spin down
ampoGS += 1.0,"Adn",b,"F",b+1,"Adagdn",b+1;
ampo += -1.0,"Cdagup",b,"Cup",b+1; //Hopping for spin up
ampo += -1.0,"Cup",b,"Adagup",b+1;
ampo += -1.0,"Cdagdn",b,"Cdn",b+1; //Hopping for spin down
ampo += -1.0,"Cdn",b,"Cdagdn",b+1;
ampoGS += U,"Nupdn",b; //Interactions between upspin/downspin electrons on the same site
ampo += U,"Nupdn",b; //Interactions between upspin/downspin electrons on the same site
}
ampoGS += U,"Nupdn",N;
ampo += U,"Nupdn",N;
auto HGS = IQMPO(ampoGS);
auto sweeps = Sweeps(5); //number of sweeps is 5
sweeps.maxm() = 10,20,100,100,200; //gradually increase states kept
sweeps.cutoff() = 1E-10; //desired truncation error
double Finalgroundstateenergydensity = dmrg(psioriginal,HGS,sweeps,{"Quiet=",true})/N; // psioriginal is now the ground state of the System.
/* auto H = IQMPO(ampo);
// auto sweeps = Sweeps(5); //number of sweeps is 5
// sweeps.maxm() = 10,20,100,100,200; //gradually increase states kept
// sweeps.cutoff() = 1E-10; //desired truncation error
// double Finalgroundstateenergydensity = dmrg(psioriginal,H,sweeps,{"Quiet=",true})/N; // psioriginal is now the ground state of the System.
std::cout << Finalgroundstateenergydensity << std::endl;