Hi,

I have codes in which I perform DMRG and then compute < S^2 >.

Before I updated ITensor to the latest version, these codes were giving correct results and now they are not. I know that the previous results were correct because I benchmarked with simple models for which I can do exact diagonalization.

One simple example that may help to identify this apparent bug is the following:

Model: DMRG for S=1 antiferromagnetic Heisenberg dimer.

The groundstate should be a singlet.

I compute < Sz > and I get 0. That's good.

I compute < S^2 > and I get different than 0. That's not right.

I defined the S^2 operator as follows:

```
function S2_op(Nsites,sites)
ampo = AutoMPO()
for i in 1:Nsites
for j in 1:Nsites
ampo += 1.,"S+",i,"S-",j
ampo += 1.,"Sz",i,"Sz",j
end
ampo += -1.,"Sz",i
end
S2 = MPO(ampo,sites)
return S2
end
```

I also tried to define it in a different manner: (EDIT: this is wrong)

```
function S2_op_other(Nsites,sites)
ampo = AutoMPO()
for i in 1:Nsites
ampo += 1.,"S+*S-",i
ampo += 1.,"Sz*Sz",i
ampo += -1.,"Sz",i
end
S2 = MPO(ampo,sites)
return S2
end
```

Both cases give values that are not 0. Also, they give different results, which I think is also telling. In the first case, I get 2/3. In the second case, I get 4. None of these numbers makes any physical sense.

The exact same thing happened when I tried to create an MPO of the Sz^2 operator and compute its average.

Best,

Gonçalo Catarina