Time evolution of MPO

Question

Hello,

I would like to time evolve a quantum spin system described by a density matrix (MPO). For the time evolution, I use the second order Trotter-Suzuki expansion of the time-evolution operator. This boils down to applying a set of gates to the density matrix rho as: gatesrho(gates)^\dagger.

I am wondering if there is a simple trick to use the routine gateTEvol for an MPO although it is implemented for the time evolution of an MPS.

When looking at tevol.h, where gateTEvol is defined, it seems to me that most routines used are also applicable to MPO's. Is it correct that, when rewriting tevol.h for MPO's, the main changes would be in the definition of the tensor AA or do I miss something important?

I would be very grateful for your advice,

Walter

Answer 1

Hi Walter,
Thanks for the question. So the short answer, at least for the C++ version of ITensor, is no, at least if your question is whether there is a way to just call gateTEvol and change only a few lines of code to switch it to time evolve an MPO instead of an MPS.

But our Julia version of ITensor can do exactly what you are asking (time evolve either an MPS or an MPO), and even better, can automatically handle acting with gates beyond nearest-neighbor sites. Also it's easier to make the gates and arrays of them to pass to the function (which in Julia is called apply instead of gateTEvol). So I would strongly suggest you consider that option unless it's really necessary for you to stick with C++. May I ask if that's the case or what is your reason for preferring the C++ version at this time?

Back to the C++ version, you could use gateTEvol to time-evolve your MPO with a little bit of work though. One thing you could do is to treat the "bra" sites of your MPO as "ancillas" or doubled sites and "fold" these indices over to view your MPO as a pure state MPS with twice as many indices. Then you can use SVD's to split your MPO tensors into two tensors so that you end up with an MPS that's formally identical to your MPO. Finally you can replace the primed "bra" sites with new site indices so that it follows the expected MPS conventions of how sites work in ITensor. If you go this route, you will also need to use swap gates in order to act on pairs of sites j,j+2 because the "ket" sites will be spaced two sites apart from each other and similar for the "bra" sites. Hope that all makes sense.

Or you could, as you suggest, use the code inside of gateTEvol as a starting point to writing your own custom code for MPOs, which is what I would most recommend if you do decide to stick with the C++ version.

Please let me know which option sounds best to you, or if you have any more questions.

Miles

Comments

Dear Miles,

thank you very much for your prompt response. I chose C++ because I had some experience with this language. Apart from that I wrongly assumed that the Julia version of ITensor has the same functionality as the C++ version. Knowing now that the Julia version has the function that I need to time evolve an MPO, I will consider changing to Julia in the future.

Meanwhile I chose the last option that you mentioned above and wrote a function based on gateTEvol. Preliminary tests of this code seem to suggest that it works well. I include this code below. The main modifications concern (as suspected) the ITensor AA.

Walter
============================================================
void gateTEvolMPO(std::vector<BondGate> const& gatelist,
          Real ttotal,
          Real tstep,
          MPO & rho,
          Args args)
    {
    const bool verbose = args.getBool("Verbose",false);

    const int nt = int(ttotal/tstep+(1e-9*(ttotal/tstep)));
    if(fabs(nt*tstep-ttotal) > 1E-9)
        {
        Error("Timestep not commensurate with total time");
        }

    if(verbose)
        {
        printfln("Taking %d steps of timestep %.5f, total time %.5f",nt,tstep,ttotal);
        }

    rho.position(gatelist.front().i1());

    Real tsofar = 0;
    int i1, i2, ni1, ni2;
    ITensor AA;
    for(auto tt : range1(nt))
        {
        auto g = gatelist.begin();
        while(g != gatelist.end())
            {
            i1 = g->i1();
            i2 = g->i2();
            AA = rho(i1)*rho(i2)*prime(g->gate());
            AA *= swapTags(conj(g->gate()),"Site,1","Site,0");
            AA.replaceTags("Site,1","Site,0");
            AA.replaceTags("Site,2","Site,1");

            ++g;
            if(g != gatelist.end())
                {
                //Look ahead to next gate position
                ni1 = g->i1();
                ni2 = g->i2();
                //SVD AA to restore MPO form
                //before applying current gate
                if(ni1 >= i2)
                    {
                    rho.svdBond(i1,AA,Fromleft,args);
                    rho.position(ni1); //does no work if position already ni1
                    }
                else
                    {
                    rho.svdBond(i1,AA,Fromright,args);
                    rho.position(ni2); //does no work if position already ni2
                    }
                }
            else
                {
                //No next gate to analyze, just restore MPO form
                rho.svdBond(i1,AA,Fromright,args);
                }
            }

        tsofar += tstep;

        args.add("TimeStepNum",tt);
        args.add("Time",tsofar);
        args.add("TotalTime",ttotal);
        }
    if(verbose)
        {
        printfln("\nTotal time evolved = %.5f\n",tsofar);
        }

    } // gateTEvolMPO

---

Hi Walter,
Thanks for that update and glad you modified the code and that it appears to be working as you want! Please let us know if you have more questions, and if so, please post them as new questions in the forum as you have been doing.

Thanks,
Miles