on purification scheme

I am trying to find the thermal state using a purification scheme for the following model

H{=}\sum_{j}^{N-1}{\bf S}_{j}.{\bf S}_{j{+}1}+ h \sum_{j}^{N} S^z_{j}

The following implementation gives the correct benchmark results only for two values of h=0, 1

 gates = ITensor[]      
 for j=1:N-1
 s1 = s[j]
 s2 = s[j+1] 
 hj =       op("Sz",s1) * op("Sz",s2) +
     1/2 * op("S+",s1) * op("S-",s2) +
     1/2 * op("S-",s1) * op("S+",s2) -
     h * op("Sz",s1) * op("Id",s2) 

 Gj = exp(-1.0* tau/2 * hj)
push!(gates,Gj)
end

hj = - h * op("Id",s[N-1]) * op("Sz",s[N]) 
Gj = exp(-1.0 * tau/2 * hj)
push!(gates,Gj)

append!(gates,reverse(gates)) 

@miles Is the above implementation of trotter gates correct for the H ?. Could you please help me with this ?

For implementing the purification method for the 1D Heisenberg model, we actually give a fully working example in this example file in the ITensorMPS.jl repo:
https://github.com/ITensor/ITensorMPS.jl/blob/main/examples/finite_temperature/purification.jl
The only difference from your case is that you would need to add a magnetic field term. To make it fit the pattern of that code, you could add a term like

\frac{h}{2} S^z_i + \frac{h}{2} S^z_{i+1}

to each two-site local Hamiltonian operator.

Also, please do not address questions directly to one of us, but instead post your questions in a general way. There are others here who are also answering questions frequently.