Lanczos module: finite-temperature calculations

Here we use the transverse-field Ising chain with periodic boundary conditions as an example of how to use QuSpin’s Lanczos module to perform finite-temperature calculations.

We consider the transverse-field Ising model, governed by the Hamiltonian:

\[H(s) = -s\sum_{i=0}^{L-1} \sigma_i^z\sigma_{i+1}^z - (1-s)\sum_{i=0}^{L-1} \sigma_i^x\]

with periodic boundary conditions. This system is not ordered at finite temperature and, therefore, we expect to see a paramagnetic phase for any finite temperature. We are interested in computing the finite-temperature expectation value of the squared magnetization:

\[\langle M^2 \rangle_\beta = \frac{Tr\left(e^{-\beta H/2} M^2 e^{-\beta H/2}\right)}{Tr\left(e^{-\beta H} \right)},\qquad M = \frac{1}{L}\sum_{i=0}^{L-1}\sigma^z_j,\]

where \(\beta\) is the inverse temperature.

The example script below demonstrates the use of both the finite-temperature (FTLM_static_iteration) and low-temperature (LTLM_static_iteration) Lanczos methods, which gives the user the opportunity to compare the approximate methods at high and low temperatures. Following the discussion under quspin.tools.lanczos.LTLM_static_iteration and quspin.tools.lanczos.FTLM_static_iteration, we make use of the approximation:

\[\langle O\rangle_\beta \approx \frac{\overline{\langle O\rangle_r}}{\overline{\langle I\rangle_r}},\]

where the notation \(\overline{\;\cdot\;}\) denotes the average over randomly drawn states \(|r\rangle\), \(\langle O\rangle_r=\langle r| e^{-\beta H/2} O e^{-\beta H/2}|r\rangle\), and \(I\) is the identity operator.

The first part of the script below define various helper functions and classes that are useful for the calculation. The class lanczos_wrapper is a simple class that has all the necessary attributes to calculate the Lanczos basis. The idea here is that we would like to calculate the Lanczos basis for a particular value of the parameter \(s\) without having to calculate the Hamiltonian at that fixed value every time. This calculation is accomplished by using quspin’s quantum_operator class combined with the wrapper class: the quantum operator stores both the Ising and the transverse-field Hamiltonians as two separate terms. Then, the wrapper class provides an interface that calls the dot method of the quantum operator with the parameters specified, which is required by quspin’s Lanczos methods.

To apply finite-temperature Lanczos methods, we need to average over random states. The loop in the code below is where this calculation for random state vectors occurs. First. a random vector is generated. Then, using that vector, the Lanczos basis is created. With that Lanczos basis, the calculation of the FTLM and LTLM is performed, and we store the results in a list. Finally, the results are averaged, using the bootstrap helper function to calculate the error bars. The results are plotted and compared against exact resultsobtained using exact diagonalization.

Script

download script

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
from __future__ import print_function, division
#
import sys,os
os.environ['KMP_DUPLICATE_LIB_OK']='True' # uncomment this line if omp error occurs on OSX for python 3
os.environ['OMP_NUM_THREADS']='1' # set number of OpenMP threads to run in parallel
os.environ['MKL_NUM_THREADS']='1' # set number of MKL threads to run in parallel
#
quspin_path = os.path.join(os.getcwd(),"../../")
sys.path.insert(0,quspin_path)
#######################################################################
#                            example 21                               #	
# This example shows how to use the `Lanczos` submodule of the        #
# `tools` module to compute finite temperature expectation values     #
# using `FTLM_statc_iteration` and `LTLM_statiic_iteration`.          #
#######################################################################
from quspin.basis import spin_basis_1d
from quspin.operators import hamiltonian,quantum_operator
from quspin.tools.lanczos import lanczos_full,lanczos_iter,FTLM_static_iteration,LTLM_static_iteration
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
#
np.random.seed(1203901) # fix seed
#
#
def bootstrap_mean(O_r,Id_r,n_bootstrap=100):
	"""
	Uses boostraping to esimate the error due to sampling.

	O_r: numerator
	Id_r: denominator
	n_bootstrap: bootstrap sample size

	"""
	O_r = np.asarray(O_r)
	Id_r = np.asarray(Id_r)
	#
	avg = np.nanmean(O_r,axis=0)/np.nanmean(Id_r,axis=0)
	n_Id = Id_r.shape[0]
	#n_N = O_r.shape[0]
	#
	i_iter = (np.random.randint(n_Id,size=n_Id) for i in range(n_bootstrap))
	#
	bootstrap_iter = (np.nanmean(O_r[i,...],axis=0)/np.nanmean(Id_r[i,...],axis=0) for i in i_iter)
	diff_iter = ((bootstrap-avg)**2 for bootstrap in bootstrap_iter)
	err = np.sqrt(sum(diff_iter)/n_bootstrap)
	#
	return avg,err
#
def get_operators(L):
	"""
	Generates hamiltonian for TFIM, see quspin tutorial papers to learn more aabout this

	"""
	# create basis
	basis = spin_basis_1d(L,pauli=True)
	# site-coupling lists
	J_list = [[-1.0,i,(i+1)%L] for i in range(L)]
	h_list = [[-1.0,i] for i in range(L)]
	M_list = [[1.0/L,i] for i in range(L)]
	# create magnetization-squared operator
	M = hamiltonian([["z",M_list]],[],basis=basis,dtype=np.float64)
	M2 = M**2
	# create parameter-dependent Hamiltonian using quantum_oprator
	ops_dict = dict(J=[["zz",J_list]],h=[["x",h_list]])
	H = quantum_operator(ops_dict,basis=basis,dtype=np.float64)
	#
	return M2,H
#
class lanczos_wrapper(object):
	"""
	Class that contains minimum requirments to use Lanczos. 
	
	Using it is equired, since the dot and dtype methods of quantum_operator objects take more parameters 
	
	"""
	#
	def __init__(self,A,**kwargs):
		"""
		A: array-like object to assign/overwrite the dot and dtype objects of
		kwargs: any optional arguments used when overwriting the methods

		"""
		self._A = A
		self._kwargs = kwargs
	#
	def dot(self,v,out=None):
		"""
		Calls the `dot` method of quantum_operator with the parameters fixed to a given value.

		"""
		return self._A.dot(v,out=out,pars=self._kwargs)
	#
	@property
	def dtype(self):
		"""
		The dtype attribute is required to figure out result types in lanczos calculations.

		"""
		return self._A.dtype
#
#
##### define system parameters #####
#
L = 10 # system size
m = 50 # dimensio of Krylov space
s = 0.5 # transverse-field Ising model parameter: H = sZZ + (1-s)X 
# 
N_samples = 100 # of samples to approximate thermal expectation value with
#
T = np.logspace(-3,3,51,base=10) # temperature vector
beta = 1.0/(T+1e-15) # inverse temperature vector
#
##### get operators #####
#
M2,H = get_operators(L)
# crate wrapper for quantum_operator
H_wrapped = lanczos_wrapper(H,J=s,h=(1-s))
# calculate ground state energy to use as shift that will prevent overflows (i.e. numerical instabilities)
[E0] = H.eigsh(k=1,which="SA",pars=dict(J=s,h=1-s),return_eigenvectors=False)
#
##### finite temperature methods #####
# 
# preallocate lists to store results from iterations
M2_FT_list = []
M2_LT_list = []
Z_FT_list = []
Z_LT_list = []
#
# allocate memory for lanczos vectors
out = np.zeros((m,H.Ns),dtype=np.float64)
#
# calculate iterations
for i in range(N_samples):
	# generate normalized random vector
	r = np.random.normal(0,1,size=H.Ns)
	r /= np.linalg.norm(r)
	# get lanczos basis
	E,V,lv = lanczos_full(H_wrapped,r,m,eps=1e-8,full_ortho=True)
	# E,V,lv = lanczos_full(H_wrapped,r,m,eps=1e-8,full_ortho=False)
	# E,V,lv = lanczos_iter(H_wrapped,r,m,eps=1e-8)
	# shift energy to avoid overflows
	E -= E0
	# calculate iteration
	results_FT,Id_FT = FTLM_static_iteration({"M2":M2},E,V,lv,beta=beta)
	results_LT,Id_LT = LTLM_static_iteration({"M2":M2},E,V,lv,beta=beta)
	# save results to a list
	M2_FT_list.append(results_FT["M2"])
	Z_FT_list.append(Id_FT)
	M2_LT_list.append(results_LT["M2"])
	Z_LT_list.append(Id_LT)
#
# calculating error bars on the expectation values
m2_FT,dm2_FT = bootstrap_mean(M2_FT_list,Z_FT_list)
m2_LT,dm2_LT = bootstrap_mean(M2_LT_list,Z_LT_list)
#
##### calculating exact results from full diagonalization #####
#
dim_cutoff=2000 # Hilbert space dimension cutoff  
if H.Ns < dim_cutoff: # Hilbert space is not too big to diagonalize on a laptop
	#
	# adding more points for smooth line
	T_new = np.logspace(np.log10(T.min()),np.log10(T.max()),10*len(T))
	beta_new = 1.0/(T_new+1e-15)
	#
	# full diagonaization of H
	E,V = H.eigh(pars=dict(J=s,h=1-s))
	# shift energy to avoid overflows
	E -= E[0]
	# get boltzmann weights for each temperature
	W = np.exp(-np.outer(E,beta_new))
	# get diagonal matrix elements for trace
	O = M2.matrix_ele(V,V,diagonal=True) 
	# calculate trace
	O = np.einsum("j...,j->...",W,O)/np.einsum("j...->...",W)
#
#
##### plot results #####
#
# setting up plot and inset
h=4.2 # figure aspect ratio parameter
f,ax = plt.subplots(figsize=(1.5*h,h))
axinset = inset_axes(ax, width="45%", height="65%", loc="upper right")
axs = [ax,axinset]
#
# plot results for FTLM and LTLM.
for a in axs:
	a.errorbar(T,m2_LT,dm2_LT,marker=".",label="LTLM",zorder=-1)
	a.errorbar(T,m2_FT,dm2_FT,marker=".",label="FTLM",zorder=-2)
	#
	if H.Ns < dim_cutoff: # hilbert space is not too big to diagonalize on a laptop
		a.plot(T_new,O,label="exact",zorder=0)
	#
	a.set_xscale("log")
#
# adding space for inset by expanding x limits.
xmin,xmax = ax.get_xlim()
ax.set_xlim((xmin,10*xmax))
ax.legend(loc="lower left")
#
# inset adjustment to zoom in low-temp limit.
xmin,xmax = axinset.get_xlim()
#
a = -0.6
m = np.logical_and(T>=xmin,T<=10**(a))
axinset.set_xlim((xmin,10**(a+0.1)))
ymin = min(m2_LT[m].min(),m2_FT[m].min())
ymax = max(m2_LT[m].max(),m2_FT[m].max())
ywin = ymax-ymin
boundy = 0.1*ywin
axinset.set_ylim((ymin-boundy,ymax+boundy))
#
# display plot
f.tight_layout()
#plt.show()
#