物理学における強束縛模型を解くパッケージQuadraticHamiltonians.jlを以前作りましたので、紹介します。
こちらは、固体物理のハミルトニアンを与えた時に、その平均場近似を行なって平均場を計算することができます。たとえば、超伝導ギャップ方程式を解いたりすることができます。

問題

ハミルトニアンとして、

\hat{H} \equiv \sum_{i j} \sum_{\alpha \beta} H_{ij} c_{i \alpha}^{\dagger} c_{j \beta}

を考えたり、

\hat{H}_{\rm BdG} \equiv \sum_{i j} \sum_{\alpha \beta} H_{ij} c_{i \alpha}^{\dagger} c_{j \beta} + \sum_{i j} \sum_{\alpha \beta} \bar{H}_{ij} c_{i \alpha} c_{j \beta}^{\dagger} + + \sum_{i j} \sum_{\alpha \beta} \Delta_{ij} c_{i \alpha}^{\dagger} c_{j \beta}^{\dagger} +  \sum_{i j} \sum_{\alpha \beta} \bar{\Delta}_{ij} c_{i \alpha} c_{j \beta}

を考えます。こちらは超伝導状態を記述するハミルトニアンです。

インストール

インストールは他のパッケージと同じで、

add QuadraticHamiltonians

です。

ハミルトニアンの作り方

4サイトのハミルトニアンは

julia> H = Hamiltonian(4)
---------------------------------
Hamiltonian: 
Num. of sites: 4
Num. of internal degree of freedom: 1
H = 
---------------------------------

とします。次に、演算子を定義します。

julia> c1 = FermionOP(1)
+C_{1,1}

julia> c2 = FermionOP(2)
+C_{2,1}

ここでは、フェルミオンの消滅演算子を定義しました。
これを使って、第二量子化表示でハミルトニアンに項を足すことができます。

julia> H += 2.0*c1'*c1 + 3*c2'*c2 + 2*c1'*c2 + 2*(c1'*c2)'
---------------------------------
Hamiltonian: 
Num. of sites: 4
Num. of internal degree of freedom: 1
H = +2.0C_{1,1}^+C_{1,1} +2.0C_{2,1}^+C_{1,1} +2.0C_{1,1}^+C_{2,1} +3.0C_{2,1}^+C_{2,1} 
---------------------------------

ハミルトニアンの行列要素は以下の二つの方法で定義できます。

H[2, 1] = 21
H[c3', c1] = 42

定義したハミルトニアンは行列として扱うことができ、ベクトルと掛け算ができます。

julia> x = rand(4)
4-element Vector{Float64}:
 0.18622544459032153
 0.36202722147302
 0.9038045102532898
 0.8497382867289013

julia> H*x
4-element Vector{Float64}:
 1.096505332126683
 1.4585325535997031
 0.0
 0.0

行列要素が欲しければ、

a = H[1,2]

したり、演算子で

c1 = FermionOP(1)
c2 = FermionOP(2)
a = H[c1',c2]

で指定したりできます。

もし、ハミルトニアン行列要素が複素数になり得る場合には、

julia> H = Hamiltonian(ComplexF64,4)
---------------------------------
Hamiltonian: 
Num. of sites: 4
Num. of internal degree of freedom: 1
H = 
---------------------------------

のようにすると、行列要素を複素数にしてもよくなります。
これを使うと、

julia> H += 2.0*c1'*c1 + 3*c2'*c2 + exp(0.2*im)*c1'*c2 + (exp(0.2*im)*c1'*c2)' 
---------------------------------
Hamiltonian: 
Num. of sites: 4
Num. of internal degree of freedom: 1
H = +2.0C_{1,1}^+C_{1,1} +(0.9800665778412416 - 0.19866933079506122im)C_{2,1}^+C_{1,1} +(0.9800665778412416 + 0.19866933079506122im)C_{1,1}^+C_{2,1} +3.0C_{2,1}^+C_{2,1} 
---------------------------------

のように、複素数をハミルトニアンに入れることができます。

多軌道模型

上の例ではサイトの数だけを指定していますが、内部自由度を取り扱うことができます。

julia> H = Hamiltonian(4,num_internal_degree=2)
---------------------------------
Hamiltonian: 
Num. of sites: 4
Num. of internal degree of freedom: 2
H = 
---------------------------------

たとえば、このようにすると、内部自由度が2になっており、これはスピンなどを扱えることを意味しています。演算子は

julia> c1up = FermionOP(1,1)
+C_{1,1}

julia> c1down = FermionOP(1,2)
+C_{1,2}

julia> c2down = FermionOP(2,2)
+C_{2,2}

などで定義できます。内部自由度を表すインデックスが追加されています。

超伝導

超伝導の平均場ハミルトニアンを使いたい場合には、

julia> H = Hamiltonian(4,isSC=true)
---------------------------------
Hamiltonian: 
Num. of sites: 4
Num. of internal degree of freedom: 1
Superconducting state
H = 
---------------------------------

のように、isSC=trueとします。
これにより、

julia> H += -1*(c1'*c1 - c1*c1') + 0.2*c2'c2' + (0.2*c2'*c2')'
---------------------------------
Hamiltonian: 
Num. of sites: 4
Num. of internal degree of freedom: 1
Superconducting state
H = -1.0C_{1,1}^+C_{1,1} +0.2C_{2,1}C_{2,1} +C_{1,1}C_{1,1}^+ +0.2C_{2,1}^+C_{2,1}^+ 
---------------------------------

のように超伝導平均場ハミルトニアンを定義することができ、
行列も

julia> H_sp = construct_matrix(H)
8×8 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 -1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅    ⋅    ⋅    ⋅    ⋅   0.2   ⋅    ⋅ 
   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅ 
   ⋅   0.2   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅

のようにできています。

強束縛模型の例

ハミルトニアンの例を示します。たとえば、

using QuadraticHamiltonians
function test2()
    N = 16
    ham = Hamiltonian(N)
    t = -1
    for i = 1:N
        ci = FermionOP(i)
        j = i + 1
        j += ifelse(j > N, -N, 0)
        cj = FermionOP(j)
        ham += t * ci' * cj

        j = i - 1
        j += ifelse(j < 1, N, 0)

        cj = FermionOP(j)
        ham += t * ci' * cj
    end
    display(ham)
    ham_matrix = construct_matrix(ham)
    display(ham_matrix)
end
test2()

とすると、出力は

---------------------------------
Hamiltonian: 
Num. of sites: 16
Num. of internal degree of freedom: 1
H = -1.0C_{1,1}^+C_{2,1} -1.0C_{1,1}^+C_{16,1} -1.0C_{2,1}^+C_{3,1} -1.0C_{2,1}^+C_{1,1} -1.0C_{3,1}^+C_{4,1} -1.0C_{3,1}^+C_{2,1} -1.0C_{4,1}^+C_{5,1} -1.0C_{4,1}^+C_{3,1} -1.0C_{5,1}^+C_{6,1} -1.0C_{5,1}^+C_{4,1} -1.0C_{6,1}^+C_{7,1} -1.0C_{6,1}^+C_{5,1} -1.0C_{7,1}^+C_{8,1} -1.0C_{7,1}^+C_{6,1} -1.0C_{8,1}^+C_{9,1} -1.0C_{8,1}^+C_{7,1} -1.0C_{9,1}^+C_{10,1} -1.0C_{9,1}^+C_{8,1} -1.0C_{10,1}^+C_{11,1} -1.0C_{10,1}^+C_{9,1} -1.0C_{11,1}^+C_{12,1} -1.0C_{11,1}^+C_{10,1} -1.0C_{12,1}^+C_{13,1} -1.0C_{12,1}^+C_{11,1} -1.0C_{13,1}^+C_{14,1} -1.0C_{13,1}^+C_{12,1} -1.0C_{14,1}^+C_{15,1} -1.0C_{14,1}^+C_{13,1} -1.0C_{15,1}^+C_{16,1} -1.0C_{15,1}^+C_{14,1} -1.0C_{16,1}^+C_{1,1} -1.0C_{16,1}^+C_{15,1} 
---------------------------------
16×16 SparseArrays.SparseMatrixCSC{Float64, Int64} with 32 stored entries:
⎡⠪⡢⡀⠀⠀⠀⠀⠈⎤
⎢⠀⠈⠪⡢⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠪⡢⡀⠀⎥
⎣⡀⠀⠀⠀⠀⠈⠪⡢⎦

となります。ハミルトニアンを定義する際に、生成演算子だけを取り回していればよいという意味で、わかりやすくなっているかと思います。

2次元s波超伝導

次に、超伝導の例です。

μ = -1.5
Nx = 128
Ny = 128
N = Nx * Ny
Δ = 0.5
ham = Hamiltonian(N, isSC=true)
t = -1
hops = [(+1, 0), (-1, 0), (0, +1), (0, -1)]
for ix = 1:Nx
    for iy = 1:Ny
        i = (iy - 1) * Nx + ix
        ci = FermionOP(i)
        for (dx, dy) in hops
            jx = ix + dx
            jx += ifelse(jx > Nx, -Nx, 0)
            jx += ifelse(jx < 1, Nx, 0)
            jy = iy + dy
            jy += ifelse(jy > Ny, -Ny, 0)
            jy += ifelse(jy < 1, Ny, 0)
            j = (jy - 1) * Nx + jx
            cj = FermionOP(j)
            ham += -1 * (ci' * cj - cj * ci')
        end
        ham += -μ * (ci' * ci - ci * ci')
        ham += Δ * ci' * ci' + Δ * ci * ci
    end
end

となります。

平均場の計算方法

以上のような方法で定義されたハミルトニアンを用いて、平気場の計算を行います。
平均場の計算は

T = 0.1
m = Meanfields_solver(ham, T)

のように簡単にsolverを定義します。ここで、hamは上で定義したハミルトニアンで、Tは温度です。計算手法は縮約シフト共役勾配法です。

もし、Chebyshev多項式展開法で計算したければ、

m = Meanfields_solver(ham, T, method="Chebyshev", nmax=200)

LK-Chebyshevの場合は、

m = Meanfields_solver(ham, T, method="Chebyshev", isLK=true, nmax=200)

とします。
たとえば、平均場$\langle c_1 c_1 \rangle$は

c1 = FermionOP(1)
Gij0 = calc_meanfields(m, c1, c1) #<c1 c1>

と計算できます。
ハミルトニアンをアップデートしたい場合は、

update_hamiltonian!(m, value, ci', ci')
update_hamiltonian!(m, value, ci, ci)

とします。使っているハミルトニアンを取り出したい場合には、

ham = get_hamiltonian(m)

とします。

グリーン関数や局所状態密度の計算

ハミルトニアンが生成消滅演算子の二種類の積の和で定義されている場合、グリーン関数は

\begin{align}
G_{ij}(z) \equiv \langle c_i^{\dagger} [z \hat{I} - \hat{H}]^{-1} c_j \rangle
\end{align}

で計算できます。
ハミルトニアン行列を使うと、

\begin{align}
G_{ij}(z) = \left[ (z \hat{I} - \hat{H})^{-1} \right]_{ij}
\end{align}

のように書くことができます。
局所状態密度は

\begin{align}
N_i(\omega) = -\frac{1}{\pi} {\rm Im} \: \lim_{\eta \rightarrow 0+} G_{ii}(\omega + i \eta) 
\end{align}

で計算できます。

このパッケージを使うと、グリーン関数は

z = 0.1 + 0.2im
i = 1
c1 = FermionOP(i)
Gij = calc_greenfunction(ham, z, c1', c1)

で計算できます。ここで、RSCG法を使っています。

RSCG法は複数の周波数を一気に計算できますので、

num = 200
zs = zeros(ComplexF64, num)
ene = range(-6, 6, length=num)
eta = 0.05
for i = 1:num
    zs[i] = ene[i] + im * eta
end
ix = Nx ÷ 2
iy = ix
i = (iy - 1) * Nx + ix
c1 = FermionOP(i)
Gij = calc_greenfunction(ham, zs, c1', c1)

とすることが可能です。

例

s波超伝導体

s波超伝導体の例を示します。

using QuadraticHamiltonians
using Plots
using BenchmarkTools

function main()
    μ = -1.5
    Nx = 128
    Ny = 128
    N = Nx * Ny
    Δ = 0.5
    ham = Hamiltonian(N, isSC=true)
    t = -1
    hops = [(+1, 0), (-1, 0), (0, +1), (0, -1)]
    for ix = 1:Nx
        for iy = 1:Ny
            i = (iy - 1) * Nx + ix
            ci = FermionOP(i)
            for (dx, dy) in hops
                jx = ix + dx
                jx += ifelse(jx > Nx, -Nx, 0)
                jx += ifelse(jx < 1, Nx, 0)
                jy = iy + dy
                jy += ifelse(jy > Ny, -Ny, 0)
                jy += ifelse(jy < 1, Ny, 0)
                j = (jy - 1) * Nx + jx
                cj = FermionOP(j)
                ham += -1 * (ci' * cj - cj * ci')
            end
            ham += -μ * (ci' * ci - ci * ci')
            ham += Δ * ci' * ci' + Δ * ci * ci
        end
    end

    num = 200
    zs = zeros(ComplexF64, num)
    ene = range(-6, 6, length=num)
    eta = 0.05
    for i = 1:num
        zs[i] = ene[i] + im * eta
    end
    ix = Nx ÷ 2
    iy = ix
    i = (iy - 1) * Nx + ix
    c1 = FermionOP(i)
    Gij = calc_greenfunction(ham, zs, c1', c1)
    plot(ene, -imag.(Gij) / π)
    savefig("ldos.png")

    T = 0.1
    m = Meanfields_solver(ham, T, method="Chebyshev", nmax=200)
    c1up = FermionOP(1)
    c1down = FermionOP(1)
    Gij0 = calc_meanfields(m, c1up, c1down) #<c1up c1down>
    @btime calc_meanfields($m, $c1up, $c1down) #<c1up c1down>
    println("Chebyshev: ", Gij0)
    m2 = Meanfields_solver(ham, T)
    Gij0 = calc_meanfields(m2, c1up, c1down) #<c1up c1down>
    @btime calc_meanfields($m2, $c1up, $c1down) #<c1up c1down>
    println("RSCG ", Gij0)

    m3 = Meanfields_solver(ham, T, method="Chebyshev", isLK=true, nmax=200)
    Gij0 = calc_meanfields(m3, c1up, c1down) #<c1up c1down>
    @btime calc_meanfields($m3, $c1up, $c1down) #<c1up c1down>
    println("Chebyshev with LKvectors: ", Gij0)

end

出力結果は、

The solver is the Chebyshev polynomial method
  27.004 ms (10 allocations: 1.00 MiB)
Chebyshev: -0.1761175169590693
The solver is the RSCG
num. of Matsubara freqs. 24
  33.507 ms (6412 allocations: 37.76 MiB)
RSCG -0.17610594510770422 - 4.9649942910811176e-15im
The solver is the Chebyshev polynomial method
  20.113 ms (42 allocations: 2.13 MiB)
Chebyshev with LKvectors: -0.17611751734430567

となります。そして、局所状態密度は下図のようになります。ちゃんと超伝導ギャップがあいていることがわかります。

自己無撞着計算

ギャップ方程式を自己無撞着に解くことで平均場を求めます。

using QuadraticHamiltonians

function make_hamiltonian(Nx, Ny, μ, Δs)
    N = Nx * Ny
    ham = Hamiltonian(N, isSC=true)
    t = -1
    hops = [(+1, 0), (-1, 0), (0, +1), (0, -1)]
    for ix = 1:Nx
        for iy = 1:Ny
            i = (iy - 1) * Nx + ix
            ci = FermionOP(i)
            for (dx, dy) in hops
                jx = ix + dx
                jx += ifelse(jx > Nx, -Nx, 0)
                jx += ifelse(jx < 1, Nx, 0)
                jy = iy + dy
                jy += ifelse(jy > Ny, -Ny, 0)
                jy += ifelse(jy < 1, Ny, 0)
                j = (jy - 1) * Nx + jx
                cj = FermionOP(j)
                ham += -1 * (ci' * cj - cj * ci')
            end
            ham += -μ * (ci' * ci - ci * ci')
            ham += Δs[i] * ci' * ci' + Δs[i] * ci * ci
        end
    end
    return ham
end

function update!(m, Δs)
    N, _ = size(m.hamiltonian)
    for i = 1:(N÷2)
        ci = FermionOP(i)
        update_hamiltonian!(m, Δs[i], ci', ci')
        update_hamiltonian!(m, Δs[i], ci, ci)
    end
end

function main()
    Nx = 24
    Ny = 24
    μ = -0.2
    Δ0 = 0.5
    Δs = Δ0 * ones(Nx * Ny)
    Δsnew = similar(Δs)
    T = 0.02
    U = -2

    ham = make_hamiltonian(Nx, Ny, μ, Δs)
    m = Meanfields_solver(ham, T)

    for it = 1:100
        @time Gs = map(i -> calc_meanfields(m, FermionOP(i), FermionOP(i)), 1:Nx*Ny)
        Δsnew .= real(U * Gs)
        res = sum(abs.(Δsnew .- Δs)) / sum(abs.(Δs))
        println("$(it)-th $(Δsnew[1]) $res")
        if res < 1e-4
            break
        end
        update!(m, Δsnew)
        Δs .= Δsnew
    end

end
main()

出力結果は、

The solver is the RSCG
num. of Matsubara freqs. 34
  1.357527 seconds (8.97 M allocations: 1.225 GiB, 7.59% gc time, 46.78% compilation time)
1-th 0.43392863318462976 0.13214274366572432
  0.827926 seconds (5.99 M allocations: 1.104 GiB, 5.64% gc time)
2-th 0.40181627897040467 0.07400376230129402
  0.881787 seconds (6.16 M allocations: 1.144 GiB, 5.09% gc time)
3-th 0.38493556921008637 0.04201103156997562
  0.884848 seconds (6.32 M allocations: 1.180 GiB, 5.40% gc time)
4-th 0.37568124638585454 0.024041167724526982
  0.891520 seconds (6.36 M allocations: 1.192 GiB, 5.27% gc time)
5-th 0.3704883477225273 0.013822612441196765
  0.889258 seconds (6.38 M allocations: 1.215 GiB, 4.45% gc time)
6-th 0.3675358557663468 0.007969231662928523
  0.904341 seconds (6.44 M allocations: 1.212 GiB, 5.19% gc time)
7-th 0.36584446964213757 0.004601915465106033
  0.946217 seconds (6.52 M allocations: 1.231 GiB, 5.29% gc time)
8-th 0.36487138389829876 0.0026599130537329862
  0.903315 seconds (6.49 M allocations: 1.224 GiB, 4.95% gc time)
9-th 0.36431012784136424 0.001538238218849528
  0.905315 seconds (6.52 M allocations: 1.232 GiB, 4.69% gc time)
10-th 0.3639859317169275 0.0008898496406517458
  0.904538 seconds (6.52 M allocations: 1.232 GiB, 4.80% gc time)
11-th 0.3637985531159851 0.0005148569471265923
  0.917017 seconds (6.56 M allocations: 1.240 GiB, 4.78% gc time)
12-th 0.36369015468817556 0.0002979221928351472
  0.945926 seconds (6.54 M allocations: 1.236 GiB, 5.33% gc time)
13-th 0.36362744231704996 0.00017239402356417755
  0.922917 seconds (6.50 M allocations: 1.225 GiB, 5.31% gc time)
14-th 0.3635911839621331 9.976804879798009e-5

となります。

スピン自由度がある場合

スピン自由度がある場合の例です。

using QuadraticHamiltonians
using Plots
using BenchmarkTools
function main2()
    μ = -1.5
    Nx = 128
    Ny = 128
    N = Nx * Ny
    Δ = 0.5
    ham = Hamiltonian(N, num_internal_degree=2, isSC=true)
    t = -1
    hops = [(+1, 0), (-1, 0), (0, +1), (0, -1)]
    for ix = 1:Nx
        for iy = 1:Ny
            i = (iy - 1) * Nx + ix
            for ispin = 1:2
                ci = FermionOP(i, ispin)
                jspin = ispin
                for (dx, dy) in hops
                    jx = ix + dx
                    jx += ifelse(jx > Nx, -Nx, 0)
                    jx += ifelse(jx < 1, Nx, 0)
                    jy = iy + dy
                    jy += ifelse(jy > Ny, -Ny, 0)
                    jy += ifelse(jy < 1, Ny, 0)
                    j = (jy - 1) * Nx + jx
                    cj = FermionOP(j, jspin)
                    ham += -1 * (ci' * cj - cj * ci')
                end
                ham += -μ * (ci' * ci - ci * ci')
            end
            ciup = FermionOP(i, 1)
            cidown = FermionOP(i, 2)
            ham += Δ * ciup' * cidown' + Δ * cidown * ciup
            ham += -Δ * cidown' * ciup' - Δ * ciup * cidown
        end
    end

    num = 200
    zs = zeros(ComplexF64, num)
    ene = range(-6, 6, length=num)
    eta = 0.05
    for i = 1:num
        zs[i] = ene[i] + im * eta
    end
    ix = Nx ÷ 2
    iy = ix
    i = (iy - 1) * Nx + ix
    c1 = FermionOP(i, 1)
    Gij = calc_greenfunction(ham, zs, c1', c1)
    plot(ene, -imag.(Gij) / π)
    savefig("ldos_spin.png")

    T = 0.1
    m = Meanfields_solver(ham, T, method="Chebyshev", nmax=200)
    c1up = FermionOP(1, 1)
    c1down = FermionOP(1, 2)
    Gij0 = calc_meanfields(m, c1up, c1down) #<c1up c1down>
    @btime calc_meanfields($m, $c1up, $c1down) #<c1up c1down>
    println("Chebyshev: ", Gij0)

    m2 = Meanfields_solver(ham, T)
    Gij0 = calc_meanfields(m2, c1up, c1down) #<c1up c1down>
    @btime calc_meanfields($m2, $c1up, $c1down) #<c1up c1down>
    println("RSCG ", Gij0)

    m3 = Meanfields_solver(ham, T, method="Chebyshev", isLK=true, nmax=200)
    Gij0 = calc_meanfields(m3, c1up, c1down) #<c1up c1down>
    @btime calc_meanfields($m3, $c1up, $c1down) #<c1up c1down>
    println("Chebyshev with LKvectors: ", Gij0)

end
main2()

出力は、

The solver is the Chebyshev polynomial method
  49.120 ms (10 allocations: 2.00 MiB)
Chebyshev: 0.1761175169590693
The solver is the RSCG
num. of Matsubara freqs. 24
  63.067 ms (6412 allocations: 75.01 MiB)
RSCG 0.17610594510770403 + 4.955324870183428e-15im
The solver is the Chebyshev polynomial method
  21.877 ms (42 allocations: 4.25 MiB)
Chebyshev with LKvectors: 0.1761175173443056

です。

トポロジカルs波超伝導体

Zeeman磁場とRashbaスピン軌道相互作用があるs波超伝導体はトポロジカル超伝導になることが知られています。それをやってみます。

using QuadraticHamiltonians

function make_TSC_hamiltonian(Nx, Ny, μ, Δs, h, α, isPBC)
    N = Nx * Ny
    ham = Hamiltonian(ComplexF64, N, num_internal_degree=2, isSC=true)
    hops = [(+1, 0), (-1, 0), (0, +1), (0, -1)]
    for ix = 1:Nx
        for iy = 1:Ny
            i = (iy - 1) * Nx + ix
            for ispin = 1:2
                ci = FermionOP(i, ispin)
                σy = ifelse(ispin == 1, -im, im)
                σx = 1
                σz = ifelse(ispin == 1, 1, -1)

                for (dx, dy) in hops
                    jx = ix + dx
                    if isPBC
                        jx += ifelse(jx > Nx, -Nx, 0)
                        jx += ifelse(jx < 1, Nx, 0)
                    end
                    jy = iy + dy
                    if isPBC
                        jy += ifelse(jy > Ny, -Ny, 0)
                        jy += ifelse(jy < 1, Ny, 0)
                    end

                    if 1 <= jx <= Nx && 1 <= jy <= Ny
                        j = (jy - 1) * Nx + jx
                        jspin = ispin
                        cj = FermionOP(j, jspin)
                        ham += -1 * (ci' * cj - cj * ci')

                        jspin = ifelse(ispin == 1, 2, 1)
                        cj = FermionOP(j, jspin)

                        if dy == 0
                            if dx == 1
                                ham += (α / (2im)) * (ci' * cj - cj * ci') * σy
                            else
                                ham += -1 * (α / (2im)) * (ci' * cj - cj * ci') * σy
                            end
                        elseif dx == 0
                            if dy == 1
                                ham += (α / (2im)) * (ci' * cj - cj * ci') * σx
                            else
                                ham += -1 * (α / (2im)) * (ci' * cj - cj * ci') * σx
                            end
                        end
                    end
                end
                ham += (-μ - h * σz) * (ci' * ci - ci * ci')
            end
            ciup = FermionOP(i, 1)
            cidown = FermionOP(i, 2)
            ham += Δs[i] * ciup' * cidown' + Δs[i] * cidown * ciup
            ham += -Δs[i] * cidown' * ciup' - Δs[i] * ciup * cidown
        end
    end


    return ham
end

function update!(m, Δs)
    N, _ = size(m.hamiltonian)
    for i = 1:(N÷4)
        ciup = FermionOP(i, 1)
        cidown = FermionOP(i, 2)
        update_hamiltonian!(m, -Δs[i], ciup', cidown')
        update_hamiltonian!(m, Δs[i], cidown', ciup')
        update_hamiltonian!(m, Δs[i], ciup, cidown)
        update_hamiltonian!(m, -Δs[i], cidown, ciup)
    end
end



function main()
    Nx = 16
    Ny = 16
    μ = 3.5
    Δ0 = 3
    Δs = Δ0 * ones(Nx * Ny)
    Δsnew = similar(Δs)
    T = 0.01
    U = -5.6
    h = 1
    α = 1

    isPBC = false
    ham = make_TSC_hamiltonian(Nx, Ny, μ, Δs, h, α, isPBC)
    m = Meanfields_solver(ham, T)

    for it = 1:100
        @time Gs = map(i -> calc_meanfields(m, FermionOP(i, 1), FermionOP(i, 2)), 1:Nx*Ny)
        Δsnew .= real(U * Gs)
        res = sum(abs.(Δsnew .- Δs)) / sum(abs.(Δs))
        println("$(it)-th $(Δsnew[1]) $res")
        if res < 1e-4
            break
        end
        update!(m, Δsnew)
        Δs .= Δsnew
    end

end
main()

出力結果は、

The solver is the RSCG
num. of Matsubara freqs. 40
  1.079232 seconds (6.06 M allocations: 707.212 MiB, 16.58% gc time, 64.22% compilation time)
1-th -1.873023776160689 1.6271730164203106
  0.601234 seconds (4.04 M allocations: 846.823 MiB, 4.71% gc time)
2-th -1.4418945612222904 0.20965023882070202
  0.918474 seconds (5.74 M allocations: 1.242 GiB, 4.02% gc time)
3-th -1.2113675935614676 0.13563888454719064
  1.213072 seconds (7.08 M allocations: 1.605 GiB, 4.38% gc time)
4-th -1.0683848898677126 0.09725266043470575
  1.519615 seconds (8.41 M allocations: 1.988 GiB, 4.20% gc time)
5-th -0.9719173169890015 0.07520905143665893
  1.774157 seconds (9.26 M allocations: 2.261 GiB, 4.52% gc time)
6-th -0.9031937676953646 0.0614030514203067
  1.960240 seconds (9.93 M allocations: 2.480 GiB, 4.87% gc time)
7-th -0.8522993398753158 0.052088904542545035
  2.113638 seconds (10.48 M allocations: 2.672 GiB, 4.82% gc time)
8-th -0.8134640044396014 0.04536935661275386
  2.238472 seconds (10.90 M allocations: 2.814 GiB, 5.02% gc time)
9-th -0.7830935035525877 0.040207500802570885
  2.127559 seconds (10.55 M allocations: 2.707 GiB, 4.57% gc time)
10-th -0.7588369948792518 0.03600999926600642
  2.334611 seconds (11.38 M allocations: 2.977 GiB, 4.48% gc time)
11-th -0.7391008228944542 0.032434049930203365
  2.442818 seconds (11.64 M allocations: 3.064 GiB, 4.63% gc time)
12-th -0.7227749038266035 0.02928518035552862
  2.496474 seconds (11.85 M allocations: 3.135 GiB, 4.41% gc time)
13-th -0.709069155872567 0.026456030788754803
  2.483085 seconds (11.97 M allocations: 3.181 GiB, 3.93% gc time)
14-th -0.6974110355699145 0.02388799764731334
  2.495291 seconds (12.05 M allocations: 3.216 GiB, 3.89% gc time)
15-th -0.6873791072315169 0.0215479761914521
  2.542803 seconds (12.16 M allocations: 3.257 GiB, 4.01% gc time)
16-th -0.6786586888786034 0.019415153943493132
  2.587446 seconds (12.25 M allocations: 3.295 GiB, 4.06% gc time)
17-th -0.671011679282689 0.017474092388203903
  2.626588 seconds (12.35 M allocations: 3.335 GiB, 4.09% gc time)
18-th -0.6642555683053231 0.01571143747168237
  2.681020 seconds (12.46 M allocations: 3.373 GiB, 4.22% gc time)
19-th -0.6582486676185052 0.014114520535718768
  2.712846 seconds (12.44 M allocations: 3.368 GiB, 4.94% gc time)
20-th -0.6528795616656254 0.012670927816306929
  2.687108 seconds (12.47 M allocations: 3.385 GiB, 4.42% gc time)
21-th -0.648059496125448 0.01136845977752193
  2.698632 seconds (12.62 M allocations: 3.432 GiB, 4.16% gc time)
22-th -0.643716834307313 0.010195258886917814
  2.774037 seconds (12.60 M allocations: 3.431 GiB, 4.52% gc time)
23-th -0.6397929998696529 0.00913994023009275
  2.804256 seconds (12.72 M allocations: 3.478 GiB, 4.61% gc time)
24-th -0.6362394570260124 0.00819171955905381
  2.775846 seconds (12.89 M allocations: 3.538 GiB, 4.08% gc time)
25-th -0.63301548265852 0.007340492709279601
  2.812354 seconds (12.98 M allocations: 3.568 GiB, 4.28% gc time)
26-th -0.6300864856968793 0.006576883266048058
  2.825387 seconds (13.03 M allocations: 3.593 GiB, 3.97% gc time)
27-th -0.6274227436671878 0.005892252532124953
  2.830010 seconds (13.08 M allocations: 3.613 GiB, 3.93% gc time)
28-th -0.6249984477604995 0.005278692772073217
  2.877537 seconds (13.10 M allocations: 3.617 GiB, 4.20% gc time)
29-th -0.6227909814572898 0.004729001813477821
  2.814024 seconds (13.14 M allocations: 3.633 GiB, 3.63% gc time)
30-th -0.6207803405763253 0.004236642241088458
  2.844506 seconds (13.15 M allocations: 3.640 GiB, 3.96% gc time)
31-th -0.6189487129414846 0.0037957051360124226
  3.005532 seconds (13.16 M allocations: 3.642 GiB, 4.95% gc time)
32-th -0.6172801374355674 0.003400858384159586
  2.982854 seconds (13.17 M allocations: 3.646 GiB, 4.76% gc time)
33-th -0.6157602345110169 0.0030473007844825273
  2.929872 seconds (13.14 M allocations: 3.638 GiB, 4.32% gc time)
34-th -0.6143759966380645 0.0027307205386831934
  2.833398 seconds (13.16 M allocations: 3.644 GiB, 3.74% gc time)
35-th -0.6131156094137203 0.0024472404597514195
  2.948779 seconds (13.18 M allocations: 3.653 GiB, 4.51% gc time)
36-th -0.6119683239104643 0.002193392792552506
  2.935843 seconds (13.21 M allocations: 3.664 GiB, 4.45% gc time)
37-th -0.6109243302390619 0.001966063888845468
  2.944080 seconds (13.22 M allocations: 3.668 GiB, 4.50% gc time)
38-th -0.6099746645043274 0.0017624669009970566
  3.032381 seconds (13.22 M allocations: 3.668 GiB, 4.98% gc time)
39-th -0.6091111276751108 0.001580108345611027
  2.912560 seconds (13.27 M allocations: 3.687 GiB, 3.91% gc time)
40-th -0.6083262105712027 0.001416755095863366
  3.070945 seconds (13.28 M allocations: 3.694 GiB, 4.85% gc time)
41-th -0.6076130419061138 0.0012704129923773262
  2.959192 seconds (13.28 M allocations: 3.693 GiB, 4.50% gc time)
42-th -0.6069653199286934 0.00113929572631485
  3.039093 seconds (13.31 M allocations: 3.700 GiB, 4.66% gc time)
43-th -0.6063772737453785 0.0010218056690919832
  3.049923 seconds (13.32 M allocations: 3.704 GiB, 4.75% gc time)
44-th -0.6058436172654165 0.0009165157705025486
  3.028901 seconds (13.31 M allocations: 3.704 GiB, 4.65% gc time)
45-th -0.6053595096728797 0.0008221475370220531
  3.014559 seconds (13.32 M allocations: 3.707 GiB, 4.55% gc time)
46-th -0.6049205193354682 0.0007375584043463334
  3.013027 seconds (13.32 M allocations: 3.709 GiB, 4.51% gc time)
47-th -0.604522590849413 0.0006617273730522861
  3.065978 seconds (13.34 M allocations: 3.713 GiB, 4.85% gc time)
48-th -0.6041620163787921 0.000593740534829054
  2.998548 seconds (13.33 M allocations: 3.712 GiB, 4.67% gc time)
49-th -0.6038354056749607 0.0005327789447378642
  3.008589 seconds (13.34 M allocations: 3.715 GiB, 4.53% gc time)
50-th -0.6035396613195223 0.00047811220217767336
  2.984136 seconds (13.34 M allocations: 3.717 GiB, 4.27% gc time)
51-th -0.603271956206971 0.0004290844742470022
  2.986102 seconds (13.34 M allocations: 3.717 GiB, 4.31% gc time)
52-th -0.6030297098728377 0.0003851109026020435
  3.008851 seconds (13.35 M allocations: 3.718 GiB, 4.27% gc time)
53-th -0.6028105687238502 0.0003456655189107658
  2.926974 seconds (13.35 M allocations: 3.720 GiB, 3.93% gc time)
54-th -0.6026123890946519 0.0003102796361298741
  3.013989 seconds (13.36 M allocations: 3.721 GiB, 4.47% gc time)
55-th -0.6024332167560249 0.00027853279022410063
  2.925957 seconds (13.35 M allocations: 3.719 GiB, 4.07% gc time)
56-th -0.6022712741026642 0.00025004790660750666
  2.990428 seconds (13.35 M allocations: 3.722 GiB, 4.25% gc time)
57-th -0.6021249422723195 0.00022448825972236396
  2.993714 seconds (13.35 M allocations: 3.719 GiB, 4.50% gc time)
58-th -0.6019927502404482 0.00020155103761652473
  2.972320 seconds (13.35 M allocations: 3.719 GiB, 4.31% gc time)
59-th -0.6018733605571956 0.0001809666539966145
  2.998017 seconds (13.36 M allocations: 3.723 GiB, 4.52% gc time)
60-th -0.6017655584431989 0.00016249232065718839
  2.958357 seconds (13.35 M allocations: 3.718 GiB, 4.23% gc time)
61-th -0.6016682408866258 0.00014591021496670423
  2.954279 seconds (13.36 M allocations: 3.721 GiB, 4.15% gc time)
62-th -0.6015804069224595 0.0001310259234340552
  2.987971 seconds (13.36 M allocations: 3.723 GiB, 4.38% gc time)
63-th -0.6015011472832793 0.00011766420738475466
  2.967308 seconds (13.35 M allocations: 3.722 GiB, 4.31% gc time)
64-th -0.6014296388255697 0.00010566954846444353
  2.952684 seconds (13.35 M allocations: 3.721 GiB, 4.15% gc time)
65-th -0.6013651365577349 9.490118655278245e-5

となります。

並列計算

並列計算も可能です。その場合には、

using Distributed
@everywhere using QuadraticHamiltonians

とDistributedパッケージを使います。
上で述べたトポロジカルs波超伝導の場合には、

function main()
    Nx = 16
    Ny = 16
    μ = 3.5
    Δ0 = 3
    Δs = Δ0 * ones(Nx * Ny)
    Δsnew = similar(Δs)
    T = 0.01
    U = -5.6
    h = 1
    α = 1
    isPBC = false

    ham = make_TSC_hamiltonian(Nx, Ny, μ, Δs, h, α, isPBC)
    m = Meanfields_solver(ham, T)

    for it = 1:100
        @time Gs = pmap(i -> calc_meanfields(m, FermionOP(i, 1), FermionOP(i, 2)), 1:Nx*Ny)
        Δsnew .= real(U * Gs)
        res = sum(abs.(Δsnew .- Δs)) / sum(abs.(Δs))
        println("$(it)-th $(Δsnew[1]) $res")
        if res < 1e-4
            break
        end
        update!(m, Δsnew)
        Δs .= Δsnew
    end

end
main()

とします。pmapが並列処理を行なっています。平均場の計算は各サイトで独立にできるために、その部分を並列化しています。

出力結果は、

julia -p 4 --project=../ self.jl
The solver is the RSCG
num. of Matsubara freqs. 40
  1.603179 seconds (1.33 M allocations: 126.028 MiB, 0.79% gc time, 25.78% compilation time)
1-th -1.873023776160689 1.6271730164203106
  0.192392 seconds (42.37 k allocations: 38.468 MiB, 2.00% gc time)
2-th -1.4418945612222904 0.20965023882070202
  0.270357 seconds (42.38 k allocations: 38.491 MiB, 0.99% gc time)
3-th -1.2113675935614676 0.13563888454719064
  0.344826 seconds (42.37 k allocations: 38.471 MiB, 0.30% gc time)
4-th -1.0683848898677126 0.09725266043470575
  0.424078 seconds (42.38 k allocations: 38.480 MiB, 0.29% gc time)
5-th -0.9719173169890015 0.07520905143665893
  0.482996 seconds (42.39 k allocations: 38.503 MiB, 0.26% gc time)
6-th -0.9031937676953646 0.0614030514203067
  0.524416 seconds (42.37 k allocations: 38.472 MiB, 0.27% gc time)
7-th -0.8522993398753158 0.052088904542545035
  0.567269 seconds (42.38 k allocations: 38.495 MiB, 0.22% gc time)
8-th -0.8134640044396014 0.04536935661275386
  0.602301 seconds (42.38 k allocations: 38.482 MiB)
9-th -0.7830935035525877 0.040207500802570885
  0.576101 seconds (42.38 k allocations: 38.477 MiB, 0.19% gc time)
10-th -0.7588369948792518 0.03600999926600642
  0.635709 seconds (42.38 k allocations: 38.475 MiB, 0.20% gc time)
11-th -0.7391008228944542 0.032434049930203365
  0.659751 seconds (42.38 k allocations: 38.492 MiB, 0.15% gc time)
12-th -0.7227749038266035 0.02928518035552862
  0.669727 seconds (42.38 k allocations: 38.491 MiB, 0.14% gc time)
13-th -0.709069155872567 0.026456030788754803
  0.679683 seconds (42.38 k allocations: 38.487 MiB, 0.14% gc time)
14-th -0.6974110355699145 0.02388799764731334
  0.687933 seconds (42.37 k allocations: 38.461 MiB, 0.14% gc time)
15-th -0.6873791072315169 0.0215479761914521
  0.697942 seconds (42.38 k allocations: 38.479 MiB, 0.14% gc time)
16-th -0.6786586888786034 0.019415153943493132
  0.723906 seconds (42.38 k allocations: 38.476 MiB)
17-th -0.671011679282689 0.017474092388203903
  0.726147 seconds (42.37 k allocations: 38.471 MiB, 0.16% gc time)
18-th -0.6642555683053231 0.01571143747168237
  0.731059 seconds (42.38 k allocations: 38.483 MiB, 0.16% gc time)
19-th -0.6582486676185052 0.014114520535718768
  0.726339 seconds (42.38 k allocations: 38.487 MiB, 0.16% gc time)
20-th -0.6528795616656254 0.012670927816306929
  0.730802 seconds (42.38 k allocations: 38.482 MiB, 0.15% gc time)
21-th -0.648059496125448 0.01136845977752193
  0.744700 seconds (42.37 k allocations: 38.462 MiB, 0.15% gc time)
22-th -0.643716834307313 0.010195258886917814
  0.739585 seconds (42.37 k allocations: 38.477 MiB, 0.14% gc time)
23-th -0.6397929998696529 0.00913994023009275
  0.749548 seconds (42.38 k allocations: 38.476 MiB, 0.16% gc time)
24-th -0.6362394570260124 0.00819171955905381
  0.756281 seconds (42.37 k allocations: 38.457 MiB, 0.17% gc time)
25-th -0.63301548265852 0.007340492709279601
  0.780807 seconds (42.39 k allocations: 38.504 MiB)
26-th -0.6300864856968793 0.006576883266048058
  0.787450 seconds (42.38 k allocations: 38.480 MiB, 0.13% gc time)
27-th -0.6274227436671878 0.005892252532124953
  0.782271 seconds (42.37 k allocations: 38.471 MiB, 0.14% gc time)
28-th -0.6249984477604995 0.005278692772073217
  0.775206 seconds (42.38 k allocations: 38.481 MiB, 0.14% gc time)
29-th -0.6227909814572898 0.004729001813477821
  0.781968 seconds (42.38 k allocations: 38.474 MiB, 0.14% gc time)
30-th -0.6207803405763253 0.004236642241088458
  0.786387 seconds (42.37 k allocations: 38.471 MiB, 0.14% gc time)
31-th -0.6189487129414846 0.0037957051360124226
  0.784794 seconds (42.38 k allocations: 38.479 MiB, 0.12% gc time)
32-th -0.6172801374355674 0.003400858384159586
  0.788380 seconds (42.38 k allocations: 38.498 MiB, 0.12% gc time)
33-th -0.6157602345110169 0.0030473007844825273
  0.796208 seconds (42.38 k allocations: 38.482 MiB)
34-th -0.6143759966380645 0.0027307205386831934
  0.783867 seconds (42.37 k allocations: 38.469 MiB, 0.12% gc time)
35-th -0.6131156094137203 0.0024472404597514195
  0.784699 seconds (42.38 k allocations: 38.487 MiB, 0.12% gc time)
36-th -0.6119683239104643 0.002193392792552506
  0.797728 seconds (42.38 k allocations: 38.476 MiB, 0.12% gc time)
37-th -0.6109243302390619 0.001966063888845468
  0.811322 seconds (42.37 k allocations: 38.466 MiB, 0.12% gc time)
38-th -0.6099746645043274 0.0017624669009970566
  0.798703 seconds (42.38 k allocations: 38.475 MiB, 0.18% gc time)
39-th -0.6091111276751108 0.001580108345611027
  0.804311 seconds (42.38 k allocations: 38.484 MiB, 0.14% gc time)
40-th -0.6083262105712027 0.001416755095863366
  0.802179 seconds (42.37 k allocations: 38.463 MiB, 0.14% gc time)
41-th -0.6076130419061138 0.0012704129923773262
  0.793821 seconds (42.38 k allocations: 38.480 MiB)
42-th -0.6069653199286934 0.00113929572631485
  0.803291 seconds (42.38 k allocations: 38.483 MiB, 0.12% gc time)
43-th -0.6063772737453785 0.0010218056690919832
  0.799485 seconds (42.37 k allocations: 38.476 MiB, 0.14% gc time)
44-th -0.6058436172654165 0.0009165157705025486
  0.797902 seconds (42.38 k allocations: 38.472 MiB, 0.14% gc time)
45-th -0.6053595096728797 0.0008221475370220531
  0.797214 seconds (42.37 k allocations: 38.472 MiB, 0.14% gc time)
46-th -0.6049205193354682 0.0007375584043463334
  0.797498 seconds (42.38 k allocations: 38.481 MiB, 0.13% gc time)
47-th -0.604522590849413 0.0006617273730522861
  0.800408 seconds (42.38 k allocations: 38.492 MiB, 0.12% gc time)
48-th -0.6041620163787921 0.000593740534829054
  0.795930 seconds (42.38 k allocations: 38.491 MiB, 0.11% gc time)
49-th -0.6038354056749607 0.0005327789447378642
  0.794780 seconds (42.37 k allocations: 38.463 MiB, 0.14% gc time)
50-th -0.6035396613195223 0.00047811220217767336
  0.793870 seconds (42.38 k allocations: 38.479 MiB)
51-th -0.603271956206971 0.0004290844742470022
  0.797897 seconds (42.38 k allocations: 38.477 MiB, 0.11% gc time)
52-th -0.6030297098728377 0.0003851109026020435
  0.801658 seconds (42.38 k allocations: 38.478 MiB, 0.13% gc time)
53-th -0.6028105687238502 0.0003456655189107658
  0.797855 seconds (42.38 k allocations: 38.487 MiB, 0.13% gc time)
54-th -0.6026123890946519 0.0003102796361298741
  0.797622 seconds (42.37 k allocations: 38.457 MiB, 0.10% gc time)
55-th -0.6024332167560249 0.00027853279022410063
  0.801957 seconds (42.37 k allocations: 38.468 MiB, 0.13% gc time)
56-th -0.6022712741026642 0.00025004790660750666
  0.801176 seconds (42.37 k allocations: 38.462 MiB, 0.10% gc time)
57-th -0.6021249422723195 0.00022448825972236396
  0.797505 seconds (42.38 k allocations: 38.482 MiB, 0.10% gc time)
58-th -0.6019927502404482 0.00020155103761652473
  0.796861 seconds (42.38 k allocations: 38.479 MiB)
59-th -0.6018733605571956 0.0001809666539966145
  0.803191 seconds (42.38 k allocations: 38.485 MiB, 0.13% gc time)
60-th -0.6017655584431989 0.00016249232065718839
  0.805552 seconds (42.38 k allocations: 38.483 MiB, 0.10% gc time)
61-th -0.6016682408866258 0.00014591021496670423
  0.800285 seconds (42.37 k allocations: 38.467 MiB, 0.11% gc time)
62-th -0.6015804069224595 0.0001310259234340552
  0.804897 seconds (42.38 k allocations: 38.469 MiB, 0.12% gc time)
63-th -0.6015011472832793 0.00011766420738475466
  0.796194 seconds (42.37 k allocations: 38.472 MiB, 0.12% gc time)
64-th -0.6014296388255697 0.00010566954846444353
  0.804151 seconds (42.38 k allocations: 38.481 MiB, 0.12% gc time)
65-th -0.6013651365577349 9.490118655278245e-5

となります。