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:
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:
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:
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¶
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()
#
|