Juliaで量子スピン系のハミルトニアンを構成する

はじめに

次のようなハミルトニアン

H = \sum_{i,j}\left( 
J_x\sigma^x_i\sigma^x_{j} 
+J_y\sigma^y_i\sigma^y_{j}
+J_z\sigma^z_i\sigma^z_{j}
\right)

をJuliaで具体的に構成する方法について考えます。ここで、$ \sigma^x, \sigma^y , \sigma^z$ はパウリ行列を表します。また下付き添字は、スピンが位置するサイトの番号を表します。$N$を全サイト数とすると、 $\sigma^z_i$ は

\sigma^z_i = 
\sigma^0 \otimes \cdots \otimes \sigma^0 \otimes \sigma^z
\otimes \sigma^0 \otimes \cdots \otimes \sigma^0

で与えられる $2^N \times 2^N$ 行列です。ここで、$\otimes$ はクロネッカー積、$\sigma^0$ は $2\times 2$ の単位行列です。

パウリ行列の定義

パウリ行列とは、

\sigma^x = 
\begin{pmatrix}
0 & 1 \\ 1 & 0
\end{pmatrix}, \quad
\sigma^y = 
\begin{pmatrix}
0 & -i \\ i & 0
\end{pmatrix}, \quad
\sigma^z = 
\begin{pmatrix}
1 & 0 \\ 0 & -1
\end{pmatrix}

というものでした。具体的に成分が与えられているので、これをプログラムすることは難しくありません。

using SparseArrays

σ⁰ = sparse([1.0+0.0im 0.0; 0.0 1.0+0.0im])
σˣ = sparse([0.0 1.0+0.0im; 1.0+0.0im 0.0])
σʸ = sparse([0.0 -1.0im; 1.0im 0.0])
σᶻ = sparse([1.0+0.0im 0.0; 0.0 -1.0+0.0im])

Juliaではimで虚数単位を表します。ここでは後のことを見越して、疎行列型としてパウリ行列を定義しました。Juliaでは変数名にUnicode文字を使用できるので、σᶻのようなギリシャ文字や上付き添字などが使えます。

念の為、代数関係が成り立っていることを確認しておきます。

julia> σᶻ*σᶻ == σ⁰
true

julia> σˣ*σʸ - σʸ*σˣ == 2.0im*σᶻ
true

クロネッカー積

行列$A,,B$のクロネッカー積とは

A \otimes B =
\begin{pmatrix}
a_{11}B & \dots & a_{1N}B \\
\vdots & \ddots & \vdots \\
a_{N1}B & \dots & a_{NN}B
\end{pmatrix}

で定義される演算です。$a_{ij}$ は $A$ の $i,j$ 成分を表します。Juliaではkron(A,B)によって計算することが可能です。あるいは、

⊗(A, B) = kron(A, B)

のように⊗を定義することで、A⊗Bとより直観的に書くこともできます。

クロネッカー積の基本的な性質を確認しておきます。

julia> A, B, C, D = rand(2, 2), rand(2, 2), rand(2, 2), rand(2, 2)
julia> A⊗(B+C) ≈ A⊗B + A⊗C
true
julia> (B+C)⊗A ≈ B⊗A + C⊗A
true
julia> (A⊗B)⊗C ≈ A⊗(B⊗C)
true
julia> (A⊗B)*(C⊗D) ≈ (A*C)⊗(B*D)
true

≈は右辺と左辺が機械精度の範囲で一致するかどうかを判定する演算子で、このようなテストをするときに重宝します。Juliaに対応したエディタ上で\approxと入力すると出ます。

ハミルトニアンの構成法1（ベタ打ち, 1次元）

以上の準備のもと、1次元Ising模型のハミルトニアン

H = \sum_{i=1}^N 
\sigma^z_i\sigma^z_{i+1}

を構成してみます。$N=3,4$ のときの具体形を手計算で求めておくと

\mathrm{diag}(3,-1,-1,-1,-1,-1,-1,3) \\
\mathrm{diag}(4,0,0,0,0,-4,0,0,0,0,-4,0,0,0,0,4)

となります。$\mathrm{diag}(\dots)$ は対角成分がその引数で与えられ、それ以外の成分が0であるような行列を表します。

さて、Julia上ではハミルトニアンを

H3 = σᶻ⊗σᶻ⊗σ⁰ + σ⁰⊗σᶻ⊗σᶻ + σᶻ⊗σ⁰⊗σᶻ
H4 = σᶻ⊗σᶻ⊗σ⁰⊗σ⁰ + σ⁰⊗σᶻ⊗σᶻ⊗σ⁰ + σ⁰⊗σ⁰⊗σᶻ⊗σᶻ + σᶻ⊗σ⁰⊗σ⁰⊗σᶻ

と書くことができます。もはや計算メモを書くのと全く変わりません。手計算の結果と一致することは、

julia> using LinearAlgebra
julia> H3 == Diagonal([3,-1,-1,-1,-1,-1,-1,3])
true
julia> H4 == Diagonal([4,0,0,0,0,-4,0,0,0,0,-4,0,0,0,0,4])
true

により確認できます。

ハミルトニアンの構成法2（一般のN, 1次元）

一般の $N$ に対応するため、クロネッカー積を再帰的に計算できるようにしておきます。

function ⊗(A,B,C...)
    A = kron(A, B)
    for Ci in C
        A = kron(A, Ci)
    end
    return A
end

C...は可変個の引数を取れることを表します。これは次のように動作します。

julia> ⊗(σᶻ,σᶻ,σ⁰,σ⁰) == σᶻ⊗σᶻ⊗σ⁰⊗σ⁰
true
julia> ⊗([σᶻ,σᶻ,σ⁰,σ⁰]...) == σᶻ⊗σᶻ⊗σ⁰⊗σ⁰
true

あとは、[σᶻ,σᶻ,σ⁰,σ⁰] に相当する配列を自由に渡せるようにすればOKです。例えば次が考えられます。

function set_spins(N, sites, σs)
    list_mats = fill(σ⁰, N)
    for (site, σ) in zip(sites, σs)
        list_mats[site] = σ
    end
    return list_mats
end

julia> set_spins(4, [1,2], [σᶻ,σᶻ]) == [σᶻ,σᶻ,σ⁰,σ⁰]
true
julia> set_spins(4, [1,4], [σᶻ,σˣ]) == [σᶻ,σ⁰,σ⁰,σˣ]
true

これらを用いると、Isingハミルトニアンは

function HIsing(N)
    H = ⊗(set_spins(N, [N,1], [σᶻ,σᶻ])...)
    for i in 1:N-1
        H += ⊗(set_spins(N, [i,i+1], [σᶻ,σᶻ])...)
    end
    return H
end

で構成できます。前節で作ったものと比較すると、

julia> HIsing(3) == H3
true
julia> HIsing(4) == H4
true

となります。一般のスピン系の場合も同様に構成することができます。冒頭のHeisenbergハミルトニアンは

function HHeisenberg(N)
    H = ⊗(set_spins(N, [N,1], [σˣ,σˣ])...)
    H += ⊗(set_spins(N, [N,1], [σʸ,σʸ])...)
    H += ⊗(set_spins(N, [N,1], [σᶻ,σᶻ])...)
    for i in 1:N-1
        H += ⊗(set_spins(N, [i,i+1], [σˣ,σˣ])...)
        H += ⊗(set_spins(N, [i,i+1], [σʸ,σʸ])...)
        H += ⊗(set_spins(N, [i,i+1], [σᶻ,σᶻ])...)
    end
    return H
end

で与えられます。

量子系である以上、取りうる $N$ の大きさにはタイトな限界がありますが、1次元ハイゼンベルク模型ならば、パソコン上でも $N=20$ 程度は工夫なしにメモリに乗せることができると思います。

2次元正方格子の場合

ここでは、$N \times N$の正方格子上のIsing模型を考えます。格子に以下のように番号を振ったとしましょう。

SiteIndeces_2d(N) = transpose(reshape(1:N^2, N,N))

julia> SiteIndeces_2d(3)
3×3 Transpose{Int64,Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}}}:
 1  2  3
 4  5  6
 7  8  9

正方格子上で周期条件のもとにIsing模型を考えると、$N=3$の場合、ボンドは(1-2)(2-3)(3-1)(4-5)(5-6)(6-4)(7-8)(8-9)(9-7)(1-4)(4-7)(7-1)(2-5)(5-8)(8-2)(3-6)(6-9)(9-3)の計18本あります。一般の$N$に対応するためには、これらのボンドを数え上げるコードを書く必要があります。落ち着いて考えれば難しくありませんが、ちょっと面倒だし、どこかでミスをしそうです。そこでJuliaのCartesianIndexを利用することを考えます。CartesianIndexは以下のように動作します。

A = [1 2 3; 4 5 6; 7 8 9]
A[CartesianIndex(1,1)] == 1
A[CartesianIndex(1,2)] == 2
A[CartesianIndex(1,3)] == 3
A[CartesianIndex(2,1)] == 4
A[CartesianIndex(2,2)] == 5
A[CartesianIndex(2,3)] == 6
A[CartesianIndex(3,1)] == 7
A[CartesianIndex(3,2)] == 8
A[CartesianIndex(3,3)] == 9

これに加えて、周期境界条件に対応するため、

function pbc(index, N)
    i = index % N
    i > 0 ? i : i+N
end

pbc.(-3:6, 3) == [3,1,2,3,1,2,3,1,2,3] # true

というものを用意しておきます。これらを用いて、

function bonds(N)
    A = SiteIndeces_2d(N)
    v = [(A[CartesianIndex(i,j)], A[CartesianIndex(i,pbc(j+1,N))]) for i in 1:N for j in 1:N]
    h = [(A[CartesianIndex(i,j)], A[CartesianIndex(pbc(i+1,N),j)]) for i in 1:N for j in 1:N]
    return vcat(v, h)
end

とすれば、すべての近接2サイトを結ぶボンドのリストが得られます。

julia> bonds(3)
18-element Array{Tuple{Int64,Int64},1}:
 (1, 2)
 (2, 3)
 (3, 1)
 (4, 5)
 (5, 6)
 (6, 4)
 (7, 8)
 (8, 9)
 (9, 7)
 (1, 4)
 (2, 5)
 (3, 6)
 (4, 7)
 (5, 8)
 (6, 9)
 (7, 1)
 (8, 2)
 (9, 3)

これで、2次元正方格子上のIsingハミルトニアンを構成する準備ができました。

function HIsing_2d(N)
    H = spzeros(2^(N*N), 2^(N*N))
    for b in bonds(N)
        H += ⊗(set_spins(N*N, [b[1],b[2]], [σᶻ,σᶻ])...)
    end
    return H
end

おまけ（もっと直感的に書く方法）

上のような方法でも相当直感的で可読性の高いコードが書けると思いますが、せっかく上付き/下付き添字を使ってコードが書けるのだから、σᶻ₂と書くだけで自動的にσ⁰⊗σᶻ⊗σ⁰...と解釈してくれないものだろうかという欲望が湧いてきます。同様のことが「Julia で下付き文字をインデックスに変換するマクロを書いてみた。」で考察されていたので、これを参考にマクロを使って所望の機能が実現できないか考えてみます。

σᶻ₂という記号が使いたいからといって、いきなり

julia> σᶻ₂
UndefVarError: σᶻ₂ not defined

と書いても、当然そんな記号は定義されていないと怒られてしまいます。しかし例えば、

julia> macro to_str(x) string(:($x)) end
julia> @to_str σᶻ₂
"σᶻ₂"

のようにマクロを使って評価することはできます。そこで、

function sub_to_num(sub)
    if sub == '₁'
        return 1
    elseif sub == '₂'
        return 2
        (省略)
    else
        UndefVarError()
    end
end
# e.g. sub_to_num('₁') == 1

function to_siteindex(array)
    n = length(array)
    sum([array[i]*10^(n-i) for i in 1:n])
end
# e.g. to_siteindex(sub_to_num.(['₁','₂'])) == 12

macro _spin(S, N)
    indeces = []
    spins = []
    for str in split(string(:($S)), "σ")[2:end]
        str = [s for s in "σ"*str]
        index = to_siteindex(sub_to_num.(str[3:end]))
        push!(indeces, index)
        @assert index <= N
        spin = eval(Meta.parse(prod(str[1:2])))
        push!(spins, spin)
    end
    ⊗(set_spins(N, indeces, spins)...)
end

のようなことをすれば@_spin σᶻ₂σᶻ₃ 4がσ⁰⊗σᶻ⊗σᶻ⊗σ⁰と等価になります。ただしこのままでは、2項以上の和を書くときに

(@_spin σᶻ₁σᶻ₂ 3) + (@_spin σᶻ₂σᶻ₃ 3) + (@_spin σᶻ₃σᶻ₁ 3)

などと書かねばならず煩わしいので、この式自身を

macro spin(S, N)
    str = "(@_spin " * replace(string(:($S)), "+"=>string(N)*")+(@_spin") * " " * string(N)* ")"
    Meta.parse(str)
end

というマクロで展開するようにすれば、

# N=3 Ising
H3 = @spin σᶻ₁σᶻ₂+σᶻ₂σᶻ₃+σᶻ₃σᶻ₁ 3

# N=4 Ising
H4 = @spin σᶻ₁σᶻ₂+σᶻ₂σᶻ₃+σᶻ₃σᶻ₄+σᶻ₄σᶻ₁ 4

でイジングハミルトニアンが書けます。

ということで、一応マクロにマクロを重ねれば、コードの見た目を整えることはいくらでも可能と言えそうです。ただ実質はspin("σᶻ₁σᶻ₂+σᶻ₂σᶻ₃+σᶻ₃σᶻ₁", 3)というような関数を定義して実行するのと変わらないので、

spin("σᶻ₁σᶻ₂+σᶻ₂σᶻ₃+σᶻ₃σᶻ₁", 3)
@spin σᶻ₁σᶻ₂+σᶻ₂σᶻ₃+σᶻ₃σᶻ₁ 3

のどちらを綺麗に感じるかという趣味の問題です。

まとめ

Juliaで量子スピン系のハミルトニアンを構成する方法を見てきました。サンプルコードはGistで公開しています。またここで定義した関数たちをパッケージとしてまとめたものをGitHubで公開しています。パッケージを自作したい人にとって参考になると思います。パッケージはブラウザ上でお試し使用することができます。（ただし、ビルドにちょっと時間がかかります。）