見出し画像

厳密対角化によるHaldane gap in Julia

研究テーマ探しから現実逃避をするためJuliaを学んでいます。スピン1 Haldane gapを計算するコードをChatGPTに作らせました。

以下コードなどはすべてJupyterに載せてあります。ダウンロードが嫌な人は続きをご覧ください。

using LinearAlgebra
using Arpack
using Plots

const ⊗ = kron

# Pauli 行列
function pauli()
    Sx = 1/sqrt(2)*[0 1 0; 1 0 1; 0 1 0]
    Sy = 1/sqrt(2)*[0 -im 0; im 0 -im; 0 im 0]
    Sz = [1 0 0; 0 0 0; 0 0 1]
    return Sx, Sy, Sz
end

# Heisenberg相互作用の局所項
function heisenbergint(Sx,Sy,Sz)
    return Sx ⊗ Sx + Sy ⊗ Sy + Sz ⊗ Sz 
end

# Hamiltonian
function hamiltonian(L::Int)
    Sx, Sy, Sz = pauli()
    H=zeros(ComplexF64, 3^L,3^L)
    for i in 1:(L-1)
        H += I(3^(i-1)) ⊗ heisenbergint(Sx,Sy,Sz) ⊗ I(3^(L-i-1))
    end
    return H
end

gaps = []
range=3:8
for L in range
    gap=eigs(hamiltonian(L),nev=2,which=:SR)[1]
    @time gaps = append!(gaps, gap[2]-gap[1])
end
@show gaps

L=6
@time eigenvalues=eigs(hamiltonian(L), nev=3^L,which=:SR)[1]
@show eigenvalues

# gaps
plot(range, real(gaps), yscale=:log10)
xlabel!("length")
ylabel!("Haldane gap")
title!("Haldane gaps of Spin-1 AF Heisenberg chain")

# energy spectrum
plt=plot(yshowaxis=false,xlim=(0,2))
xlabel!("energy")
title!("Spin-1 AF Heisenberg chain energy spectrum")

for j in 1:3^L-5
    plot!(fill(abs(eigenvalues[j]),2),[0,.5],label="", lc=:black)
end

display(plt)

いいなと思ったら応援しよう!