Unable to reproduce DMRG calculations using ITensor

DMRG and Numerical Methods
,
1

Hi,
I am a beginner to DMRG, and for the learning purpose, I have started with the paper(Phase diagram of the extended Bose–Hubbard model; DOI 10.1088/1367-2630/14/6/065012), which I intend to reproduce using ITensor. The model extended Bose-Hubbard Hamiltonian, which has the form

H=-J\sum_ib_i^{\dagger}b_i+h.c+\frac{U}{2}\sum_in_i(n_i-1)+V\sum_in_in_{i+1}.

The calculation has been done with J=1, the maximum number of particles per site n_{max}=3 and the average filling \bar{n}=1. The neutral gap and charge gap are the easiest to calculate here. For this purpose, I have used 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;

const int precision = 10;
const complex<double> z(0, 1);

const double hop = -1;
const int N = 50;      // number of particle
const int L = 50;      // number of lattice sites
const int max_occ = 3; // maximum occupancy per site

const int max_value = 10;
const double h = 0.01;
const int max_size = (max_value / h) + 1;

class hamiltonian
{
private:
    int L, N, max_occ;

    double U_ONSITE, U_NN;

    MPS psi;
    MPO H;
    Sweeps sweeps;

public:
    Real en0, en1;
    MPS psi0, psi1;

    // number of lattice sites
    // number of particles
    // maximum occupancy
    // onsite interaction
    // nearest neighour interaction

    hamiltonian(int m, int n, int o, double j, double k) : U_ONSITE(j), U_NN(k), L(m), N(n), max_occ(o)
    {
        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;
            }
            // 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;
            }
        }
        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]));
        }
        psi = MPS(state);
    }

    void ground_state()
    {
        sweeps = Sweeps(20);
        sweeps.maxdim() = 10, 20, 30, 40, 50, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 330, 360, 400;
        sweeps.niter() = 6;
        sweeps.noise() = 0.0;
        sweeps.cutoff() = 1E-8;

        // Run the DMRG calculation
        auto [en_g, psi_g] = dmrg(H, psi, sweeps, {"Quiet=", true});
        en0 = en_g;
        psi0 = psi_g;
    }

    void excited_state()
    {
        sweeps = Sweeps(50);
        sweeps.maxdim() = 10, 10, 10, 10, 10, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 160, 160, 160, 160, 160, 300, 300, 600, 600, 1200, 1200, 2400, 2400, 4800, 4800;
        sweeps.niter() = 6;
        sweeps.noise() = 0.0;
        sweeps.cutoff() = 1E-8;

        auto wfs = vector<MPS>(1);
        wfs.at(0) = psi0;
        auto [en_e, psi_e] = dmrg(H, wfs, psi, sweeps, {"Quiet=", true, "Weight=", 50.0});
        en1 = en_e;
        psi1 = psi_e;
    }
};

// number of lattice sites
// number of particles
// maximum occupancy
// onsite interaction
// nearest neighour interaction

class nutral_gap : public hamiltonian
{

public:
    double ng;
    nutral_gap(int m, int n, int o, double j, double k) : hamiltonian(m, n, o, j, k)
    {
        ground_state();
        excited_state();
        ng = en1 - en0;
    }
};

class charge_gap
{
public:
    double cg;
    charge_gap(int m, int n, int o, double j, double k)
    {
        hamiltonian Hl(m, n, o, j, k);
        hamiltonian Hll(m - 1, n - 1, o, j, k);
        hamiltonian HlL(m + 1, n + 1, o, j, k);
        Hl.ground_state();
        Hll.ground_state();
        HlL.ground_state();

        cg = HlL.en0 + Hll.en0 - (2 * Hl.en0);
    }
};

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;
    cin >> c_point;

    double U_ONSITE = 5.0;
    double U_NN = u_nn[c_point];

    nutral_gap NG(L, N, max_occ, U_ONSITE, U_NN);
    charge_gap CG(L, N, max_occ, U_ONSITE, U_NN);

    double ng = NG.ng;
    double cg = CG.cg;

    cout << setprecision(precision) << "ng: " << ng << endl;
    cout << setprecision(precision) << "cg: " << cg << endl;
    cout << "U_ONSITE: " << U_ONSITE << endl;
    cout << "U_NN: " << U_NN << endl;

    // const int column_width = 20;

    // ofstream file;

    // file.open("nutral_gap/" + to_string(c_point) + ".txt");
    // if (file.is_open())
    // {
    //     file << left << setw(column_width) << "V" << left << setw(column_width) << "U_ONSITE" << left << setw(column_width) << "U_NN";
    //     file << left << setw(column_width) << "nutral_gap";
    //     file << left << setw(column_width) << "C_POINT" << endl;

    //     file << left << setw(column_width) << V << left << setw(column_width) << U_ONSITE << left << setw(column_width) << U_NN;
    //     file << left << setw(column_width) << ng;
    //     file << left << setw(column_width) << c_point << endl;
    // }
    // file.close();

    // file.open("charge_gap/" + to_string(c_point) + ".txt");
    // if (file.is_open())
    // {
    //     file << left << setw(column_width) << "V" << left << setw(column_width) << "U_ONSITE" << left << setw(column_width) << "U_NN";
    //     file << left << setw(column_width) << "charge_gap";
    //     file << left << setw(column_width) << "C_POINT" << endl;

    //     file << left << setw(column_width) << V << left << setw(column_width) << U_ONSITE << left << setw(column_width) << U_NN;
    //     file << left << setw(column_width) << cg;
    //     file << left << setw(column_width) << c_point << endl;
    // }
    // file.close();
    return 0;
}

For each different value of c_point, one can obtain each point of the graph. The U in the Hamiltonian is same U_ONSITE in code, and similarly, V in the Hamiltonian is same as U_NN in code. What I obtained is the following(left-hand image).

ng
ng1252×481 206 KB

I have no idea why the results do not resemble Figure 6 of that paper(right-hand image)!

The charge gap is also not matching, and it got negative values. I have done all the calculations for 50 particles in 50 lattice sites. A comment or a suggestion would be a great help to me.

2

I can see why you have a question, but it’s hard for me to say or guess why you aren’t getting the answer you’d expect, even when looking at your code.

We keep a list of techniques you can try to ensure that DMRG calculations are converged:
https://itensor.github.io/ITensors.jl/dev/faq/DMRG.html
These techniques can all be very helpful to consider carefully.

Another “under-rated” technique to understand DMRG convergence is to make real-space plots of the properties of your state. For example, if you plot the density of particles of your excited states, you may (or may not) see that the density is not as symmetric as you thought, and could even show evidence that a particle has formed a kind of “bound state” with the open edges of the system. In that case, you could get better results by starting with a different initial state or by changing the Hamiltonian near the boundaries of your system to make it repel particles from the boundary. So all sorts of such analyses can be needed for difficult-to-converge systems or systems that exhibit confusing or surprising behavior.