Hello,
(specific questions about implementing the TEBD method are at the bottom of this post)
since it’s been asked many times without a concrete answer, I am requesting the developers to include a complete example for calculating the structure factor of a well known model: the nearest-neighbor antiferromagnetic Heisenberg model on the square lattice. Specifically, a recent paper showing this model is Fig.4 of arXiv:2209.00739v2
The DSSF is defined to be the Fourier transform of the correlation function
In the paper, the structure factor is approximated for a finite system as
where
-
they take a square lattice of size N = L_x L_y with cylindrical boundary condition
-
\boldsymbol{q}_i run over a path around the BZ given by \big\{ K \to X \to M \to K \to \Gamma \to X \big\}.
-
G(\boldsymbol{r}, t) = \langle \psi_0 | \: \boldsymbol{\hat{S}}_{\boldsymbol{r}_n} (t) \: \boldsymbol{\hat{S}}_c (0) \: | \psi_0 \rangle is the non-equal-time, two-point (multi-site), spin-spin correlation function. For a rotationally invariant Hamiltonian, like \hat{H} = J \sum_{\langle ij \rangle} \boldsymbol{\hat{S}}_i \cdot \boldsymbol{\hat{S}}_j, we can focus on a single component of \boldsymbol{\hat{S}} so that up to a constant
where, in turn, \hat{H} |\psi_0 \rangle = E_0 |\psi_0 \rangle.
- Moreover, \boldsymbol{r} = x \, \boldsymbol{a}_1 + y \, \boldsymbol{a}_2 is the position vector with basis \boldsymbol{a}_1 = (1,0), \boldsymbol{a}_2 = (0,1) and c \equiv N/2 is defined as the center site is taken to be the origin. Hence, for sites i,n=1, 2\dots, N, if we measure x_i and y_i from the lower-left corner of the lattice, we have (x_n, y_n) = (x_i, y_i) - (x_c, y_c).
Remark: since my other goal is to apply the answer to 2D anisotropic Hamiltonians, I construct the Hamiltonian in each lattice direction separately. Also, I include every part of my code in the implicit hope of people catching mistakes and/or pointing out other efficient ways.
Attempt at DSSF code:
using ITensors
using HDF5
function dist(Start, End)
xS = Start[1]
yS = Start[2]
xE = End[1]
yE = End[2]
distance = sqrt((xE - xS)^2 + (yE - yS)^2 )
return distance
end
function find_2dPath(Start, End, N)
xS = Start[1]
yS = Start[2]
xE = End[1]
yE = End[2]
Path = [(x,y) for (x,y) in zip(LinRange(xS,xE,N), LinRange(yS,yE,N))]
return Path
end
function apply_op(ϕ::MPS, opname::String, sites, siteidx::Int)
ϕ = copy(ϕ) # Make a copy of the original state
orthogonalize!(ϕ, siteidx)
new_ϕj = op(opname,sites[siteidx]) * ϕ[siteidx] # Apply the local operator
noprime!(new_ϕj)
ϕ[siteidx] = new_ϕj
return ϕ
end
let
Nx = 4
Ny = 4
N = Nx*Ny
J = 1.0
PBC = true
# ----------
sites = siteinds("S=1/2", N; conserve_qns = true)
lattice = square_lattice(Nx, Ny; yperiodic = PBC)
println()
println("******************************")
println("DSSF of the Heisenberg model")
println("Nx=", Nx,", Ny=", Ny, ", N=", N, ", J=", J)
println()
## ===========================================================
## Part 1: the g.s. energy of H
## Output: E0, and psi_gs
## Define the Heisenberg spin Hamiltonian
os = OpSum()
for b in lattice
os .+= 0.5, "S+", b.s1, "S-", b.s2
os .+= 0.5, "S-", b.s1, "S+", b.s2
os .+= J, "Sz", b.s1, "Sz", b.s2
end
## construct the Hamiltonian
H = MPO(os,sites)
## Make an array of strings which alternate between
## "Up" and "Dn" on odd and even sites.
state = [isodd(n) ? "Up" : "Dn" for n=1:N]
## construct a random MPS psi0 with physical indices sites
psi0 = randomMPS(sites,state,20)
sweeps = Sweeps(7)
maxdim!(sweeps,100,300,300,800,800,800,1000)
cutoff!(sweeps,1E-8)
# Run the DMRG algorithm
energy,phi_gs = dmrg(H,psi0,sweeps)
println("E0 = ", energy)
println("E0/N = ", energy/N)
=#
## ===========================================================
## Part 2: TEBD, Trotter gates for time evolution.
## Input: uses psi_gs from Part 1
## Output: gates and corrs
## Make gates (1,2),(2,3),(3,4),...
gates = ITensor[]
for j=1:N-1
s1 = sites[j]
s2 = sites[j+1]
hj = J * op("Sx",s1) * op("Sx",s2)
+ J * op("Sy",s1) * op("Sy",s2)
+ J * op("Sz",s1) * op("Sz",s2)
Gj = exp(-1.0im * tau/2 * hj)
push!(gates, Gj)
end
append!(gates,reverse(gates)) ## in reverse order: (N,N-1),(N-1,N-2),...
#@show gates
# Initialize psi to be a product state (alternating up and down)
psi = productMPS(sites, n -> isodd(n) ? "Up" : "Dn")
c = div(N, 2) # center site
## --------------------
## Measure the correlation
corrs = []
nsteps = 100 # number of time steps for time evolution
t = 0.0
psi_t0 = apply_op(psi, "Sz", sites, c)
phi = apply(gates, phi; maxdim=maxdim, cutoff=cutoff)
function measure_corr(j::Int)
new_psi = apply_op(psi, "Sz", sites, j)
return inner(psi,new_psi)
end
push!(corrs,measure_corr.(collect(1:N)))
#@show corrs
## ===========================================================
## Part 3: DSSF
## Input: numkPts, kPath, tMax, dt, omegaMax, omegaPts
## Also uses the result corrs from Part 2
## Output: DSSF(kPath, omega) an array of numkPts by
## --------------------
## Part 3a: k-Path around symmetry points of the square lattice BZ
## Path: K --> X --> M --> K --> G --> X
sqKpt = (pi/2.0, pi/2.0)
sqXpt = (0.0, pi)
sqMpt = (pi, pi)
sqGpt = (0.0, 0.0)
numkPts = 101 # tot no. of k-points
kPath = [sqKpt,sqXpt,sqMpt,sqKpt,sqGpt,sqXpt] # path around BZ1
numPaths = length(kPath) -1
TotDist = 0.0
for i in 1:numPaths
TotDist += dist(kPath[i],kPath[i+1])
end
nPts = Any[] #number of k-points along each segment of the path
for i in 1:numPaths
push!(nPts, Int(round((dist(kPath[i],kPath[i+1])/TotDist)* numkPts)))
end
temp = Any[]
for i in 1:numPaths
push!(temp, find_2dPath(kPath[i],kPath[i+1], nPts[i]) )
end
kPts = collect(reduce(vcat,temp)) #array-of-tuples of size (numkPts,)
## --------------------
tMax = 6.0
dt = 0.1
t = LinRange(0.0, tMax, Int(floor(tMax/dt) + 1) )
omegaMax = 6.0
omegaPts = 61
omega = LinRange(0.0,omegaMax, omegaPts)
## calculate the DSSF
DSSF = Array{Float64}(undef, numkPts, omegaPts)
for qi in 1:numkPts
qVec = collect(kPts[qi])
for wi in 1:omegaPts
tSum = 0.0
for ti in 1:tNum
rSum = 0.0
for rj in 1:N
rVec = SqRvec(rj,Ny) - SqRvec(c,Ny) #move r to center
factor1 = cos(dot(qVec,rVec))
factor2 = cos(omega[wi]*t[ti]) * real(corrs[rj,rj])
factor3 = sin(omega[wi]*t[ti]) * imag(corrs[rj,rj])
rSum += factor1 * (factor2 - factor3)
end
tSum += rSum
end
DSSF[qi, wi] = tSum/(pi*sqrt(N))
end
end
## wite DSSFzz to a file
fileDSSFzz = h5open("DSSFzz_"*string(Nx)*"x"*string(Ny)*".h5", "w")
write(fileDSSFzz, "Dssfz", DSSF)
close(fileDSSFzz)
end
Python snippet to plot the result
pi = np.pi
#------------------
Nx, Ny = 4, 4
N = Nx*Ny
Jex = 1.0
Jnum = f'{round(Jex*100):04d}'
num_kPts = 101 # number of k-points along path
omegaPts = 61 # number of frequency points
# -------------------------------------------------------------
# Plot the DSSF along a path in the first BZ
nodeFile = open('knodes.txt', 'r')
data = [line for line in nodeFile.read().splitlines()]
nodeFile.close()
k_node = np.array(data,dtype=int)
# Read the DSSF hdf5 file
DSSFFile = "DSSFzz_" + str(Nx) + "x" + str(Ny) +".h5"
with h5py.File(DSSFFile , "r") as f:
Group_key = list(f.keys())[0]
DSSF = list(f[Group_key]) # Get the data
zz = np.array(DSSF) # convert list to array
# --------------------------------
fig, ax = plt.subplots(figsize=(7.5,3))
TitleTxt = r'$S^{zz} (\boldsymbol{q}, \omega)$ for $N=$' \
+ str(Nx) + r"$\times$" + str(Ny) + ", $J =$" + str(Jex)
ax.set_title(TitleTxt,fontsize=10)
#
kVals = np.linspace(0,num_kPts,num_kPts)
wVals = np.linspace(0,6.0,omegaPts)
xx,yy = np.meshgrid(kVals,wVals)
z_min, z_max = np.abs(zz).min(), np.abs(zz).max()
#
pcm = ax.pcolormesh(xx,yy,zz,cmap='viridis',\
norm=colors.Normalize(vmin = 1e-2, vmax = z_max), \
shading='auto',lw=0,rasterized=True)
# set limit of the plot to the limits of the data
ax.set_xticks(k_node)
label=(r'$K$',r'$X$', r'$M$', r'$K$', r'$\Gamma$', r'$X$')
ax.set_xticklabels(label)
ax.set_xlabel(r'$\boldsymbol{q}$',fontsize=10)
ax.set_ylabel(r'$\omega$' ,fontsize=10)
ax.tick_params(direction='out')
for n in range(len(k_node)): # add vertical lines at node positions
ax.axvline(x=k_node[n],linewidth=0.5, linestyle='dashed', color='k')
fig.colorbar(pcm, ax=ax)
plt.show()
plt.close()
# -------------------------------------------------------------
Questions
I am new to both Julia and iTensor
-
How do I store/pass/call the parameters and results of the DMRG code (“Part 1” of the code above)? That is, how do I share the parameters and outputs
E0andpsi_gswithin the `solve_GS1 function? -
Is the time dependence in G(\boldsymbol{r},t) controlled in Part 3 or Part 2?
-
Is there a method to choosing t-range (tMax?) and is there a way to keep track of the truncation error in time-evolution? That is, any closed-form lower-bounds rather than finite-size scaling?
-
Why does this work
hj = J * op("Sz",s1) * op("Sz",s2) +
1/2 * op("S+",s1) * op("S-",s2) +
1/2 * op("S-",s1) * op("S+",s2)
but the following does not?
hj = J * op("Sx",s1) * op("Sx",s2) +
J * op("Sy",s1) * op("Sy",s2) +
J * op("Sz",s1) * op("Sz",s2)
That is, when creating the Trotter gates, are s1 and s2 the same objects as b.s1 and b.s2 from the construction of the OpSum() of the Hamiltonian?
I merged the code I found in these posts:
(1) Computing correlation function
(2) How to calculate dynamical spin structure factor S(k,w) of a spin chain?
(3) How to evaluate two-point correlator matrix element
(4) Unsuccesful TEBD for Time Independent Fermi-Hubbard Model
(5) ITensor Code Examples
I understand implementing the whole thing may take time but any immediate help you can give on answering the questions above is greatly appreciated!