I noticed that the finite-temperature results I obtain using the purification method do not match ED calculations; however, they do match if i change the temperature by a factor of 2 (i.e., the ED results for a temperature T seem to correspond to purification results at 'temperature' T/2). I wrote a minimal code that reproduces this (including also an ED snippet as a benchmark):
Python: Select all
import numpy as np
from tenpy.networks.purification_mps import PurificationMPS
from tenpy.models.spins import SpinChain
from tenpy.algorithms import dmrg, mpo_evolution
from tenpy.algorithms.purification import PurificationApplyMPO
import scipy.sparse as sparse
from numpy.linalg import eig
import tenpy
tenpy.show_config()
L = 4
chi = 10
Jz = 1
hx = 0.5
dt = 0.05
beta = 2.
print("\nConsider the Ising chain H = Jz*Sz*Sz - hx*Sx with Jz =", Jz, ", hx =", hx,"and length L =", L)
################### ED ###################
sx = sparse.csr_matrix(np.array([[0.,1.],[1.,0.]]))
sz = sparse.csr_matrix(np.array([[1.,0.],[0.,-1.]]))
Sx = []
Sz = []
for i_site in range(L):
if i_site==0:
X=sx
Z=sz
else:
X=sparse.csr_matrix(np.eye(2))
Z=sparse.csr_matrix(np.eye(2))
for j_site in range(1,L):
if j_site==i_site:
X=sparse.kron(X,sx,'csr')
Z=sparse.kron(Z,sz,'csr')
else:
X=sparse.kron(X,np.eye(2),'csr')
Z=sparse.kron(Z,np.eye(2),'csr')
Sx.append(0.5*X)
Sz.append(0.5*Z)
H = -hx*Sx[0]
for i in range(L-1): H += Jz*Sz[i]*Sz[i+1] - hx*Sx[i+1]
e,v = np.linalg.eigh(H.todense())
V = v; Vd = v.conj().T; D = np.diag(e) ### note that H = V.dot(D.dot(Vd))
print("\nFor inverse temperature beta =", beta/2, ", the thermal state obtained with ED gives us")
U = V.dot(np.diag(np.exp(-(beta/2)*e)).dot(Vd))
print("<Sx> =", np.real([np.trace(Sx[i]*U)/np.trace(U) for i in range(L)]))
print("For inverse temperature beta =", beta, ", the thermal state obtained with ED gives us")
U = V.dot(np.diag(np.exp(-beta*e)).dot(Vd))
print("<Sx> =", np.real([np.trace(Sx[i]*U)/np.trace(U) for i in range(L)]))
################### purification and tenpy ###################
model_params = dict(L=L,bc_MPS='finite',conserve=None,S=0.5,Jx=0,Jy=0,Jz=Jz,hx=hx)
M = SpinChain(model_params)
psi = PurificationMPS.from_infiniteT(M.lat.mps_sites(),bc='finite')
psi.canonical_form()
steps = int(0.5*beta/dt)
dt = 0.5*beta/steps ### makes sure that dt is such that an integer multiple gives the desired beta
par = {
'N_steps':1,
'compression_method': 'SVD', ### same results for 'SVD' and 'variational'
'trunc_params': {"chi_max": chi,"svd_min": 1e-10,},
'dt':-1j*dt,
}
evol = mpo_evolution.ExpMPOEvolution(psi,M,par)
evol.calc_U(-1j*dt,order=2)
U = evol._U_MPO
print("\nWe will now try to derive the same result using the purification method (using tenpy)")
print("We perform imaginary time-evolution with dt =", dt)
option = {'chi_list':{0:chi}}
for i in range(steps):
evol1 = PurificationApplyMPO(psi,U[0],option)
evol2 = PurificationApplyMPO(psi,U[1],option)
evol1.run()
evol2.run()
print("After having evolved by", steps, "steps, we obtain an effective beta =", 2*steps*dt)
psi.canonical_form()
print("<Sx> =", psi.expectation_value("Sx"), "\n")
We see that the ED results at beta = 1 match the purification results at beta = 2.
Note: in the purification approach, I time-evolved over an imaginary time of steps*dt = 1. As explained in the documentation ( https://tenpy.readthedocs.io/en/latest/ ... cation_mps ), this should correspond to beta = 2. (More generally, when using the purification approach, one should only need to time-evolve with exp(-beta H/2) to obtain an effective inverse temperature beta when calculating correlations.)
If anyone has any ideas what is going on, that would be really appreciated!
