Greetings everyone,
I intend to calculate the time dependent correlation function with TDVP algorithm. I applied the 2-site TDVP. The result is a bit strange, since the correlation function is oscillating (negative and positive values). Is this normal?
C_{jj}(t) = \bra{\Psi} \hat{j}(t)\hat{j}(0) \ket{\Psi} = \bra{\Psi} e^{i\hat{H}t}\hat{j}(0)e^{-i\hat{H}t} \hat{j}(0) \ket{\Psi}
\\
\mathcal{D}_s(t) = \frac{1}{\chi_0}\int_0^{t_{max}}C_{jj}(t) dt
where j is the spin current operator. The code below is for a system of 2D spin -1/2 particles, in the presence of random disorder. For the case of Ly=1, I get perfect result as for the case of 1D chain.
# ==========================================================================
# ==========================================================================#
using ITensors
using LinearAlgebra
using Random
using Plots
using Printf
using Distributions
using Statistics
using ITensorMPS
using ITensorTDVP
# ==========================================================================
# Parameters
include("geometry.jl")
const J = 1.0
const Δ = 1.0
const W = 0.0
const Lx = 8
const Ly = 1
const N = Lx * Ly
const tmax = 100.0
const dt = 0.01
const η = 0.1
const χmax = 512
rng = MersenneTwister(1234) # For reproducibility
#===========================================================================#
# ========================= Lattice geometry and site indices ==============#
function coord_to_index(x::Int, y::Int, Lx::Int, Ly::Int)
return (y - 1) * Lx + x
end
sites = siteinds("S=1/2", N, conserve_qns=true)
#println("Lattice site numbering (y goes up, x goes right):")
for y in Ly:-1:1
for x in 1:Lx
idx = coord_to_index(x, y, Lx, Ly)
end
end
# ==========================================================================#
# ========================= Hamiltonian MPO construction ===================#
function build_hamiltonian_mpo(sites, Lx, Ly, J, Δ, W, xperiodic = true, yperiodic = false, rng = Random.default_rng())
hos = OpSum()
# On-site disorder
if W > 0.0
for y in 1:Ly
for x in 1:Lx
i = coord_to_index(x,y,Lx,Ly)
w_i = rand(rng, Uniform(-W, W))
hos += w_i, "Sz", i
end
end
end
# x-bonds
for y in 1:Ly
for x in 1:(Lx - 1)
i = coord_to_index(x, y, Lx, Ly)
j = coord_to_index(x + 1, y, Lx, Ly)
hos += (J/2), "S+", i, "S-", j
hos += (J/2), "S-", i, "S+", j
hos += Δ * J, "Sz", i, "Sz", j
end
if xperiodic
i = coord_to_index(Lx, y, Lx, Ly)
j = coord_to_index(1, y, Lx, Ly)
hos += (J/2), "S+", i, "S-", j
hos += (J/2), "S-", i, "S+", j
hos += Δ * J, "Sz", i, "Sz", j
end
end
# y-bonds
if Ly > 1
for x in 1:Lx
for y in 1:(Ly-1)
i = coord_to_index(x, y, Lx, Ly)
j = coord_to_index(x, y+1, Lx, Ly)
hos += (J/2), "S+", i, "S-", j
hos += (J/2), "S-", i, "S+", j
hos += Δ * J, "Sz", i, "Sz", j
end
if yperiodic
i = coord_to_index(x, Ly, Lx, Ly)
j = coord_to_index(x, 1, Lx, Ly)
hos += (J/2), "S+", i, "S-", j
hos += (J/2), "S-", i, "S+", j
hos += Δ * J, "Sz", i, "Sz", j
end
end
end
# Convert to MPO and return!
return MPO(hos, sites)
end
H = build_hamiltonian_mpo(sites, Lx, Ly, J, Δ, W, true, false, rng)
# ==========================================================================#
# ========================= Spin current operator MPO construction ===========#
function build_spin_current_mpo(sites, Lx, Ly, J, xperiodic = true)
jos = OpSum()
for y in 1:Ly
for x in 1:(Lx - 1)
i = coord_to_index(x, y, Lx, Ly)
j = coord_to_index(x + 1, y, Lx, Ly)
# j^s_x = (iJ/2) [S^+_i S^-_j - S^-_i S^+_j]
jos += (im * J / 2), "S+", i, "S-", j
jos += (-im * J / 2), "S-", i, "S+", j
end
if xperiodic
i = coord_to_index(Lx, y, Lx, Ly)
j = coord_to_index(1, y, Lx, Ly)
jos += (im * J / 2), "S+", i, "S-", j
jos += (-im * J / 2), "S-", i, "S+", j
end
end
return MPO(jos, sites)
end
Jx = build_spin_current_mpo(sites, Lx, Ly, J, true)
# ==========================================================================#
# ========================= Random initial state and time evolution ===========#
# Random state |Ψ⟩ with bond dimension 1 (product state)
println("Generating random initial product state |Ψ⟩ (χ=1) in Sz=0 sector...")
#1. Create a random product state with equal numbers of up and down spins
half_N = div(N, 2)
basis_state = vcat(fill("↑", half_N), fill("↓", N - half_N))
#2. Shuffle the state to create a random configuration
shuffle!(rng, basis_state)
#3. Convert the random basis state to an MPS with bond dimension 1
Ψ = randomMPS(rng,sites, basis_state;linkdims=1)
normalize!(Ψ)
#4. Check the total Sz to ensure it's zero (should be zero by construction)
#sz_values = expect(Ψ, "Sz")
#total_sz = sum(sz_values)
#println("Total Sz of initial state:", total_sz)
#normm = norm(Ψ)
#println("Initial state norm:", normm)
#println("inner product of Ψ with itself:", inner(Ψ, Ψ))
ψ0 = random_mps(sites, basis_state; linkdims = 10)
# ==========================================================================#
edmrg, ψdmrg = dmrg(H, ψ0; nsweeps=10, maxdim=[10, 20, 50,50,100, 100], cutoff=1e-10)
println("Ground state energy from DMRG: ", edmrg)
# =============== Acting with the current operator in the random state =====#
println("inner Ψ with itself before applying Jx:", inner(Ψ, Ψ))
println("Acting with the current operator Jx on |Ψ⟩ to get |Φ⟩ = J_{x}|Ψ⟩...")
Φ = apply(Jx,Ψ;alg="naive",cutoff=1e-8)
#orthogonalize!(Φ)
normalize!(Φ)
println("max bond dimension of |Φ⟩ after applying Jx:", maximum(linkdims(Φ)))
println("overlap of |Φ⟩ with itself:",inner(Φ,Φ))
# ===========================================================================#
# ========================= Time evolution of |Φ⟩ and |Ψ⟩ under H ===========#
# ========================= Define TDVP evolution parameters ================#
tdvp_kwargs = (
nsite = 2, # use the 2 site algorithm
nsweeps = 1, # number of full sweeps per time step dt
maxdim = χmax, # allow bond dimension to grow to χmax
cutoff = 1e-10, # SVD truncation error threshold
normalize = true # prserve physical norm during evolution (we will normalize manually after each step
)
# ========================= Data containers for observables =================#
times = Float64[]
Ct_values = ComplexF64[]
Ct_abs = Float64[]
# ========================== Record t=0 data ================================#
# we use Ψ' to make sure physical indices match the MPO legs
# C(0) = ⟨Ψ|J_x|Φ⟩ = ⟨Ψ|J_x²|Ψ⟩/||J_x|Ψ⟩|| = ||J_x|Ψ⟩||
# This MUST be real and positive since J_x is Hermitian
C0 = inner(Ψ', Jx, Φ)
push!(times, 0.0)
push!(Ct_values, C0)
push!(Ct_abs, abs(C0))
println("Initial correlation C(0) = ⟨Ψ|Jx|Φ⟩: ", C0)
println(" Re[C(0)] = ", real(C0), " (should be positive)")
println(" Im[C(0)] = ", imag(C0), " (should be ~0 for Hermitian J)")
if real(C0) < 0
@warn "C(0) has negative real part - check operator definition!"
end
if abs(imag(C0)) > 1e-10 * abs(real(C0))
@warn "C(0) has significant imaginary part - check operator Hermiticity!"
end
# ========================= Time evolution loop =============================#
println("Starting TDVP time evolution ")
# Number of time steps
nsteps = Int(round(tmax / dt))
for step in 1:nsteps
t_current = step * dt
# Evolve |Φ⟩: d/dt = -i H |Φ⟩
global Φ = tdvp(H,-im*dt, Φ;tdvp_kwargs...)
# Evolve |Ψ⟩: d/dt = -i H |Ψ⟩
global Ψ = tdvp(H,-im*dt, Ψ;tdvp_kwargs...)
# Compute correlation C(t) = ⟨Ψ(t)|Jx|Φ(t)⟩
C_t = inner(Ψ', Jx, Φ)
push!(times, t_current)
push!(Ct_values, C_t)
push!(Ct_abs, abs(C_t))
@printf("t = %.2f | Re[C(t)] = %12.6e | |C(t)| = %12.6e | Max χ = %d\n", t_current, real(C_t), abs(C_t), maxlinkdim(Φ))
end
# ==========================================================================#
# ========================= Compute Integral (DC Response) =================#
# We compute the running integral of Re[C(t)] using the trapezoidal rule.
# The long-time saturation value relates to the DC conductivity/diffusion.
re_C = real.(Ct_values)
running_integral = zeros(Float64, length(times))
# Trapezoidal numerical integration with damping factor exp(-η*t) for Kubo formula
# This regularizes the integral and gives the DC conductivity response
for i in 2:length(times)
# Include exponential damping: ∫ e^{-ηt} Re[C(t)] dt
damped_C_i = re_C[i] * exp(-η * times[i])
damped_C_im1 = re_C[i-1] * exp(-η * times[i-1])
# Area of trapezoid = base * average height
running_integral[i] = running_integral[i-1] + dt * 0.5 * (damped_C_i + damped_C_im1)
end
final_integral = running_integral[end]
println("\n=======================================================")
println("Kubo formula integral with damping η = ", η)
@printf("∫₀^%.1f e^{-ηt} Re[C(t)] dt = %.6e\n", tmax, final_integral)
println("This quantity is related to DC spin conductivity.")
println("For physical systems, this should be positive.")
if final_integral < 0
@warn "Final integral is negative - may indicate numerical issues or insufficient evolution time."
end
println("=======================================================")
# Plot 1: Real and Imaginary parts
# ==========================================================================#
# ========================= Visualization & Export =========================#
# Create the plots
# ==========================================================================#
# ========================= Visualization & Export =========================#
# Create the plots
p1 = plot(times, real.(Ct_values), label="Re[C(t)]", lw=2, color=:blue)
plot!(p1, times, imag.(Ct_values), label="Im[C(t)]", lw=2, color=:red)
ylabel!(p1, "Correlation C(t)")
title!(p1, "Spin Current Correlation")
p2 = plot(times, Ct_abs, label="|C(t)|", lw=2, color=:black, fill=(0, 0.2, :black))
ylabel!(p2, "|C(t)|")
# Plot 3: The Running Integral
p3 = plot(times, running_integral, label="∫ Re[C(t')] dt'", lw=2, color=:forestgreen)
ylabel!(p3, "Integral")
xlabel!(p3, "Time (t)")
title!(p3, "Running Integral (Conductivity)")
# Combine and SAVE (Updated layout to 3 rows)
full_plot = plot(p1, p2, p3, layout=(3,1), size=(800, 1000))
savefig(full_plot, "correlation_results.png")
println("Plot saved as 'correlation_results.png'")
# ========================= Save Raw Data ==================================#
# Updated to include the running integral in the CSV
open("correlation_data.csv", "w") do io
println(io, "time,re_ct,im_ct,abs_ct,integral_re_ct") # Header
for i in 1:length(times)
@printf(io, "%.6f,%.6e,%.6e,%.6e,%.6e\n",
times[i], real(Ct_values[i]), imag(Ct_values[i]), Ct_abs[i], running_integral[i])
end
end
println("Data saved as 'correlation_data.csv'")
# Keep the window open (only needed if running from terminal)
# gui()
# readline() # This pauses the script so the window stays open
Can anyone help? Am I applying correctly the TDVP algorithm?