Quick code review: Possible misunderstanding of MPO/MPS

Question

The below (< 50 lines) of code is the simplest way I could illustrate my possible misunderstanding. To be clear, this is not the actual problem I am trying to solve in my research (the solution is trivial); I'm just trying to find my way around ITensor and MPO/MPS.

Here is my current understanding of what this code should do: I construct an MPS of all "occupied" sites and two MPOs. The first is the Hamiltonian, which is supposed to subtract particles from occupied sites to yield empty sites, at a simple rate of 1. At each timestep, I apply the Hamiltonian times a small timestep @@h@@. This is the explicit Euler solution to the equation:

$$
\frac{\mathrm{d} \mathbf{y}}{\mathrm{d} t} = A \mathbf{y}
$$

The observable I want to record is the number of occupied sites, so I use the second MPO, constructed from the number operator, and apply it as @@\langle \psi, B \psi \rangle@@.

Examining the raw output, or even plotting the results, is not encouraging, because the total number of occupied states is increasing, not decreasing. I'm sure I've just misunderstood something along the way, but I can't think of any particcular reason why the method I've described wouldn't yield (something close to) the analytic solution,

$$
y(t) = 3 e^{-t}.
$$

As always, help is greatly appreciated.

#include <itensor/all_basic.h>
#include <itensor/all_mps.h>
#include <iostream>

int main(int argc, char *argv[])
{
    int N = 3;
    auto sites = itensor::Spinless(N);

    auto ampo = itensor::AutoMPO(sites);
    auto ampo2 = itensor::AutoMPO(sites);
    auto state = itensor::InitState(sites);

    for (auto j : itensor::range1(N))
    {
        ampo  += "A", j;
        ampo2 += "N", j;
        state.set(j, "Occ");
    }

    auto mpo = itensor::toMPO<itensor::ITensor>(ampo, { "Exact=", true });
    auto mpo2 = itensor::toMPO<itensor::ITensor>(ampo2, { "Exact=", true });
    auto mps = itensor::MPS(state);

    double h = 1.e-6, t = 0;

    for (auto i : itensor::range(1000))
    {
        mps.plusEq(h * itensor::applyMPO(mpo, mps, { "Exact=", true }));
        t += h;

        auto norm = itensor::overlap(mps, mps);
        auto count = itensor::overlap(mps, mpo2, mps);
        std::cout << std::scientific << t << " " << norm << " " << count
                  << std::defaultfloat << std::endl;
    }

    return 0;
}

Comments

Hi, so we think this is due to a bug in ITensor that we know about, but haven't fixed yet. We'll try to fix it by middle of next week and let you know -

Also, by the way, if you put the line: using namespace itensor; above your main, just after the #include statements, you won't have to write itensor:: before every function call.

Best regards,
Miles

---

Thanks for letting me know, Miles! Any workarounds for this bug, or is it best to wait for the fix?

---

Here is a workaround for you to make the MPO's you were trying to:

    template<typename Tensor>
    MPOt<Tensor>
    makeOpSumMPO(SiteSet const& sites,
                 std::string opname,
                 std::vector<Real> coefs = std::vector<Real>())
        {
        using IndexT = typename Tensor::index_type;
        auto N = sites.N();
        auto M = MPOt<Tensor>(sites);

        if(coefs.empty())
            {
            coefs = std::vector<Real>(1+N);
            for(auto n : range1(N)) coefs.at(n) = 1.0;
            }

        auto F = flux(sites.op(opname,1));

        std::vector<IndexT> links(N+1);
        for(auto n : range(N+1))
            {
            links.at(n) = IQIndex(format("l%d",n),Index("R"),QN(),
                                                  Index("L"),F);
            }

        for(auto j : range1(N))
            {
            auto& A = M.Aref(j);
            IndexT ll = links.at(j-1);
            IndexT rl = links.at(j);
            IndexT s = sites(j);
            Real coef = coefs.at(j);
            A = Tensor(dag(ll),dag(s),prime(s),rl);
            A += setElt(dag(ll)(1))*sites.op("Id",j)*setElt(rl(1));
            A += setElt(dag(ll)(2))*sites.op("Id",j)*setElt(rl(2));
            A += coef*setElt(dag(ll)(2))*sites.op(opname,j)*setElt(rl(1));
            }

        auto L = Tensor(links.at(0));
        L.set(2, 1.0);
        M.Aref(1) *= L;

        auto R = Tensor(dag(links.at(N)));
        R.set(1, 1.0);
        M.Aref(N) *= R;

        return M;
        }


You can use it like this:

int main()
    {
    int N = 10;
    auto sites = Spinless(N);

    auto M = makeOpSumMPO<IQTensor>(sites,"A");

    PrintData(M);

    return 0;
    }

---

That function will make an MPO (or IQMPO, depending on whether you put ITensor or IQTensor as the template parameter) which is a sum of single-site operators. You can put different coefficients on each site by passing a std::vector "coefs" of coefficients. We'll be working on making sure the AutoMPO system can do this, but for now this ought to be a good solution for you, and was quicker for me to make than the time it will take to fix this AutoMPO issue. (Btw the issue with AutoMPO, if you want to know, boils down to something like how I defined the IQIndex inside that function to have a quantum number "F" which is the "flux" of the operator being placed at each site. Right now the AutoMPO system doesn't quite handle that detail correctly, but it does work perfectly for MPOs of things like Hamiltonians.)