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^z*k (the Fourier transform of S^z*i), 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