前回の記事で$S=1/2$のHeisenbergモデルの行列を構成して固有値を求める話を書きました。
今回はビット演算子を用いて、直接行列を構成しないで、ある状態$\ket{\psi}$に対して$H\ket{\psi}$を計算する方法を紹介します。

ビットシフト演算について

まずはビットシフト演算子について説明します。
例えば10進法で5は2進法では$101_{(2)}$となります。例えば各桁を右に1だけずらすと$10_{(2)}$が得られます。これをビットシフトと呼びます（本当は論理シフトとか算術シフトとか言いますが、ここでは深入りしません）。この計算をJuliaでは次のように行うことができます。

a = 5
@show b = a >> 1
#b = a >> 1 = 2

AND演算に関しては&を使うことでできます。
例えば5を2進法で表示したときに、2の位が1なのかどうかを判別するには以下のようにします。

a = 5
b = a >> 1 & 1

当然この結果としては0です。

次にビット反転について説明しましょう。ビットを反転させる（0を1に、あるいは1を0に）するには、排他的論理和⊻を使います。
同じ例を用いて、$5=101_{(2)}$の2の位の0を1に変えてみましょう。⊻は\veebarで入力できます。
1 << 1は先ほどのビットシフトの逆で左に1桁だけずらします。

a = 5
@show c = a ⊻ (1 << 1)
#c = a ⊻ 1 << 1 = 7

さて、以上の道具を使ってスピン系の計算をしてみましょう。

Heisenbergモデルのハミルトニアンと基底

$S=1/2$のHeisenbergモデルのハミルトニアンは以下のように書けるのでした。

H = J\sum_{i=1}^{N-1}\vec{S}_i\cdot\vec{S}_{i+1} = J\sum_{i=1}^{N-1}\left(\frac12 S_i^+S_{i+1}^-+\frac12 S_i^-S_{i+1}^+ +S_i^zS_{i+1}^z\right)

以下では例として$N＝4$を考えます。
そこで、状態$\ket{\psi}$を次のようなベクトル表示にします。

\begin{align}
\ket{\psi} &= x_0 \ket{\uparrow\uparrow\uparrow\uparrow}+x_1\ket{\uparrow\uparrow\uparrow\downarrow}+\cdots \\
&=\pmatrix{
x_0 \\ x_1 \\ \vdots \\ x_{15}}
\end{align}

ここでのアイデアはハミルトニアン行列を構成することではなく、直接$H\ket{\psi}$を計算することです。
これらのそれぞれの基底を10進法では0,1,2,...,15に対応させましょう。ビット列で表現すれば、

\begin{align}
\ket{\uparrow\uparrow\uparrow\uparrow} \to \ket{0000}\\
\ket{\uparrow\uparrow\uparrow\downarrow} \to \ket{0001}
\end{align}

のように対応させています。ここで、基底とサイトの関係については$\ket{\sigma_N\sigma_{N-1}\cdots\sigma_1}$とします（前回の記事や、多くの他の文献の場合には逆順になっていることが多いことに注意してください）。

例えば$S_1^zS_2^z\ket{0101}$を計算したければ、

new_psi = zeros(Float32, Int(2^N))
state = 5
s1 = state >> 0 & 1
s2 = state >> 1 & 1
new_psi[state+1] += 0.25 * ((-1)^s1) * (-1)^s2

とします。これで計算結果のベクトルにすることができます。Juliaは1-indexであることに注意して、state+1としています。
ではスピンフリップの項について、$(S_1^+S_2^-+S_1^-S_2^+)\ket{0101}$の計算は次のようにする。

if s1 != s2
    new_state = state ⊻ (1 << 0) ⊻ (1 << 1)
    new_psi[new_state+1] += 1.0
end

以上の計算をすべてのサイト間で行えば$H\ket{\psi}$を計算できるようになります。
磁場の項$S_1^z\ket{0101}$なども計算することが可能です。

new_psi[state+1] += 0.5 * ((-1)^s1)

で十分です。

全体のコード

全体のコードをまとめると次のようになります。

function hpsi(psi::Vector{Float32}, N::Int, J::Float32)
    """
    Heisenbergモデルのハミルトニアンとベクトルの積 H|ψ⟩ を計算
    psi: 入力状態ベクトル
    N: サイト数
    J: 交換相互作用定数
    """
    new_psi = zeros(Float32, length(psi))  # 出力状態ベクトル

    for state in 0:length(psi)-1
        x = psi[state+1]
        if x == 0.0
            continue  #係数が0ならスキップ
        end
        # スピン間の相互作用を計算
        for i in 0:(N - 2)
            j = i + 1  # 隣接サイト

            # 現在のスピン配置をビット列で解析
            spin_i = (state >> i) & 1
            spin_j = (state >> j) & 1

            # SzSz の寄与
            sz_energy = J * 0.25 * ((-1)^spin_i) * ((-1)^spin_j)
            new_psi[state+1] += sz_energy * x

            # SxSx + SySy の寄与 (スピンフリップ)
            if spin_i != spin_j
                # スピンフリップした新しい状態を生成
                new_state = state ⊻ (1 << i) ⊻ (1 << j)  # ビットフリップ
                new_psi[new_state+1] += J * 0.5 * x
            end
        end
    end
    return new_psi
end

試しに$H\ket{0101}$を計算してみましょう。

H\ket{0101} = -0.75\ket{0101} + 0.5\ket{0110}+0.5\ket{0011}+0.5\ket{1001}

になるはずです。

N = 4
J = Float32(1.0)
psi = zeros(Float32, 2^N)
psi[6] = 1.0

new_psi = hpsi(psi, N, J)
for i in 1:2^N
    if new_psi[i] != 0.0
        println("$(i-1) $(new_psi[i])")
    end
end
#3 0.5
#5 -0.75
#6 0.5
#9 0.5

このように、正しく計算できていることがわかります。

（追記）冪乗法で基底状態を求める

せっかくなので、冪乗法で基底状態を求めることもやってみましょう。

\begin{align}
{\bf y}_k&= H{\bf x}_k\\
{\bf x}_{k+1} &= {\bf y}_k/||{\bf y}_k||
\end{align}

というふうに定めます。このとき、絶対値が最も大きい固有値は

\lambda_0 = \lim_{k\to \infty} {\bf y}_k^\dagger{\bf x}_k

で与えられます。
基底状態が欲しいので、単位行列の作用分だけ加えます。具体的にはhpsiの中身で

#new_psi = zeros(Float32, length(psi))  # 出力状態ベクトル
new_psi = -1 * psi

とします。
ハミルトニアン$H$を繰り返しかける関数を次のように定義しましょう。この関数では${\bf x}_k$を入力として、
${\bf y}_k$,

${\bf x}_{k+1}$

を出力にしています。

function hkpsi(psi::Vector{Float32}, N::Int, J::Float32)
        nextpsi = hpsi(psi, N, J)
        norm = sqrt(nextpsi'*nextpsi)
        nextpsi_normed = nextpsi/norm
    return nextpsi, nextpsi_normed
end

例えば$N=10$として、基底状態のエネルギーを求めるには次のようにします。
ついでに前回の記事で作ったmake_HamをArpackを用いて固有値を求める計算と比較してみます。

function main()
    N = 10
    J = Float32(1.0)
    psi = zeros(Float32, 2^N)
    ind =2+2^3+2^5+2^7+2^9
    psi[ind] = 1.0
    psin, psin_normed= hkpsi(psi, N, J)
    for i in 1:200
        psi = psin_normed
        psin, psin_normed = hkpsi(psi, N, J)
    end
    return psin'*psi + 1
end
@time e = main()
println(e)

using SparseArrays
using Arpack
function hikaku()
    H = 0.25* make_Ham(10, 0.0)
    return eigs(H, nev=1, which=:SR)[1][1]
end

@time e = hikaku()
println(e)

#0.003054 seconds (1.21 k allocations: 1.599 MiB)
#-4.2580347
#0.485044 seconds (6.45 M allocations: 287.661 MiB, 4.02% gc time, 79.80% compilation time)
#-4.258035207282885

Arpackの方は疎行列でLanzcos法を用いて計算しているはずなので、時間がかかるのも仕方ないですが、それでも2桁の差が出ています。
ちなみに、$N=16$としてみると、疎行列を構成して計算する手法では全くダメになってしまいます。計算時間がだいぶ変わりますね。

  0.338363 seconds (1.21 k allocations: 100.775 MiB, 8.75% gc time)
-6.9116955
384.566716 seconds (22.37 G allocations: 1.139 TiB, 10.97% gc time, 0.11% compilation time)
-6.911737145575105

まとめ

基底状態のエネルギーを求めるには行列を作らない方がよい。