この記事は？

緑Julia本「1週間で学べる！ Julia数値計算プログラミング」の6日目に入れるために書いたのですが、あまりにも分量が多くなってしまって泣く泣くカットした部分です。発売から2年経ちましたし、供養のために公開します。多いので、二つの記事に分けます。

多体電子系の厳密対角化：疎行列

ここまで固体物理で登場する様々な問題を解いてきましたが、最後に量子多体問題を厳密に解く方法について紹介します。この最後のトピックスは、Juliaを使って実際の研究の現場で取り扱うような問題をどのように解くか、という例の紹介です。ですので、書かれている内容は若干難しいかもしれませんが、行っていることは線形代数と対角化がメインです。細かな流れを追わなくても構いません。その場合には数式がどのようにコードに落とされていくのか、という観点で読んでみてください。

今回扱う方法は、厳密対角化と呼ばれます。本来であれば電子系の多体問題を扱うには「第二量子化」を用いた生成消滅演算子を用いたハミルトニアンを使うのが便利です。しかし、第二量子化による生成消滅演算子を導入するのにはこの本の趣旨からは外れてしまいます。そこで、生成消滅演算子を「あるルールに従う演算」とみなして、行列として表現することにします。生成消滅演算子を行列として表現できればハミルトニアンも行列として表現できることになり、対角化して固有値を得ることができます。この節では、上で登場した強束縛模型に電子間斥力相互作用が入った模型であるハバード模型と呼ばれる模型の厳密な基底状態の波動関数を数値計算で求めてみます。その際、Juliaの多彩な機能（疎行列、多重ディスパッチ、独自型等）を利用してコードを作成します。多重ディスパッチと独自型を組み合わせることで、他のプログラミング言語でのオブジェクト指向プログラミングと同等の自由度とメンテナンス性のあるプログラミングが可能となります。

1サイト問題

電子にはスピンという自由度がありまして、アップスピンとダウンスピンという状態があります。そこで、スピンが上向きの電子（アップスピン電子）の波動関数を$\psi_{\uparrow}(x)$、スピンが下向きの電子（ダウンスピン電子）の波動関数を$\psi_{\downarrow}(x)$としましょう。さらに、連続的な空間だと大変ですので、有限の数の格子点を考えます。これは、上で登場した強束縛模型の考え方で、一つの格子点が一つの原子に対応していると考えることができます。
格子点が一つしかない時、 格子点1つにアップスピン状態とダウンスピン状態がありますから、電子は最大2個入ることができます。
電子が2個ある時、この2電子の波動関数を$\phi(\sigma_1,\sigma_2)$とします。$\sigma_1$と$\sigma_2$は$\uparrow$か$\downarrow$の2通りを取ることができます。しかし、電子はフェルミオンですので、粒子の区別をつけることができず、同じ量子状態を取ることもできませんので、$\phi(\sigma_1,\sigma_2)$は実は$\phi(\uparrow,\downarrow)$という1種類の状態しか取ることができません。この時のシュレーディンガー方程式は

\left( \epsilon_{\uparrow} + \epsilon_{\downarrow} + U \right) \phi(\sigma_1,\sigma_2) = E \phi(\uparrow,\downarrow)

となります。ここで、$\epsilon_{\sigma}$は相互作用がない時のアップスピンあるいはダウンスピンが持つエネルギーです。そして、$U$はクーロン相互作用です。アップスピンとダウンスピンが同じ格子点にある場合、お互いにクーロン斥力の相互作用が働きます。
電子が1個ある時、1電子の波動関数を$\phi(\sigma_1)$とします。電子は1個しかありませんので、$\sigma_1$は$\uparrow$と$\downarrow$のどちらも取ることができます。もちろん、相互作用する相手がいませんのでクーロン相互作用も効きません。この時のシュレーディンガー方程式は

\left( \epsilon_{\uparrow} \delta_{\sigma_1,\uparrow} + \epsilon_{\downarrow} \delta_{\sigma_1,\downarrow} \right) \phi(\sigma_1) = E \phi(\sigma_1)

です。
電子が0個の時、波動関数はありませんが、便宜上$\psi_0$という波動関数を考えてみます。この時のシュレーディンガー方程式は

\left(0 \right) \psi_0 = E \psi_0

となります。
電子数ごとにシュレーディンガー方程式を立ててみましたが、$\epsilon_{\uparrow}$など同じものが出てきていますので、一つの方程式にまとめることはできないのでしょうか。ここで${\vec \phi} = (\psi_0, \phi(\uparrow),\phi(\downarrow),\phi(\uparrow,\downarrow))^{T}$というベクトルを考えますと、

\left( \begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 0 & \epsilon_{\uparrow} & 0 & 0 \\
0 & 0 & \epsilon_{\downarrow} & 0 \\
 0 & 0 & 0 & \epsilon_{\uparrow} + \epsilon_{\downarrow} + U 
 \end{array} \right) {\vec \phi} = E {\vec \phi}

という形になりまして、$4 \times 4$の行列に関する固有値問題になります。この行列は対角行列ですから対角化をするまでもなく固有値が求まっており、あまり行列で表示した恩恵がありませんね。
ここで登場したベクトル${\vec \phi}$ですが、今後サイト数が増えたときにわかりやすくするために

\psi_0 \rightarrow |00 \rangle

\phi(\uparrow) \rightarrow |01 \rangle

\phi(\downarrow) \rightarrow |10 \rangle

\phi(\uparrow,\downarrow) \rightarrow |11 \rangle

というように0と1で表すことにします。このように書くと波動関数を二進法の数字一つで表すことができます。今は電子が最大2個ですので二進法の桁数が2になっています。

生成消滅演算子の行列表示

1サイトしかない問題は数値計算するまでもないほぼ自明な問題になっていました。本当に考えたい問題はサイト数が多い問題です。そのような問題のシュレーディンガー方程式を行列の形で表現すれば数値計算で解くことができます。しかし、上でやったように毎回シュレーディンガー方程式を書き下してその行列の表示を求めるのは面倒です。何か良い方法はないでしょうか。
ここで、生成消滅演算子というものを導入することで面倒な手続きを回避してみましょう。生成消滅演算子とは生成演算子と消滅演算子という二種類の演算子のことです。生成演算子というのは電子を生成させる演算子で、消滅演算子というのは電子を消滅させる演算子のことです。場の量子論について触れたことがない方は電子が生成する消滅すると言われてもピンとこないかもしれません。しかしここでは場の量子論について知っていなくても問題ありません。ここでは天下り的に生成消滅演算子を定義します。生成消滅演算子はフェルミオンの性質を記述するために導入された演算子でして、「同じ状態には一つしか電子は占有できない」と「粒子を入れ替えた場合には波動関数にマイナスの符号がつく」と性質を満たすように作られています。例えば、上で出てきた2電子状態の波動関数$\phi(\uparrow,\downarrow)$の電子の位置を入れ替えると$\phi(\downarrow,\uparrow)$となりますが、これは$\phi(\uparrow,\downarrow) = - \phi(\downarrow,\uparrow)$となるような波動関数になっています。今は2電子までしか考えていませんから粒子の入れ替えは単純ですが、複数の電子がいる場合には任意のどの電子2つを選んで入れ替えてもマイナスの符号をつくように波動関数を作る必要があります。これは非常に面倒ですが、生成消滅演算子を使うことでかなり見通しよく扱いやすくなります。
あるスピン状態$\sigma$の電子に関する生成演算子$c_{\sigma}^{\dagger}$と消滅演算子$c_{\sigma}$というものを考えます。この二つの演算子は

[c_{\sigma},c_{\sigma'}^{\dagger}]_+ = \delta_{\sigma \sigma'}

[c_{\sigma},c_{\sigma'}]_+ =0

[c_{\sigma}^{\dagger},c_{\sigma'}^{\dagger}]_+ =0

という条件を満たしているとします。ここで、

[A,B]_{+} \equiv AB + BA

は反交換関係と呼ばれます。これらの演算子を使って、上で定義した4つの状態を表すことにします。
まず、電子が何もいない状態$|00 \rangle$を考えます。この状態にアップスピンに関する生成演算子$c_{\uparrow}^{\dagger}$を作用させると$|01 \rangle$という状態になるとします($c_{\uparrow}^{\dagger} |00 \rangle = |01 \rangle$)。また、ダウンスピンに関する生成演算子$c_{\downarrow}^{\dagger}$を作用させると$|10 \rangle$になるとします($c_{\downarrow}^{\dagger} |00 \rangle = |10 \rangle$)。また、$|11 \rangle$という状態は、$|00 \rangle$に$c_{\uparrow}^{\dagger}$を作用させてから$c_{\downarrow}^{\dagger}$を作用させたら得られるとします（$c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |00 \rangle = |11 \rangle$）。これで三つの状態を$|00 \rangle$から定義することができました。これらの定義の他に、消滅演算子$c_{\sigma}$を$|00 \rangle$に作用させた場合には0になることも約束しておきます（$c_{\sigma} |00 \rangle = 0$）。
これらの定義と上の反交換関係を用いると、消滅演算子の挙動が決定されます。例えば、$|01 \rangle$という状態に消滅演算子$c_{\uparrow}$を作用させると$c_{\uparrow} c_{\uparrow}^{\dagger} |00 \rangle = (1 - c_{\uparrow}^{\dagger} c_{\uparrow}) |00 \rangle = |00 \rangle$となり、$|01 \rangle$からアップスピンの電子が消滅して$|00 \rangle$になります。また、同じ演算子が二つついた場合には上の反交換関係から$c_{\sigma}^{\dagger}c_{\sigma}^{\dagger}+c_{\sigma}^{\dagger}c_{\sigma}^{\dagger} = 2 c_{\sigma}^{\dagger}c_{\sigma}^{\dagger} = 0$となりますので、$c_{\uparrow}^{\dagger}|01 \rangle = 0$となります。
これでそれぞれの波動関数を生成消滅演算子を使って表すことができました。2電子状態の波動関数は$|11\rangle = c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |00 \rangle$と定義されていますが、この波動関数の2電子のスピンを入れ替えると$c_{\uparrow}^{\dagger} c_{\downarrow}^{\dagger} |00 \rangle$となります。反交換関係$c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} + c_{\uparrow}^{\dagger} c_{\downarrow}^{\dagger} =0$を使うと、$c_{\uparrow}^{\dagger} c_{\downarrow}^{\dagger} |00 \rangle = - c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |00 \rangle = - |11\rangle$となりまして、ちゃんと2つの電子の入れ替えをしたときにマイナスの符号が出ます。今は2電子までしか考えていませんが、3電子、4電子と増やしていっても、生成消滅演算子で波動関数を書いておけば、自動的にフェルミオンの性質を満たす波動関数になるため便利です。
次に、生成消滅演算子を行列で表現してみましょう。行列で表現することで数値計算で簡単に扱うことができる形になります。今、状態は$|00\rangle$、$|01\rangle$、$|10\rangle$、$|11\rangle$の4つあります。波動関数がこれらの線型結合:

|\psi \rangle = a_{00} |00\rangle + a_{01} |01\rangle + a_{10} |10\rangle + a_{11} |11 \rangle

で書かれているとしますと、この波動関数に生成演算子$c_{\uparrow}^{\dagger}$をそれぞれかけてみますと、$c_{\uparrow}^{\dagger}|00\rangle = |01 \rangle$、$c_{\uparrow}^{\dagger}|01\rangle = 0$、$c_{\uparrow}^{\dagger}|10\rangle = c_{\uparrow}^{\dagger} c_{\downarrow}^{\dagger}|00 \rangle = -|11 \rangle$、$c_{\uparrow}^{\dagger}|11\rangle = 0$となりますから、

c_{\uparrow}^{\dagger} |\psi \rangle = a_{00} |01\rangle  - a_{10} |11\rangle 

となります。よって、

C_{\uparrow}^{\dagger} 
\left( \begin{array}{c}
  a_{00} \\
  a_{01} \\
  a_{10} \\
  a_{11} 
 \end{array} \right)
 =
 \left( \begin{array}{c}
  0 \\
  a_{00} \\
  0 \\
  - a_{10} 
 \end{array} \right)
= 
\left( \begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
 0 & 0 & -1 & 0
 \end{array} \right) 
 \left( \begin{array}{c}
  a_{00} \\
  a_{01} \\
  a_{10} \\
  a_{11} 
 \end{array} \right)

となり、生成演算子の行列表示$C_{\uparrow}^{\dagger}$が求まります。同様に消滅演算子$c_{\uparrow}$をかけてみると、$c_{\uparrow}|00\rangle = 0$、$c_{\uparrow}|01\rangle = |00\rangle$、$c_{\uparrow}|10\rangle = 0$、$c_{\uparrow}|11\rangle = c_{\uparrow} c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |00 \rangle = - |10 \rangle$となりますので、その行列は

C_{\uparrow}
\left( \begin{array}{c}
  a_{00} \\
  a_{01} \\
  a_{10} \\
  a_{11} 
 \end{array} \right) = 
\left( \begin{array}{cccc}
 0 & 1 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
0 & 0 & 0 & -1 \\
 0 & 0 & 0 & 0
 \end{array} \right)
 \left( \begin{array}{c}
  a_{00} \\
  a_{01} \\
  a_{10} \\
  a_{11} 
 \end{array} \right)

となります。Juliaで書くと

julia> Cdag = [0 0 0 0
       1 0 0 0 
       0 0 0 0
       0 0 -1 0 ]
4×4 Matrix{Int64}:
 0  0   0  0
 1  0   0  0
 0  0   0  0
 0  0  -1  0
julia> C=[0 1 0 0 
       0 0 0 0 
       0 0 0 -1
       0 0 0 0]
4×4 Matrix{Int64}:
 0  1  0   0
 0  0  0   0
 0  0  0  -1
 0  0  0   0

となりますね。行列にすると$[C_{\uparrow}]^{\dagger} = C_{\uparrow}^{\dagger}$となっていることがはっきりわかります。また、$C_{\uparrow}^{\dagger} C_{\uparrow}$は

julia> Cdag*C
4×4 Matrix{Int64}:
 0  0  0  0
 0  1  0  0
 0  0  0  0
 0  0  0  1

となりますから、

C_{\uparrow}^{\dagger} C_{\uparrow} =\left( \begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 1
 \end{array} \right)

となりますね。ダウンスピンの生成演算子$c_{\uparrow}^{\dagger}$も同様にやってみると、$c_{\downarrow}^{\dagger}|00\rangle = |10 \rangle$、$c_{\downarrow}^{\dagger}|01\rangle = |11 \rangle$、$c_{\downarrow}^{\dagger}|10\rangle = 0$、$c_{\downarrow}^{\dagger}|11\rangle = 0$となりますから、

C_{\downarrow}^{\dagger} \left( \begin{array}{c}
  a_{00} \\
  a_{01} \\
  a_{10} \\
  a_{11} 
 \end{array} \right)  = 
\left( \begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 \\
 0 & 1 & 0 & 0
 \end{array} \right) \left( \begin{array}{c}
  a_{00} \\
  a_{01} \\
  a_{10} \\
  a_{11} 
 \end{array} \right) 

と書くことができまして、Juliaでは

julia> Cdagdown=[0 0 0 0
       0 0 0 0
       1 0 0 0
       0 1 0 0]
4×4 Matrix{Int64}:
 0  0  0  0
 0  0  0  0
 1  0  0  0
 0  1  0  0

となります。ですので、$C_{\downarrow}^{\dagger}C_{\downarrow}$は行列で書くと

C_{\downarrow}^{\dagger}C_{\downarrow} = 
\left( \begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 \\
 0 & 0 & 0 & 1
 \end{array} \right) 

となります。
これらの生成消滅演算子の組$C_{\sigma}^{\dagger}C_{\sigma}$の行列表示をみてわかることは、スピン状態$\sigma$がいる状態の場合には1を返し、いない場合には0を返す、という演算になっているということです。したがって、アップスピンの電子がいるときにはエネルギーが$\epsilon_{\uparrow}$増え、アップスピンとダウンスピンの両方の電子がいるときには相互作用エネルギー$U$が加えられる、という表現をこの行列を使って表すことができます。つまり、ハミルトニアンも生成消滅演算子（あるいはその行列表示）で表すことができるということです。
ハミルトニアンを

H = \sum_{\sigma} \epsilon_{\sigma} c_{\sigma}^{\dagger} c_{\sigma} + U c_{\uparrow}^{\dagger} c_{\uparrow} c_{\downarrow}^{\dagger} c_{\downarrow}

と書いてみましょう。生成消滅演算子の行列表示を代入して整理すると、上で登場した$4 \times 4$行列のハミルトニアンになることがわかります。生成消滅演算子の行列表示がわかっていれば、その組み合わせで書けるハミルトニアンは必ず行列で表現できます。行列で表現できれば固有値問題をJuliaで解ける行列の問題として定義することができます。ここで扱ったのは1サイトの問題なのでご利益がわかりにくいかもしれませんので、次の節ではサイト数を増やしてみましょう。

2サイト問題

次は格子点が2つある問題を考えてみます。1つのサイトにはアップスピンとダウンスピンの状態があり、アップスピンに入る、ダウンスピンに入る、両方に入る、両方とも入らない、の4状態がありました。2つのサイトの場合にはそれぞれに4状態がありますから、$4^2 = 16$個の状態があります。状態を書き下してみると$|0000\rangle$、$|0001\rangle$、$|0010\rangle$、$|0011\rangle$、$|0100\rangle$、$|0101\rangle$、$|0110\rangle$、$|0111\rangle$、$|1000\rangle$、$|1001\rangle$、$|1010\rangle$、$|1011\rangle$、$|1100\rangle$、$|1101\rangle$、$|1110\rangle$、$|1111\rangle$の16状態です。二進法で0から15までを4桁で書くことでこれらの状態を表現できます。ここでは、左の2つをダウンスピン、右の2つをアップスピンの状態を表しているとしましょう。つまり、$|1111\rangle = c_{2,\downarrow}^{\dagger}c_{1,\downarrow}^{\dagger}c_{2,\uparrow}^{\dagger}c_{1,\uparrow}^{\dagger}|0000 \rangle$という形で状態を定義しておきます。例えば、$|1001 \rangle$は$|1001 \rangle = c_{2,\downarrow}^{\dagger}c_{1,\uparrow}^{\dagger}|0000 \rangle$です。
1サイトの時と同じように生成消滅演算子を行列として表現してみましょう。行列表示する時のポイントは、作用させようとする生成演算子の位置です。例えば、1サイトの時には2電子状態は$|11 \rangle = c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |00\rangle$と書けていましたが、これにアップスピンに関する消滅演算子$c_{\uparrow}$を作用させると、$c_{\uparrow} |11 \rangle = c_{\uparrow} c_{\downarrow}^{\dagger} c_{\uparrow}^{\dagger} |00\rangle = - c_{\downarrow}^{\dagger} c_{\uparrow} c_{\uparrow}^{\dagger} |00\rangle$のようにマイナス符号が出ます。このマイナス符号は反交換関係$c_{\uparrow} c_{\downarrow}^{\dagger} + c_{\downarrow}^{\dagger} c_{\uparrow} =0$によって生じています。ある状態に演算子を作用させるときは、必ず一番左にその演算子がありますから、その演算子を反交換関係によって移動させる必要があります。この移動が奇数回であればマイナスが、偶数回であればプラスの符号がつきます。この符号をJuliaで計算できれば、生成消滅演算子を行列で表現できるようになります。
以後は$N$サイト問題を取り扱えるように一般的なコードを書いてみましょう。

独自型の利用とNサイト問題

それでは、Juliaでコードを書いていきましょう。まず、状態を独自の型で定義してしまいます:

mutable struct Electronbasis{N}
    state::Array{Bool,1}
    index::Int64
    function Electronbasis(N,index)
        state = i2state(index,N)
        return new{N}(state,index)
    end
end

function i2state(i,N)
    @assert 0 <= i <  2^N "index should be in  0 <= i < $(2^N)!"
    state = zeros(Bool,N)
    bit = string(i, base=2,pad=N)
    for i=1:N
        state[i]= parse(Bool,bit[N-i+1])
    end
    return state
end

ここで、Electronbasis{N}のように{}がついている場合は、この型はパラメトリック型と呼ばれます。パラメトリック型を用いますと、型がどんな型かを指定することができます。実はこれまで登場してきた配列ArrayにもArray{Bool,1}のような形で{}がついており、Array{Bool,1}であれば要素の型がBoolで1次元、Array{Float64,2}であれば要素の型がFloat64で2次元、という意味になります。今回独自定義したElectronbasis{N}では、{N}のNがパラメータです。今回の目的ではNはサイト数が入ることを想定しています。この型の中にあるstateは状態の0と1の並びを保持しておく変数、indexは状態のラベルを表します。この型の中に入っているfunctionのElectronbasis(N,index)はコンストラクタと呼ばれるものでして、どうやってこの型を作るか、ということを指定しています。今回の場合は、引数としてNとindexを取る関数でこの型を定義します、と言っています。この関数の中にはi2state(index,N)という関数が使われていますので、その挙動を説明します。@assert　0 <= i < 2^Nは引数iがこの範囲に入っていないとエラーメッセージ"index should be in 0 <= i < $(2^N)!"を返すマクロです。想定していない数字が入ってきた時にエラーで止める役割をしています。string(i, base=2,pad=N)はiという10進数の数字を桁がNの2進数の文字列に変換する関数:

julia> a = string(12, base=2,pad=5)
"01100"

です。あとはparse(Bool,bit[N-i+1])を使って、右からi番目の文字が0か1かを判定してstate[i]に格納しています。なお、Bool型はtrueかfalseを格納する型なのですが、trueは1、falseは0として言い直すことができるようになっています。
次に、作ったElectronbasis型を表示する関数を作りましょう。表示する関数としてdisplayという関数がJuliaにはあります。REPLで

julia> a = [1 2
       3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

のように行列を作るとみやすい形で表示してくれます。これはdisplay(a)という関数を呼んでいます。ですので、REPLではなくjulia hoge.jlの形で呼び出した時には

a = [1 2
    3 4]
display(a)

のようにすることで、行列をいい感じに表示してくれます。少し気をつけるべき点としては、display(a)としたときには最後に改行されませんので、改行をしたければprintln("")等を次の行に追加してください。
display(a)は変数aの型に応じていい感じに出力してくれる関数ですが、今自分で定義したElectronbasis型にはもちろん対応していません。対応させるには、多重ディスパッチを用いてdisplay(a)のaがElectronbasis型の時の挙動を定義すればよいです。Juliaの多重ディスパッチの利点として、すでにどこかで定義されている関数に新しい型の挙動を入れることが容易であることがあげられます。displayという関数は何もパッケージを読み込まなくても使える関数ですが、このような関数はBaseというパッケージに入っていることを意味しています。ですので、

function Base.display(electron::Electronbasis{N}) where N
    label = get_label(electron)
    println(label)
    return
end
function get_label(electron::Electronbasis{N}) where N
    label=""
    label*="|"
    for i=N:-1:1
        if electron.state[i]
            label*="1"
        else
            label*="0"
        end
    end
    label*=">"    
end

を書きます。これで、displayという関数の引数にElectronbasis型が来た時の挙動が定義されました。なお、ここでElectronbasis{N}とすることでElectronbasis型がパラメトリック型であることを明示しています。where NというのはNがどのような型かを指定する文法です。例えばwhere N <: IntとするとNが整数型の時だけこのdisplayの定義が使われるようにすることができます。where Nと何も型の情報を指定していない場合には、Nはなんでもよい（Any型）ことを意味します。パラメトリック型の面白いところは、ここで登場したNを関数内で普通に使えることです。例えばこの関数ではfor i=N:-1:1とすることでfor文の範囲をNを使って指定しています。
このdisplay関数を使いますと、

N = 2*2
index = 3
state = Electronbasis(N,index)
display(state)

の出力結果は|1100>となります。これでどんな状態を取り扱っているかがすぐにわかるようになりました。

さて、次はElectronbasis型に色々な機能をつけていきましょう。Pythonであればクラスに色々なメソッドをつけるのと同じような感覚で独自型に関する関数を追加していくことができます。まず、あるサイト$i$に電子がいるかどうかの関数を

function check_site(electron::Electronbasis,isite)
    return electron.state[isite]
end

と定義します。ElectronbasisのフィールドのstateはArray{Bool,1}という型ですので、要素はtrueかfalseをとります。isite番目のサイト(1サイトの系であればアップスピンの電子を1番目、ダウンスピンの電子を2番目と呼ぶことにします)の電子を消滅させる消滅演算子を作用させたときに符号がどうなるかを判定する関数を

function check_sign(electron::Electronbasis{N},isite) where N
    sign = 1
    for jsite=N:-1:isite+1
        if check_site(electron,jsite)
            sign *= -1
        end
    end
    return sign
end

で定義します。ここで、for jsite=N:-1:isite+1のループは、isiteより左側にあるサイトのうち、電子がいるサイトを数えて符号を決定しています。そして、isite番目のサイトの電子を消滅させる関数を

function erace!(electron::Electronbasis,isite)
    @assert check_site(electron,isite) "there is no electron at $(isite)-site!"
    electron.state[isite] = false
    electron.index = state2i(electron.state)
end

と定義します。@assert check_site(electron,isite)というのは、消滅させる電子がちゃんといるかどうかを判定している箇所で、電子がそのサイトにいないのにもかかわらず消滅演算子を作用させようとした場合にはエラーメッセージを出して止まるようにしています。electron.state[isite] = falseはisite番目の電子をいないことにしており、electron.index = state2i(electron.state)は新しいラベルを計算しています。このstate2iは状態が与えられた時にラベルを計算する関数で、

function state2i(state)
    N = length(state)
    ii = 0
    for i=1:N
        ii += state[i]*2^(i-1)
    end
    return ii
end

で定義しました。これでisite番目のサイトの消滅演算子を行列として作成する準備が整いました。行列を作成する関数を

function make_C(isite,nsite)
    N = 2*nsite
    ndim = 4^nsite
    mat_C = zeros(Int8,ndim,ndim)
    for index=0:ndim-1 #消した先の状態のラベル
        electron = Electronbasis(N,index)
        if check_site(electron,isite) #isiteに電子がいるかどうかチェック
            sign = check_sign(electron,isite) #消滅演算子の符号の計算
            erace!(electron,isite)
            jindex = electron.index #消えた状態のラベル
            mat_C[jindex+1,index+1] = sign #行列要素を決める
        end
    end
    return mat_C
end

と定義すればOKです。ここで、消滅演算子は$c_{\uparrow} |01 \rangle = |00 \rangle$のような形で定義されているので、$|01 \rangle$の係数$a_{01}$が$|00 \rangle$の成分に移る形になっています。したがって、$|00 \rangle$はjindexというラベル、$|01 \rangle$がindexというラベルになっています。生成演算子は消滅演算子の転置で計算できますから、これでハミルトニアンを行列表示する準備が整ったことになります。

part2に続きます。