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#
  2import sys, os
  3
  4os.environ["KMP_DUPLICATE_LIB_OK"] = (
  5    "True"  # uncomment this line if omp error occurs on OSX for python 3
  6)
  7os.environ["OMP_NUM_THREADS"] = "1"  # set number of OpenMP threads to run in parallel
  8os.environ["MKL_NUM_THREADS"] = "1"  # set number of MKL threads to run in parallel
  9#
 10
 11#######################################################################
 12#                            example 21                               #
 13# This example shows how to use the `Lanczos` submodule of the        #
 14# `tools` module to compute finite temperature expectation values     #
 15# using `FTLM_statc_iteration` and `LTLM_statiic_iteration`.          #
 16#######################################################################
 17from quspin.basis import spin_basis_1d
 18from quspin.operators import hamiltonian, quantum_operator
 19from quspin.tools.lanczos import (
 20    lanczos_full,
 21    lanczos_iter,
 22    FTLM_static_iteration,
 23    LTLM_static_iteration,
 24)
 25import numpy as np
 26import matplotlib.pyplot as plt
 27from mpl_toolkits.axes_grid1.inset_locator import inset_axes
 28
 29#
 30np.random.seed(1203901)  # fix seed
 31
 32
 33#
 34#
 35def bootstrap_mean(O_r, Id_r, n_bootstrap=100):
 36    """
 37    Uses boostraping to esimate the error due to sampling.
 38
 39    O_r: numerator
 40    Id_r: denominator
 41    n_bootstrap: bootstrap sample size
 42
 43    """
 44    O_r = np.asarray(O_r)
 45    Id_r = np.asarray(Id_r)
 46    #
 47    avg = np.nanmean(O_r, axis=0) / np.nanmean(Id_r, axis=0)
 48    n_Id = Id_r.shape[0]
 49    # n_N = O_r.shape[0]
 50    #
 51    i_iter = (np.random.randint(n_Id, size=n_Id) for i in range(n_bootstrap))
 52    #
 53    bootstrap_iter = (
 54        np.nanmean(O_r[i, ...], axis=0) / np.nanmean(Id_r[i, ...], axis=0)
 55        for i in i_iter
 56    )
 57    diff_iter = ((bootstrap - avg) ** 2 for bootstrap in bootstrap_iter)
 58    err = np.sqrt(sum(diff_iter) / n_bootstrap)
 59    #
 60    return avg, err
 61
 62
 63#
 64def get_operators(L):
 65    """
 66    Generates hamiltonian for TFIM, see quspin tutorial papers to learn more aabout this
 67
 68    """
 69    # create basis
 70    basis = spin_basis_1d(L, pauli=True)
 71    # site-coupling lists
 72    J_list = [[-1.0, i, (i + 1) % L] for i in range(L)]
 73    h_list = [[-1.0, i] for i in range(L)]
 74    M_list = [[1.0 / L, i] for i in range(L)]
 75    # create magnetization-squared operator
 76    M = hamiltonian([["z", M_list]], [], basis=basis, dtype=np.float64)
 77    M2 = M**2
 78    # create parameter-dependent Hamiltonian using quantum_oprator
 79    ops_dict = dict(J=[["zz", J_list]], h=[["x", h_list]])
 80    H = quantum_operator(ops_dict, basis=basis, dtype=np.float64)
 81    #
 82    return M2, H
 83
 84
 85#
 86class lanczos_wrapper(object):
 87    """
 88    Class that contains minimum requirments to use Lanczos.
 89
 90    Using it is equired, since the dot and dtype methods of quantum_operator objects take more parameters
 91
 92    """
 93
 94    #
 95    def __init__(self, A, **kwargs):
 96        """
 97        A: array-like object to assign/overwrite the dot and dtype objects of
 98        kwargs: any optional arguments used when overwriting the methods
 99
100        """
101        self._A = A
102        self._kwargs = kwargs
103
104    #
105    def dot(self, v, out=None):
106        """
107        Calls the `dot` method of quantum_operator with the parameters fixed to a given value.
108
109        """
110        return self._A.dot(v, out=out, pars=self._kwargs)
111
112    #
113    @property
114    def dtype(self):
115        """
116        The dtype attribute is required to figure out result types in lanczos calculations.
117
118        """
119        return self._A.dtype
120
121
122#
123#
124##### define system parameters #####
125#
126L = 10  # system size
127m = 50  # dimensio of Krylov space
128s = 0.5  # transverse-field Ising model parameter: H = sZZ + (1-s)X
129#
130N_samples = 100  # of samples to approximate thermal expectation value with
131#
132T = np.logspace(-3, 3, 51, base=10)  # temperature vector
133beta = 1.0 / (T + 1e-15)  # inverse temperature vector
134#
135##### get operators #####
136#
137M2, H = get_operators(L)
138# crate wrapper for quantum_operator
139H_wrapped = lanczos_wrapper(H, J=s, h=(1 - s))
140# calculate ground state energy to use as shift that will prevent overflows (i.e. numerical instabilities)
141[E0] = H.eigsh(k=1, which="SA", pars=dict(J=s, h=1 - s), return_eigenvectors=False)
142#
143##### finite temperature methods #####
144#
145# preallocate lists to store results from iterations
146M2_FT_list = []
147M2_LT_list = []
148Z_FT_list = []
149Z_LT_list = []
150#
151# allocate memory for lanczos vectors
152out = np.zeros((m, H.Ns), dtype=np.float64)
153#
154# calculate iterations
155for i in range(N_samples):
156    # generate normalized random vector
157    r = np.random.normal(0, 1, size=H.Ns)
158    r /= np.linalg.norm(r)
159    # get lanczos basis
160    E, V, lv = lanczos_full(H_wrapped, r, m, eps=1e-8, full_ortho=True)
161    # E,V,lv = lanczos_full(H_wrapped,r,m,eps=1e-8,full_ortho=False)
162    # E,V,lv = lanczos_iter(H_wrapped,r,m,eps=1e-8)
163    # shift energy to avoid overflows
164    E -= E0
165    # calculate iteration
166    results_FT, Id_FT = FTLM_static_iteration({"M2": M2}, E, V, lv, beta=beta)
167    results_LT, Id_LT = LTLM_static_iteration({"M2": M2}, E, V, lv, beta=beta)
168    # save results to a list
169    M2_FT_list.append(results_FT["M2"])
170    Z_FT_list.append(Id_FT)
171    M2_LT_list.append(results_LT["M2"])
172    Z_LT_list.append(Id_LT)
173#
174# calculating error bars on the expectation values
175m2_FT, dm2_FT = bootstrap_mean(M2_FT_list, Z_FT_list)
176m2_LT, dm2_LT = bootstrap_mean(M2_LT_list, Z_LT_list)
177#
178##### calculating exact results from full diagonalization #####
179#
180dim_cutoff = 2000  # Hilbert space dimension cutoff
181if H.Ns < dim_cutoff:  # Hilbert space is not too big to diagonalize on a laptop
182    #
183    # adding more points for smooth line
184    T_new = np.logspace(np.log10(T.min()), np.log10(T.max()), 10 * len(T))
185    beta_new = 1.0 / (T_new + 1e-15)
186    #
187    # full diagonaization of H
188    E, V = H.eigh(pars=dict(J=s, h=1 - s))
189    # shift energy to avoid overflows
190    E -= E[0]
191    # get boltzmann weights for each temperature
192    W = np.exp(-np.outer(E, beta_new))
193    # get diagonal matrix elements for trace
194    O = M2.matrix_ele(V, V, diagonal=True)
195    # calculate trace
196    O = np.einsum("j...,j->...", W, O) / np.einsum("j...->...", W)
197#
198#
199##### plot results #####
200#
201# setting up plot and inset
202h = 4.2  # figure aspect ratio parameter
203f, ax = plt.subplots(figsize=(1.5 * h, h))
204axinset = inset_axes(ax, width="45%", height="65%", loc="upper right")
205axs = [ax, axinset]
206#
207# plot results for FTLM and LTLM.
208for a in axs:
209    a.errorbar(T, m2_LT, dm2_LT, marker=".", label="LTLM", zorder=-1)
210    a.errorbar(T, m2_FT, dm2_FT, marker=".", label="FTLM", zorder=-2)
211    #
212    if H.Ns < dim_cutoff:  # hilbert space is not too big to diagonalize on a laptop
213        a.plot(T_new, O, label="exact", zorder=0)
214    #
215    a.set_xscale("log")
216#
217# adding space for inset by expanding x limits.
218xmin, xmax = ax.get_xlim()
219ax.set_xlim((xmin, 10 * xmax))
220ax.legend(loc="lower left")
221#
222# inset adjustment to zoom in low-temp limit.
223xmin, xmax = axinset.get_xlim()
224#
225a = -0.6
226m = np.logical_and(T >= xmin, T <= 10 ** (a))
227axinset.set_xlim((xmin, 10 ** (a + 0.1)))
228ymin = min(m2_LT[m].min(), m2_FT[m].min())
229ymax = max(m2_LT[m].max(), m2_FT[m].max())
230ywin = ymax - ymin
231boundy = 0.1 * ywin
232axinset.set_ylim((ymin - boundy, ymax + boundy))
233#
234# display plot
235#plt.tight_layout()
236plt.show()
237#