Direct Calculation of spectral functions using symmetries¶
This example demonstrates how to use the general basis function method Op_shift_sector() to compute spectral functions directly without Fourier -transforming auto-correlation functions. The idea here is to use SciPy’s iterative solvers for sparse linear systems to calculate the action of an inverted operator on a state.
The spectral function in ground state \(|0\rangle\) can be written as:
where \(\eta\) is the boradening factor for the spectral peaks. For an operator \(A\) that has quantum numbers that differ from the ground state \(|0\rangle\) we can use Op_shift_sector() to calculate the a new ket: \(|A\rangle = A|0\rangle\). With this new ket, we can construct the Hamiltonian \(H\) in the new sector, and solve the following equation:
This equation can be solved using the bi-conjugate gradient (cf. scipy.sparse.linalg.bicg) method that is implemented in SciPy. In order to solve that equation with SciPy we define a custom class (see code below) that produces the action \((\omega \pm i\eta + E_0 - H)|u\rangle\) for an arbitrary state \(|u\rangle\). Now we can use \(|x(\omega)\rangle\) to calculate the spectral function a given value of \(\omega\) as follows:
Note that this method can also be used to calculate spectral functions for any state by replacing \(|0\rangle\rightarrow|\psi\rangle\). This may be useful for calculating spectral functions away from equilibrium by using \(|\psi(t)\rangle\) that has undergone some kind of time evolution (in the Schrödinger picture).
In the example below, we look at the Heisenberg chain (\(J=1\)) with periodic boundary conditions. The periodic boundary condition allows us to use translation symmetry to reduce the total Hilbert space into blocks labeled by the many-body momentum quantum number. We limit the chain lengths to \(L=4n\) for \(n=1,2,3...\), since this is required to have the ground state in the zero-momentum sector (for other system sizes, the ground state of the isotropic Heisenberg chain has finite momentum).
Below, we calculate two different spectral functions:
where we have defined:
Because the model has full SU(2) symmetry, we expect that the two spectral functions should give the same result. We also exclude the spectral function for \(q=L/2\) (\(\pi\)-momentum) as the spectral peak is very large due to the quasi long-range antiferromagnetic order in the ground state. The variable on_the_fly can be used to switch between using a hamiltonian object to a quantum_LinearOperator object to calculate the matrix vector product. One can also change the total spin of the local spin operators by changing the variable S to compute the spectral function for higher-spin Heisenberg models.
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 | 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 26 #
# This example shows how to use the `Op-shit_sector` method of the #
# general basis class to compute spectral functions using symmetries. #
#######################################################################
from quspin.basis import spin_basis_general
from quspin.operators import hamiltonian,quantum_LinearOperator
import scipy.sparse as sp
import numexpr,cProfile
import numpy as np
import matplotlib.pyplot as plt
#
#
# define custom LinearOperator object that generates the left hand side of the equation.
#
class LHS(sp.linalg.LinearOperator):
#
def __init__(self,H,omega,eta,E0,kwargs={}):
self._H = H # Hamiltonian
self._z = omega +1j*eta + E0 # complex energy
self._kwargs = kwargs # arguments
#
@property
def shape(self):
return (self._H.Ns,self._H.Ns)
#
@property
def dtype(self):
return np.dtype(self._H.dtype)
#
def _matvec(self,v):
# left multiplication
return self._z * v - self._H.dot(v,**self._kwargs)
#
def _rmatvec(self,v):
# right multiplication
return self._z.conj() * v - self._H.dot(v,**self._kwargs)
#
##### calculate action without constructing the Hamiltonian matrix
#
on_the_fly = False # toggles between using a `hamiltnian` or a `quantum_LinearOperator` object
# chain length
L = 12
# on-site spin size
S = "1/2"
# translation transformation on the lattice sites [0,...,L-1]
T = (np.arange(L)+1)%L
# this example does not work under these conditions because ground-state sector is not q=0
if (L//2)%2 != 0:
raise ValueError("Example requires modifications for Heisenberg chains with L=4*n+2.")
if L%2 != 0:
raise ValueError("Example requires modifications for Heisenberg chains with odd number of sites.")
# construct basis
basis0 = spin_basis_general(L,S=S,m=0,pauli=False,kblock=(T,0))
# construct static list for Heisenberg chain
Jzz_list = [[1.0,i,(i+1)%L] for i in range(L)]
Jxy_list = [[0.5,i,(i+1)%L] for i in range(L)]
static = [["zz",Jzz_list],["+-",Jxy_list],["-+",Jxy_list]]
# construct operator for Hamiltonian in the ground state sector
if on_the_fly:
H0 = quantum_LinearOperator(static,basis=basis0,dtype=np.float64)
else:
H0 = hamiltonian(static,[],basis=basis0,dtype=np.float64)
# calculate ground state
[E0],psi0 = H0.eigsh(k=1,which="SA")
psi0 = psi0.ravel()
# define all possible momentum sectors excluding q=L/2 (pi-momentum) where the peak is abnormally large
qs = np.arange(-L//2+1,L//2,1)
# define frequencies to calculate spectral function for
omegas = np.arange(0,4,0.05)
# spectral peaks broadening factor
eta = 0.1
# allocate arrays to store data
Gzz = np.zeros(omegas.shape+qs.shape,dtype=np.complex128)
Gpm = np.zeros(omegas.shape+qs.shape,dtype=np.complex128)
# loop over momentum sectors
for j,q in enumerate(qs):
print("computing momentum block q={}".format(q) )
# define block
block = dict(qblock=(T,q))
#
####################################################################
#
# ---------------------------------------------------- #
# calculation but for SzSz #
# ---------------------------------------------------- #
#
# define operator list for Op_shift_sector
f = lambda i:np.exp(-2j*np.pi*q*i/L)/np.sqrt(L)
Op_list = [["z",[i],f(i)] for i in range(L)]
# define basis
basisq = spin_basis_general(L,S=S,m=0,pauli=False,**block)
# define operators in the q-momentum sector
if on_the_fly:
Hq = quantum_LinearOperator(static,basis=basisq,dtype=np.complex128,
check_symm=False,check_pcon=False,check_herm=False)
else:
Hq = hamiltonian(static,[],basis=basisq,dtype=np.complex128,
check_symm=False,check_pcon=False,check_herm=False)
# shift sectors
psiA = basisq.Op_shift_sector(basis0,Op_list,psi0)
#
### apply vector correction method
#
# solve (z-H)|x> = |A> solve for |x> using iterative solver for each omega
for i,omega in enumerate(omegas):
lhs = LHS(Hq,omega,eta,E0)
x,*_ = sp.linalg.bicg(lhs,psiA)
Gzz[i,j] = -np.vdot(psiA,x)/np.pi
#
#####################################################################
#
# ---------------------------------------------------- #
# same calculation but for S-S+ #
# ---------------------------------------------------- #
#
# divide by extra sqrt(2) to get extra factor of 1/2 when taking sandwich: needed since H = 1/2 (S^+_i S^-_j + h.c.) + S^z_j S^z_j
f = lambda i:np.exp(-2j*np.pi*q*i/L)*np.sqrt(1.0/(2*L))
Op_list = [["+",[i],f(i)] for i in range(L)]
# change S_z projection up by S for action of S+ operator
S_z_tot = 0 + eval(S)
# calculate magnetization density from S_z_tot as defined in the documentation
# m = S_z_tot / (S * N), S: local spin (as number), N: total spins
m = S_z_tot/(eval(S)*L)
basisq = spin_basis_general(L,S=S,m=m,pauli=False,**block)
# define operators in the q-momentum sector
if on_the_fly:
Hq = quantum_LinearOperator(static,basis=basisq,dtype=np.complex128,
check_symm=False,check_pcon=False,check_herm=False)
else:
Hq = hamiltonian(static,[],basis=basisq,dtype=np.complex128,
check_symm=False,check_pcon=False,check_herm=False)
# shift sectors
psiA = basisq.Op_shift_sector(basis0,Op_list,psi0)
#
### apply vector correction method
#
# solve (z-H)|x> = |A> solve for |x> using iterative solver for each omega
for i,omega in enumerate(omegas):
lhs = LHS(Hq,omega,eta,E0)
x,*_ = sp.linalg.bicg(lhs,psiA)
Gpm[i,j] = -np.vdot(psiA,x)/np.pi
#
##### plot results
#
ks = 2*np.pi*qs/L # compute physical momentum values
f,(ax1,ax2) = plt.subplots(2,1,figsize=(3,4),sharex=True)
ax1.pcolormesh(ks,omegas,Gzz.imag,shading='nearest')
ax2.pcolormesh(ks,omegas,Gpm.imag,shading='nearest')
ax2.set_xlabel('$k$')
ax2.set_ylabel('$\\omega$')
ax1.set_ylabel('$\\omega$')
ax1.set_title('$G_{zz}(\\omega,k)$')
ax2.set_title('$G_{+-}(\\omega,k)$')
#plt.show()
plt.close()
|