This may be a basic question, but I am trying to implement a thermal state calculation, and I am definitely doing something incorrectly. I think my problem lies in using AutoMPO.
At one point, I am calculating the operator S^zk (the Fourier transform of S^zi), and I build S^z_k using AutoMPO as
        auto ampo = AutoMPO(sites);
        for(auto s : lattice->getSites())
        {
            if(s.t == Lattice::physical)
            {
               ampo += sin(0.5*(s.s+1)*q),"Sz",s.s;
            }
        }
sites is of type SpinHalf with 2N sites, s.t is a tag as either "physical" or "environment", and s.s is the site number. all physical sites are odd, and ancilla sites are even.
I also need to calculate the Louivillian defined by L = HxI - IxH, with H the Heisenberg Hamiltonian. My code for this is
    auto ampo = AutoMPO(sites);
    for(auto b : lattice->getBonds())
    {
        if(b.t == Lattice::physical)
        {
            ampo += 0.5,"S+",b.s1,"S-",b.s2;
            ampo += 0.5,"S-",b.s1,"S+",b.s2;
            ampo +=     "Sz",b.s1,"Sz",b.s2;
        }
        if(b.t == Lattice::environment)
        {
            ampo += -0.5,"S+",b.s1,"S-",b.s2;
            ampo += -0.5,"S-",b.s1,"S+",b.s2;
            ampo += -1.0,"Sz",b.s1,"Sz",b.s2;
        }
    }
    L = toMPO(ampo);
My main question is, can I recycle AutoMPO for use with a thermal state or do I need to modify something else to do this? And how do I correctly implement tensoring with the identity on the ancilla degrees of freedom?
Thank you for your time,
Nick