# Dynamical spin structure factor (DSSF) for AFM Heisenberg model via TEBD

**URL:** <https://itensor.discourse.group/t/dynamical-spin-structure-factor-dssf-for-afm-heisenberg-model-via-tebd/1251>\
**Category:** DMRG and Numerical Methods\
**Tags:** julia\
**Created:** [October 27, 2023, 6:36am UTC](https://itensor.discourse.group/t/dynamical-spin-structure-factor-dssf-for-afm-heisenberg-model-via-tebd/1251 "2023-10-27T06:36:50Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![WekawTirso](https://avatars.discourse-cdn.com/v4/letter/w/a587f6/32.png) [@WekawTirso](https://itensor.discourse.group/u/WekawTirso)\
**Post date:** [October 27, 2023, 6:36am UTC](https://itensor.discourse.group/t/dynamical-spin-structure-factor-dssf-for-afm-heisenberg-model-via-tebd/1251/1 "2023-10-27T06:36:50Z")

</div>

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](https://arxiv.org/abs/2209.00739)

The DSSF is defined to be the Fourier transform of the correlation function

S(\boldsymbol{q}, \omega) = \; \int\_0^\infty \frac{dt}{\pi \sqrt{N}} \: \sum\_{\boldsymbol{r}} \: \cos{(\boldsymbol{q} \cdot \boldsymbol{r})} \times \Big\{ \cos{(\omega t)} \: \textrm{Re}[G(\boldsymbol{r}, t)] \: -\: \sin{(\omega t)} \: \textrm{Im}[G(\boldsymbol{r}, t)] \Big\}.

In the paper, the structure factor is approximated for a finite system as

S^{zz}(\boldsymbol{q}\_i, \omega\_j) \approx \; \frac{1}{\pi \sqrt{N}} \; \sum\_{t\_m =0}^{T\_{\textrm{max}} }\; \sum\_{\boldsymbol{r}\_n = 1}^{N} \: \cos{(\boldsymbol{q}\_i \cdot \boldsymbol{r}\_n)} \times \Big\{ \cos{(\omega\_j t\_m)} \: \textrm{Re}[G^{zz}(\boldsymbol{r}\_n, t\_m)] \: -\: \sin{(\omega\_j t\_m)} \: \textrm{Im}[G^{zz}(\boldsymbol{r}\_n, t\_m)] \Big\}

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

G(\boldsymbol{r}, t) \: \approx \: e^{i E\_0 t} \: \langle \psi\_0 | S^z\_{\boldsymbol{r}\_n} \: e^{- i \hat{H} t} \: S^z\_c \: | \psi\_0 \rangle

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:

```auto
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

```auto
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

1. 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 `E0` and `psi_gs` within the `solve\_GS1 function?

2. Is the time dependence in G(\boldsymbol{r},t) controlled in Part 3 or Part 2?

3. 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?

4. Why does this work

```auto
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?

```auto
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](https://itensor.org/support/3795/computing-correlation-function?show=3795#q3795)  
(2) [How to calculate dynamical spin structure factor S(k,w) of a spin chain?](https://itensor.org/support/681/how-to-calculate-dynamical-spin-structure-factor-spin-chain?show=681#q681)  
(3) [How to evaluate two-point correlator matrix element](https://itensor.org/support/1240/how-to-evaluate-two-point-correlator-matrix-element)  
(4) [Unsuccesful TEBD for Time Independent Fermi-Hubbard Model](https://itensor.org/support/3883/unsuccesful-tebd-for-time-independent-fermi-hubbard-model)  
(5) [ITensor Code Examples](https://itensor.github.io/ITensors.jl/stable/examples/ITensor.html)

I understand implementing the whole thing may take time but any immediate help you can give on answering the questions above is greatly appreciated!