So I’ve been trying to use this function: expHermitian(ITensor H, Cplx tau = 1)
but the operator that I want to exponentiate is the total Sz. Here is the code that I have so far:

auto ampo2 = AutoMPO(sites); for(auto k : range1(Nspins)) { ampo2 += "Sz",k; } auto proj = toMPO(ampo2); auto exp2H = expHermitian(proj,2);

but I get an error that the function expects an ITensor object as the first argument (as expected). Is it possible to somehow do what I am trying to do?

Good question. It’s correct what the error says, that expHermitian is only for ITensor objects, not for AutoMPO objects.

Technically we do have a system for exponentiating AutoMPO’s, but that would actually be overkill here. Fortunately the operator you want to exponentiate is a sum of single-site, commuting operators. So the exponential of the sum turns into a product of the exponentials. This means you can just prepare an MPO by setting each MPO tensor one by one to be the exponential of Sz on that site.

Here is a rough draft of a code – I haven’t tested it – but you can adapt it for your case.

auto expSz = MPO(N);
for(auto j : range1(N))
{
auto Szj = op(sites,"Sz");
expSz.ref(j) = expHermitian(Szj);
}

Of course you might want to make \exp(\tau \sum_j S^z_j) for some \tau so you can use the optional Cplx tau argument to expHermitian to do that.

Awesome, thanks a lot! (I suppose you meant auto Szj = op(sites,“Sz”, j); , right?)

But I am having another question now… If I want to add that new state to a state that I created and was real, what should I do? I tried converting it to cmplx, by multiplying with 1+0_i, but I had no luck… I am using plusEq to add them.

Hmm actually the issue is more complicated than I thought. Below I have the code I am having so far:

auto state1 = InitState(sites);
state1.set(1,"Dn");
state1.set(2,"Up");
state1.set(3,"Dn");
state1.set(4,"Up");
auto psi1 = MPS(state1);
auto expSz = MPO(N);
for(auto k : range1(N))
{
auto Szk = op(sites,"Sz", k);
expSz.ref(k) = expHermitian(Szk,num);
}
auto psi2 = applyMPO(expSz, psi1);
psi1.plusEq(psi2);

and below is the error I am getting:

L = ITensor ord=2: (2|id=286|n=1,Site,S=1/2)’ (2|id=598|l=1,Link)'1254313
{norm=1.00 (Dense Real)}
R = ITensor ord=2: (2|id=286|n=1,Site,S=1/2) (2|id=598|l=1,Link)'1254313
{norm=1.00 (Dense Real)}
From line 1021, file itensor.cc
ITensoITensor::operator+=: different index structure
ITensoITensor::operator+=: different index structure

Essentially the issue was not only with a complex number in the exponential but more general.
Any ideas?

If you look more closely, the site indices of the second MPS have a primelevel of 1. (A single prime on them.) It is a little bit subtle and easy to miss!

So the line you need here is

psi2.noPrime();

before the plusEq line.

The reason for the primes may be clear to you, but if not, if you draw the diagram for the MPO*MPS contraction, the MPO has uncontracted primed indices coming out of the top, so afterward the MPS returned has these primed indices. (We’ve cleaned this up in the Julia version a bit where we have two styles of “apply” functions, one which removes the primes for you and another which is more literal.)