Hi,
I am trying to reproduce the paper Phase diagram of the extended Bose–Hubbard model. For which I am using the following code
//test.cpp
#include <iostream>
#include <cmath>
#include <iomanip>
#include <cmath>
#include "itensor/all.h"
#include "meshgrid.hpp"
using namespace itensor;
using namespace std;
using namespace meshgen;
// Define the parameters
const Real hop = -1; // hopping parameter
const int L = 50; // number of lattice sites
const int N = 50; // number of particles
const double alpha = 0;
const int max_occ = 3; // maximum occupancy in a particular lattice site
const int max_value = 10;
const double h = 0.01;
const int max_size = (max_value / h) + 1;
int main()
{
vector<double> u_nn(max_size, 0.0);
linspace(u_nn.begin(), u_nn.end(), 0.0, h);
// input value
int c_point = 556;
// cin >> c_point;
double V = 0.0;
double U_ONSITE = 5.0;
double U_NN = u_nn[c_point];
// ITensor code
auto sites = Boson(N, {"MaxOcc", max_occ});
// Create the Hamiltonian
auto ampo = AutoMPO(sites);
for (int i = 1; i <= N; ++i)
{
// hopping term
if (i < N)
{
ampo += hop, "Adag", i, "A", i + 1;
ampo += hop, "A", i, "Adag", i + 1;
}
// potential term
ampo += V * cos(2 * M_PI * alpha * i), "N", i;
// on-site interaction term
ampo += (1.0 / 2.0) * U_ONSITE, "Adag", i, "Adag", i, "A", i, "A", i;
// nearest-neighbor interaction term
if (i < N)
{
ampo += U_NN, "N", i, "N", i + 1;
}
}
auto H = toMPO(ampo);
vector<int> box(L, 0);
for (int i = 0; i < N; i++)
{
box[i % L] += 1;
}
// Create the initial MPS
auto state = InitState(sites);
for (int i = 1; i <= L; i++)
{
state.set(i, to_string(box[i - 1]));
}
auto psi = MPS(state);
// Set up the DMRG sweeps
auto sweeps = Sweeps(10);
sweeps.maxdim() = 10000;
sweeps.cutoff() = 1E-8;
// Run the DMRG calculation
auto [en0, psi0] = dmrg(H, psi, sweeps, {"Quiet=", true});
auto wfs = vector<MPS>(1);
wfs.at(0) = psi0;
//
// Here the Weight option sets the energy penalty for
// psi1 having any overlap with psi0
//
auto [en1, psi1] = dmrg(H, wfs, psi, sweeps, {"Quiet=", true, "Weight=", 50.0});
double ng = en1 - en0;
cout << "ng: " << ng << endl;
cout << "V: " << V << endl;
cout << "U_ONSITE: " << U_ONSITE << endl;
cout << "U_NN: " << U_NN << endl;
return 0;
}
This produces different results in two different machines. On my laptop, when I am using
g++ -std=c++17 test.cpp -o test -I/home/debamalya/cpp_libraries/ITensor/ -I/home/debamalya/cpp_libraries/meshgen/include/ -L/home/debamalya/cpp_libraries/ITensor/lib -litensor -llapack -lblas
produces output
ng: 0.0212581
V: 0
U_ONSITE: 5
U_NN: 5.56
In a cluster of open PBS, when I am using
/home/debamalya/gcc_compiler/bin/g++ -std=c++17 test.cpp -o test -I/home/debamalya/cpp_libraries/ITensor/ -I/home/debamalya/cpp_libraries/meshgen/include/ -L/home/debamalya/cpp_libraries/ITensor/lib -litensor -llapack -lblas -pthread
produces output
ng: 0.0203305
V: 0
U_ONSITE: 5
U_NN: 5.56
There is a clear difference in ‘ng’. I am confused why there appears to be a difference and which one is correct!
