NeoAlchemist No.294

最終更新日: 2022-11-18

$$ \newcommand{\bm}[1]{{\boldsymbol{#1}}} $$

量子コンピューターで分子の基底状態のエネルギーを計算する

イントロ

 近年、量子コンピューターという言葉をよく耳にするようになりました。現在様々な研究機関や企業が量子コンピューターの開発に取り組んでおり、これまでのコンピューターではできなかったことをこの新しい量子コンピューターで実現できるのではないかという期待が高まっています。この記事では量子コンピューターを用いて分子の基底状態のエネルギーを計算する手法を解説していこうと思います。

今のコンピューターについて

 「ムーアの法則」と呼ばれる法則があります。これはインテルの創業者の一人であるゴードン・ムーアが1965年に出した論文で提唱した法則(予測)で、「半導体の集積率がおよそ2年ごとに2倍になる」というものです。集積回路が発明された1959年からそれほど時間が経っていない時期に行われた予測ですが、その後何年も同様の傾向が続いたため、「法則」という言葉が使われるようになりました。

 余談ですが、集積率を上げる理由は複数あります。まず集積率が上昇するとトランジスタ同士の間隔が密になり、半導体の性能が向上します。次に、歩留まり(良品率)が向上します。半導体を作る際、非常に高純度なシリコンウエハーを円盤状にスライスし、方眼紙状の区画に分割して回路のパターンをエッチングにより描き込むというプロセスを経るのですが、このときシリコンウエハーに不純物が存在すると半導体素子に多大な影響を及ぼし、不良品となってしまいます。このときある線幅で回路パターンを描き込む場合と、より大きな線幅で描き込む場合を比べたとき、前者より後者のほうがチップのサイズが大きくなってしまい、シリコンウエハー上の不純物に遭遇してチップが不良品になりやすくなります。逆に、チップのサイズが小さければ不良品率が下がります。そして、消費電力の低下も理由の一つです。これは集積率が上がると回路の線幅が小さくなり、回路を電気が流れることによって生じるジュール熱が小さくなるからです。

 このムーアの法則ですが、そろそろ限界を迎えるのではないかと言われています。なぜかというと、集積率が増加すれば増加するほど回路の線幅を小さくする必要に迫られるため、やがて原子のサイズに近づいてこれ以上の集積率増加が不可能になってしまうからです。現在最新のスマートフォンに利用されているプロセッサは線幅が5nmのものが主流になっています。原子のサイズは0.1nmほどであり、これ以上線幅が小さくなってくるとトンネル効果の影響が懸念されます。トンネル効果とは、本来エネルギー障壁の存在で粒子が通過できない領域を一定の確率で通過できてしまう現象のことです。回路の設計次第で絶縁体を介した電子の移動が発生し、意図しない方向に電流が流れてしまう危険があるため、さらなるプロセスの微細化は難しいだろうと考えられています。そうなれば、これまで通り微細化を追求するのではなく、全く異なる方法で計算能力を向上させることができないか、という考えになります。

 なお、実際は様々な技術的制約でプロセッサの回路線幅がすべて5nmになっているわけではないようです。半導体の集積率は半導体の性能に大きく影響するため、メーカーは半導体の性能を宣伝する際によく「〇〇nmのプロセスを利用した半導体です」とアピールしますが、その数字は必ずしも実態を反映していません。現在ではこうした「〇〇nm」という言葉は、半導体がどの世代の技術を利用して製造されているか、ということを表現する以上の意味を持たないと言われています。

量子コンピューターとは何か

 現在「量子コンピューター」と呼ばれているのものは2種類存在します。一つは量子アニーリング方式の量子コンピューター、もう一つは量子ゲート方式の量子コンピューターです。今回説明するのは量子ゲート方式の量子コンピューターで、量子アニーリング方式についての説明は割愛します。

 まず量子コンピューター上で情報をどのように表現するかについて話しましょう。よく知られているように、我々が普段用いているコンピューター(量子コンピューターに対して「古典コンピューター」と呼びます)ではあらゆる情報は0と1のビットによって表現されています。これに対し、量子コンピューターでは量子ビットを用いて次のような0と1の重ね合わせを表現することができます。

\[ \ket{\psi} = \alpha \ket{0} + \beta \ket{1} \]

ここで、\(\ket{0}\)と\(\ket{1}\)はそれぞれ量子ビットが0と1のときの状態を表現しており、状態\(\ket{\psi}\)の量子ビットを測定すると、確率\(\abs{\alpha}^2\)で0が測定され、確率\(\abs{\beta}^2\)で1が測定されます。ここで、確率の和は1でないといけないので\(\abs{\alpha}^2 + \abs{\beta}^2 = 1\)という規格化条件がかかります。

この量子ビットの状態は、\(\ket{0} = \begin{pmatrix} 1 & 0 \end{pmatrix}^T, \ket{1} = \begin{pmatrix} 0 & 1 \end{pmatrix}^T,\)というベクトルを用いて以下のようにも表現できます。

\[ \ket{\psi} = \alpha \begin{pmatrix} 1 \\ 0 \end{pmatrix} + \beta \begin{pmatrix} 0 \\ 1 \end{pmatrix} = \begin{pmatrix} \alpha \\ \beta \end{pmatrix} \]

 量子コンピューターが「古典コンピューターより優れている」と言われるのはこの重ね合わせにあり、古典ビットでは00, 01, 10, 11のビット列を個別に表現しなければならないのが、量子ビットではこれらのビット列をたったの2量子ビットで表現できます。2n個の情報を、n個の量子ビットの確率の重みを表す係数に保管することができるわけですから、古典コンピューターに比べて指数関数的な計算能力の向上を見込むことができるわけです。

 量子コンピューターはこの量子ビットを「ゲート」と呼ばれるもので操作します。量子ビットの操作はこのゲートをいろいろ組み合わせることで行います。古典コンピューターにおける、AND回路やOR回路のような論理回路に対応するものだと思っていただければ良いです。

 実例を見てもらったほうが理解が早いと思うので、以下のようなゲートを\(\ket{0}\)に作用させることを考えます。

\[ X = \begin{pmatrix} 0 & 1 \\ 1 & 0\end{pmatrix} \]

これはXゲートと呼ばれるゲートで、これを作用させると以下のようになります。

\[ X\ket{0} = \begin{pmatrix} 0 & 1 \\ 1 & 0\end{pmatrix} \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \end{pmatrix} \]

ここから、Xゲートはビットを反転させるゲートだということがわかりました。

 次に、以下のようなゲートを考えてみましょう。

\[ H = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1\end{pmatrix} \]

これはアダマール(Hadamard)ゲートと呼ばれるゲートで、これを作用させると以下のようになります。

\[ H\ket{0} = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1\end{pmatrix} \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 \\ 1 \end{pmatrix} = \frac{1}{\sqrt{2}} \qty(\ket{0} + \ket{1}) \]

アダマールゲートを作用させたあとのこの量子ビットを測定すると、50%の確率で0が得られ、50%の確率で1が得られます。他にも色々なゲートがあり、このようなゲートを組み合わせることで好きなように量子ビットを操作することができます。

 このようなゲートの大事な性質としてユニタリーであることが挙げられます。というのも、ベクトルにかかる行列がユニタリーでないとベクトルのノルムが変化してしまい、確率の総和が1でなくなってしまいます。規格化条件が常に満たされている必要があるので、直接実装できるゲートはユニタリーなゲートのみです。

行列\(U\)がユニタリーであるとは、次のような性質を\(U\)が満たすときのことを言う。 \[ UU^\dagger = U^\dagger U = I\] ユニタリー行列\(U\)の固有値の絶対値は1であり、\(||U\bm{x}||=||\bm{x}||\)

ゲートの操作を行う順番は、以下のような図式を用いて表現することができます。

量子回路の例
図1. 量子回路の例

これは量子ビットに対して

  1. アダマールゲートを作用
  2. Xゲートを作用
  3. 量子ビットを測定し、測定結果を保存
という順番で操作を行うことを表現しています。最後のメーターのアイコンは量子ビットの測定を意味します。このような図式を量子回路と呼びます。

 ここで、複数の量子ビットを扱う場合についても述べておきます。例えば0番目の量子ビットが1、1番目の量子ビットが0の状態をケットで表現すると、\(\ket{01}\)となります。量子ビットは0, 1, 2...と数えるのが一般的で、この記事では右から0番目、1番目、2番目…と数えるLittle Endianを採用することにします。左から0番目、1番目、2番目…(Big Endian)ではないことに注意してください。(人によって右から数えるか左から数えるかが異なります。)

 次のような量子回路を考えてみましょう。

量子回路の例2
図2. 量子回路の例2

このとき以下のような計算式に沿って量子ビットが操作されます。

\[ \begin{align*} X \otimes XH \ket{00} &= (X \otimes XH) (\ket{0} \otimes \ket{0}) \\ &= X\ket{0} \otimes XH\ket{0} \\ &= \ket{1} \otimes X \qty(\frac{1}{\sqrt{2}}) \qty(\ket{0} + \ket{1}) \\ &= \ket{1} \otimes \qty(\frac{1}{\sqrt{2}}) \qty(\ket{1} + \ket{0}) \\ &= \frac{1}{\sqrt{2}} \qty(\ket{10} + \ket{11}) \end{align*} \]

ここで出てきた\(\otimes\)という記号はテンソル積を表す記号で、特に行列間のものをクロネッカー積、ベクトル間のものを直積と呼ぶようです。以下のような計算規則にしたがって計算されます。

\[ \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix} \otimes \begin{pmatrix} 0 & 5 \\ 6 & 7 \end{pmatrix} = \begin{pmatrix} 1 \begin{pmatrix} 0 & 5 \\ 6 & 7 \end{pmatrix} & 2 \begin{pmatrix} 0 & 5 \\ 6 & 7 \end{pmatrix} \\ 3\begin{pmatrix} 0 & 5 \\ 6 & 7 \end{pmatrix} & 4\begin{pmatrix} 0 & 5 \\ 6 & 7 \end{pmatrix} \\ \end{pmatrix} = \begin{pmatrix} 0 & 5 & 0 & 10 \\ 6 & 7 & 12 & 14 \\ 0 & 15 & 0 & 20 \\ 18 & 21 & 24 & 28 \end{pmatrix} \]

パウリ行列の性質

 先程出てきたXゲートはパウリゲートと呼ばれるものの一つで、他に以下のようなゲートがあります。

\[Y=\begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix},\ Z=\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \]

この\(X,Y,Z\)をパウリ行列と呼びます。場合によっては\(X,Y,Z\)に単位行列\(I\)を加えた4つの行列をパウリ行列と呼ぶこともあり、任意の\(2 \times 2\)行列\(A\)は以下のように4つの行列の線形結合で表現できることが知られています。(証明は省略)

\[ A=a_I I+a_X X+a_Y Y+a_Z Z \]

このとき\(A\)がエルミートであれば係数\(a_I,a_X,a_Y,a_Z\)は実数です。

行列\(H\)がエルミートであるとは、次のような性質を\(H\)が満たすときのことを言う。 \[ H = H^\dagger\]

 また、パウリ行列を\(P(=X,Y,Z)\)とすると、以下の性質を満たします。

変分法

 次に、古典コンピューターで基底状態のエネルギーを求める方法について解説します。

 まずシュレディンガー方程式を考えます。

\[ H\ket{\psi} = E_0 \ket{\psi} \]

\(\ket{\psi}\)はハミルトニアン\(H\)の固有ベクトルであり、すべての固有ベクトルの中で一番値の小さい固有値\(E_0\)を持ちます。この\(\ket{\psi}\)を基底状態とし、\(E_0\)は基底状態のエネルギーを表します。左から\(\bra{\psi}(=\ket{\psi}^\dagger)\)をかけると、

\[ E_0 = \frac{\bra{\psi}H\ket{\psi}}{\braket{\psi}{\psi}} \]

となります。ここで、\(\ket{\psi}\)を固有ベクトルとは限らない任意のベクトル\(\ket{\phi}\)に置き換えてみると、

\[ E_\phi = \frac{\bra{\phi}H\ket{\phi}}{\braket{\phi}{\phi}} \]

となり、以下の不等式が成立します。(証明は省略)

\[ E_\phi \ge E_0 \]

 ここから、どのようなベクトルを持ってきてもハミルトニアンを作用させると基底状態のエネルギー以上の値が得られることがわかります。したがって、基底状態がわからなくとも、得られるエネルギーの値が小さければ小さいほど基底状態のエネルギーに近づくはずです。これを変分法と呼びます。

VQEアルゴリズム

 量子コンピューターで基底状態のエネルギーを計算するときも、先程と似た方法をとります。

  1. まず、初期状態\(\ket{\psi}\)を用意する。
  2. パラメータ\(\bm{\theta}\)の導入された量子回路\(U(\bm{\theta})\)を用意する。
  3. パラメータを設定し、状態\(U(\bm{\theta})\ket{\psi} = \ket{\psi (\bm{\theta})}\)を計算する。
  4. 量子ビットを測定し、期待値\(E_\bm{\theta} = \bra{{\psi (\bm{\theta})}} H \ket{{\psi (\bm{\theta})}}\)を計算する。
  5. 期待値\(E_\bm{\theta}\)が最小になるまで、ステップ3,4を繰り返す
以上の操作を量子コンピューター上で行い、分子の基底状態のエネルギーを求めることができます。これをVQE(Variational Quantum Eigensolver)と言います。

 ここからはVQEアルゴリズムの詳細を理解するため、ここから分子の電子状態をどのように量子コンピューター上で表現するか、そしてどのようにしてパラメータを導入した量子回路\(U(\bm{\theta})\)を用意するのかについて解説していきます。

スレーター行列式(第一量子化)

 まず、一般の分子のハミルトニアンを考えてみましょう。

\[ \hat{H}=-\sum_i^N \frac{1}{2} \nabla_i^2-\sum_A^M \frac{1}{2 M_A} \nabla_A^2-\sum_i^N \sum_A^M \frac{Z_A}{\left|\boldsymbol{r}_i-\boldsymbol{R}_A\right|}+\sum_{i>j}^N \frac{1}{\left|\boldsymbol{r}_i-\boldsymbol{r}_j\right|}+\sum_{A>B}^M \frac{Z_A Z_B}{\left|\boldsymbol{R}_A-\boldsymbol{R}_B\right|} \]

ここでは原子単位系を用いて電子の質量\(m_e\)やディラック定数\(\hbar\)などの定数が現れないようにしています。各項は以下のような意味を持ちます。

このままでは複雑なので、原子核の位置が空間に固定されていることにして(ボルン-オッペンハイマー近似)第2項と第5項を落としてしまいましょう。すると以下のようになります。

\[ \hat{H}=-\sum_i^N \frac{1}{2} \nabla_i^2-\sum_i^N \sum_A^M \frac{Z_A}{\left|\boldsymbol{r}_i-\boldsymbol{R}_A\right|}+\sum_{i>j}^N \frac{1}{\left|\boldsymbol{r}_i-\boldsymbol{r}_j\right|} \]
各文字の定義:
  • \(N\): 電子の総数
  • \(M\): 原子核の総数
  • \(M_A\): 原子核\(A\)の質量
  • \(Z_A\): 原子核\(A\)の電荷
  • \(\bm{r}_i\): 電子\(i\)の座標
  • \(\bm{R}_A\): 原子核\(A\)の座標

やっかいなのは残った第4項です。このクーロン反発項のため、シュレディンガー方程式は解析的に解くことが難しくなっています。ここで見方を変え、「ある電子が分子の中を運動するとき、他のN-1個の電子雲の中を運動する」といった捉え方をすると、電子は他の電子が作る平均化された電場の中を運動することになります。このとき\(i\)番目の電子が受ける平均場ポテンシャル\(V_\mathrm{eff} (\bm{r_i})\)は\(i\)番目の電子の1電子波動関数\(\phi_i (\bm{r}_i)\)を用いて

\[ V_\mathrm{eff} (\bm{r_i}) = \sum_{j \neq i}^N \int d\bm{r}_j \frac{1}{\left|\boldsymbol{r}_i-\boldsymbol{r}_j\right|} \abs{\phi_j (\bm{r}_j)}^2 \]

と表せます。これを平均場近似と呼びます。

 ここで、\(H=H^{(1)}+H^{(2)}\)と分解できるハミルトニアンが存在し、\(H^{(1)}\)の固有関数を\(\psi^{(1)}\)、\(H^{(2)}\)の固有関数を\(\psi^{(2)}\)とすれば、\(\psi^{(1)}\psi^{(2)}\)が\(H\)の固有関数となることが知られています。

 平均場近似を導入したあとのハミルトニアンは

\[ \hat{H}=\sum_i^N \qty(-\frac{1}{2} \nabla_i^2-\sum_A^M \frac{Z_A}{\left|\boldsymbol{r}_i-\boldsymbol{R}_A\right|}+V_\mathrm{eff} (\bm{r_i})) \]

と各電子ごとの項に分解できるので、この近似の導入でN電子波動関数を1電子波動関数に分解することができることがわかります。

 ここで、電子はフェルミ粒子のため波動関数は2電子の交換に対して反対称である必要があります。

\[ \psi\left(\cdots, \boldsymbol{r}_i, \cdots, \boldsymbol{r}_j, \cdots\right)=-\psi\left(\cdots, \boldsymbol{r}_j, \cdots, \boldsymbol{r}_i, \cdots\right) \]

単純にN電子波動関数を1電子波動関数の積で表現しただけではこの反対称性は表現できません。行列式を使うと、この反対称性を表現できます。

\[ \psi\left(\boldsymbol{r}_1, \boldsymbol{r}_2, \cdots, \boldsymbol{r}_N\right)=\frac{1}{\sqrt{N !}}\left|\begin{array}{cccc} \phi_1\left(\boldsymbol{r}_1\right) & \phi_1\left(\boldsymbol{r}_2\right) & \cdots & \phi_1\left(\boldsymbol{r}_N\right) \\ \phi_2\left(\boldsymbol{r}_1\right) & \phi_2\left(\boldsymbol{r}_2\right) & \cdots & \phi_2\left(\boldsymbol{r}_N\right) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_N\left(\boldsymbol{r}_1\right) & \phi_N\left(\boldsymbol{r}_2\right) & \cdots & \phi_N\left(\boldsymbol{r}_N\right) \end{array}\right| \]

この行列式で表した波動関数は他に以下の性質も満たします。

したがって、この行列式を用いた波動関数は多電子系の波動関数の表式として適切なのです。この波動関数をスレーター行列式と呼びます。

第二量子化

 ここまでスレーター行列式について解説してきましたが、スレーター行列式を作る際に各電子に\(\boldsymbol{r}_1,\ \boldsymbol{r}_2,\ \boldsymbol{r}_3, \cdots\)とそれぞれ異なるラベルを割り当てました。電子は量子論によれば本来区別のできない粒子であり、別々にラベルを割り当てることは少々不自然に見えます。(結果として区別できなくはなっていますが。)

 そこで、異なるアプローチを取ることにしましょう。このプロセスは第二量子化と呼ばれ、スレーター行列式の構成で行った量子化に対して「2番目」なので「第二」量子化という言葉が用いられています。また、先程のスレーター行列式を構成するプロセスは「第二量子化」に対する言葉として「第一量子化」と表現されることもあるようです。

 第二量子化の込み入ったプロセスは省略し、ここではその結果を利用する形をとります。パウリの排他原理より分子の各スピン軌道を占有する電子の数は高々1つですから、分子の基底状態と全ての励起状態を並べ、占有されていなければ「0」を、占有されていたら「1」を割り当てることですべての電子状態をビット列を用いて表現することができます。この各軌道の占有数をビット列で指定した\(\ket{0011}\)といった状態をFock状態と呼び、その中でも分子の基底状態を表す場合はHartree-Fock状態と呼びます。

 水素分子でSTO-3G基底関数系を用いた場合を考えましょう。このとき水素の1s軌道のみを考慮するので、スピン軌道は

\[ \sigma_{\mathrm{u} \downarrow}, \sigma_{\mathrm{u} \uparrow}, \sigma_{g \downarrow}, \sigma_{g \uparrow} \]

の4種類になります。このFock状態は各スピン軌道の占有数をnとして、

\[ \left|n_{\sigma_{u \downarrow}}, n_{\sigma_{u \uparrow}}, n_{\sigma_{g \downarrow}}, n_{\sigma_{g \uparrow}}\right\rangle \rightarrow\left|q_3, q_2, q_1, q_0\right\rangle \quad\left(q_i \in 0,1\right) \]

で表現でき、この場合のHartree-Fock状態は\(\ket{0011}\)です。このようにビット列で表現した電子状態をそのまま量子ビットにエンコードすることをJordan-Wigner変換と呼びます。

 第二量子化では、占有数表示を記述するため

を導入します。これを数式で表すと

\[ \begin{aligned} &a_p\left|n_{M-1}, n_{M-2}, \ldots, n_0\right\rangle =\delta_{n_p, 1}(-1)^{\sum_{i=0}^{p-1}n_i} \left|n_{M-1}, n_{M-2}, \ldots, n_p \oplus 1, \ldots, n_0\right\rangle \\ &a_p^{\dagger}\left|n_{M-1}, n_{M-2}, \ldots, n_0\right\rangle =\delta_{n_p, 0}(-1)^{\sum_{i=0}^{p-1} n_i}\left|n_{M-1}, n_{M-2}, \ldots, n_p \oplus 1, \ldots, n_0\right\rangle \end{aligned} \]

となります。ここで、\(\oplus\)はmod 2での足し算を行う排他的論理和を表しています。Jordan-Wigner変換のもとでは、以下のように表現できます。

\[ \begin{aligned} a_p &= \frac{1}{2} (X_p + iY_p) \otimes Z_{p-1} \otimes \cdots \otimes Z_0 \\ a_p^{\dagger} &= \frac{1}{2} (X_p - iY_p) \otimes Z_{p-1} \otimes \cdots \otimes Z_0 \end{aligned} \]

これは、粒子の生成消滅部分を \[\frac{1}{2} (X_p + iY_p) = \ketbra{0}{1},\ \frac{1}{2} (X_p - iY_p) = \ketbra{1}{0}\] と書け、反対称性の起源である\((-1)^{\sum_{i=0}^{p-1}n_i}\)をパウリZ行列で表現できるからです。

 また、第二量子化においてハミルトニアンは以下のように書けることが知られています。

\[ H=\sum_{p, q} h_{p q} a_p^{\dagger} a_q+\frac{1}{2} \sum_{p, q, r, s} h_{p q r s} a_p^{\dagger} a_q^{\dagger} a_r a_s \]

ここで

\[\begin{aligned} h_{p q} &=\int d \bm{r} \phi_p^*(\bm{r})\left(-\frac{\nabla^2}{2}-\sum_A^M \frac{Z_A}{\left|\bm{r}-\boldsymbol{R}_A\right|}\right) \phi_q(\bm{r}) \\ h_{p q r s} &=\int d \bm{r}_1 d \bm{r}_2 \frac{\phi_p^*\left(\bm{r}_1\right) \phi_q^*\left(\bm{r}_2\right) \phi_r\left(\bm{r}_2\right) \phi_s\left(\bm{r}_1\right)}{\left|\bm{r}_1-\bm{r}_2\right|} \end{aligned}\]

であり、これらは事前に古典コンピューターで計算できます。よって、ハミルトニアンが複数のパウリ行列の積\(P_i\)と係数\(c_i\)を用いて

\[H=\sum_i c_i P_i\]

と書けることがわかります。

エネルギー期待値の計算

 先程の式から、

\[E=\langle\psi(\boldsymbol{\theta})|H| \psi(\boldsymbol{\theta})\rangle=\sum_i c_i\left\langle\psi(\boldsymbol{\theta})\left|P_i\right| \psi(\boldsymbol{\theta})\right\rangle\]

と書けることがわかりますから、あとは波動関数で挟まれたパウリ行列の積の期待値を求めれば良いです。実際の計算は大変なので省略しますが、期待値の計算方法を簡単な例で示しておきます。

 例えば、パウリZ行列を\(\ket{\psi}\)(適当な状態の1量子ビット)で挟んだ期待値\(\langle\psi|Z| \psi\rangle\)を求めたいとき、

\[\begin{align*} E_\psi(Z) &=\langle\psi|Z| \psi\rangle \\ &=\langle\psi|(\ketbra{0}{0}-\ketbra{1}{1})| \psi\rangle \\ &=\braket{\psi}{0}\braket{0}{\psi}-\braket{\psi}{1}\braket{1}{\psi} \\ &=\abs{\braket{0}{\psi}}^2-\abs{\braket{1}{\psi}}^2 \\ &=P_{\ket{\psi}}(0)-P_{\ket{\psi}}(1) \end{align*}\]

より、\(\ket{\psi}\)の測定結果が0だったときの確率から測定結果が1だったときの確率を引くことで求められます。また、\(\langle\psi|X| \psi\rangle\)は\(X=HZH\)という関係を利用し、

\[\begin{align*} E_\psi(X) &=\langle\psi|X| \psi\rangle \\ &=\langle\psi|H Z H| \psi\rangle \\ &=\left\langle H^{\dagger} \psi|Z| H \psi\right\rangle \\ &=\langle H \psi|Z| H \psi\rangle \\ &=E_{H \psi}(Z) \\ &=P_{H|\psi\rangle}(0)-P_{H|\psi\rangle}(1) \end{align*}\]

のようにして\(H\ket{\psi}\)の測定結果が0だったときの確率から測定結果が1だったときの確率を引くことで求められます。ここでは1量子ビットの期待値の求め方を示しましたが、2個以上の量子ビットの期待値も同様にして求められることがわかると思います。

パラメータを導入した波動関数の生成

 先程のハミルトニアンの式を再度示します。

\[ H=\sum_{p, q} h_{p q} a_p^{\dagger} a_q+\frac{1}{2} \sum_{p, q, r, s} h_{p q r s} a_p^{\dagger} a_q^{\dagger} a_r a_s \]

この式をよく見ると、一電子励起演算子と二電子励起演算子の和で組み合わさっていることがわかります。したがって、Hartree-Fock状態とその励起状態で展開された波動関数

\[|\mathrm{CI}\rangle=c_0|\mathrm{HF}\rangle+\sum_{i, a} c_{i a} a_a^{\dagger} a_i|\mathrm{HF}\rangle+\sum_{i, j, a, b} c_{i j a b} a_a^{\dagger} a_b^{\dagger} a_j a_i|\mathrm{HF}\rangle+\cdots\]

の係数を変分法で最適化すれば、正しい答えに近づくことがわかります。これを配置間相互作用法(Configuration Interaction Method)と呼びます。もしあらゆる励起配置、すなわちN電子系においてすべてのN電子励起状態を考慮してパラメータを最適化すれば、基底関数の張る空間で最良の解を与えることが知られています。これをFull-CI法と呼びます。しかしながら、全ての励起状態を考慮することは計算量が膨大になってしまうため単純な分子を除いて不可能です。

 計算量を削減する方法として、一電子励起(single excitation)まで考慮するCIS法や二電子励起(double excitation)まで考慮するCISD法があります。ただ、このように途中で展開を打ち切ると、「サイズに対する無矛盾性(size consistency)が満たされない」という問題が発生します。

 例を挙げましょう。ある分子1つのエネルギーをCISD法で計算したとします。次に、同じCISD法で同じ分子を2つに増やして計算しても、エネルギーは2倍になりません。なぜなら、分子1つのときはその分子が二電子励起された状態まで考慮されますが、分子2つのときは分子2つを合わせて二電子励起したときの状態までしか考慮されないからです。

 この「サイズに対する無矛盾性」を満たす計算手法として結合クラスター法(Coupled Cluster Method)と呼ばれる方法があります。この結合クラスター法では、励起演算子\(T\)を指数の肩に載せたものをHartree-Fock状態に作用させ、パラメータを変分法によって最適化するという手法を取ります。

\[ |\mathrm{CC}\rangle =\exp (T)|\mathrm{HF}\rangle \] \[ T=\sum_{i, a} t_{i a} a_a^{\dagger} a_i +\sum_{i, j, a, b} t_{i j a b} a_a^{\dagger} a_b^{\dagger} a_j a_i+\cdots \]

指数関数の導入により、二電子励起までで打ち切っても指数関数の展開によってより高次の励起についての効果も取り込めるため、CI法と異なり「サイズに対する無矛盾性」が満たされます。先程と同様に、一電子励起まで考慮するとCCS法、二電子励起まで考慮するCCSD法となります。

 さて、この結合クラスター法を量子コンピューター上で実装しようとしたとき困ることがあります。それは、この演算子\(\exp (T)\)がユニタリーではないという問題です。この問題を解決するため、少し演算子を書き換えて以下のようにします。

\[ |\mathrm{UCC}\rangle=\exp \left(T-T^{\dagger}\right)|\mathrm{HF}\rangle \] \[ T=\sum_{i, a} t_{i a} a_a^{\dagger} a_i+\sum_{i, j, a, b} t_{i j a b} a_a^{\dagger} a_b^{\dagger} a_j a_i+\cdots \]

このとき\(\left(T-T^{\dagger}\right)^\dagger = -\left(T-T^{\dagger}\right)\)ですから、\(\exp \left(T-T^{\dagger}\right)\)はユニタリーです。これにより量子コンピューター上での実装が可能になりました。これをユニタリー結合クラスター法(Unitary Coupled Cluster Method)と呼びます。

 さて、これでようやくパラメータを導入した波動関数を用意できたわけですが、まだ問題が残っています。次の項目でこれを見ていきましょう。

Lie-Trotter-Suzuki分解

 第二量子化の項目で説明したように、励起演算子\(T\)も複数のパウリ行列の積\(P_i\)と係数\(c_i\)を用いて

\[T = \sum_i c_i P_i\]

と表すことができます。すると、

\[\exp \qty(T-T^{\dagger}) = \exp \qty(\sum_i c_i \qty(P_i - P_i^\dagger))\]

となります。

 残念ながら、この演算子を直接ゲートの組み合わせで実装することはできません。これは一般に行列の積が可換でないことに起因しています。単なる実数であれば

\[e^{a+b} = e^ae^b\]

のように積の形に分解できますが、これが\(A, B\)という非可換な行列の場合、

\[e^{A+B} = e^Ae^B\]

は成立しません。しかしながら、以下のような式は成立します。

\[{\displaystyle e^{A+B}=\lim _{n\rightarrow \infty }(e^{\frac{A}{n}}e^{\frac{B}{n}})^{n}}\]

ここから、計算の精度を上げたければ行列の係数が十分小さくなるまで分解すれば良いことがわかります。この分解をLie-Trotter-Suzuki分解と呼びます。

Rx・Ry・Rzゲート

 指数の肩にパウリ行列を載せた演算子をどのようにゲートで実装するか疑問に思われると思うので、この項目ではその方法を解説します。パウリ行列の項目で説明したように、パウリ行列を\(P(=X,Y,Z)\)とすると、パウリ行列は以下の性質を満たします。

ここから、\(\exp (-i \theta P)\)を計算すると、

\[\begin{aligned} \exp (-i \theta P) &=I-i \theta P+\frac{1}{2 !}(i \theta P)^2-\frac{1}{3 !}(i \theta P)^3+\frac{1}{4 !}(i \theta P)^4 \cdots \\ &=I-\frac{\theta^2}{2 !} I+\frac{\theta^4}{4 !} I+\cdots-i \theta P+i \frac{\theta^3}{3 !} P-i \frac{\theta^5}{5 !} P-\cdots \\ &=I\left(1-\frac{\theta^2}{2 !}+\frac{\theta^4}{4 !}+\cdots\right)-i P\left(\theta-\frac{\theta^3}{3 !}+\frac{\theta^5}{5 !}-\cdots\right) \\ &=I \cos \theta-i P \sin \theta \end{aligned}\]

となります。このような操作を行うゲートを\(P=X\)のときRxゲート、\(P=Y\)のときRyゲート、\(P=Z\)のときRzゲートと呼び、行列の形で表すと以下のようになります。

\[\begin{align*} R x(\theta)&=\exp \left(-i \frac{\theta}{2} X\right)=\left(\begin{array}{cc} \cos \left(\frac{\theta}{2}\right) & -i \sin \left(\frac{\theta}{2}\right) \\ -i \sin \left(\frac{\theta}{2}\right) & \cos \left(\frac{\theta}{2}\right) \end{array}\right) \\ Ry(\theta)&=\exp \left(-i \frac{\theta}{2} Y\right)=\left(\begin{array}{cc} \cos \left(\frac{\theta}{2}\right) & -\sin \left(\frac{\theta}{2}\right) \\ \sin \left(\frac{\theta}{2}\right) & \cos \left(\frac{\theta}{2}\right) \end{array}\right) \\ Rz(\theta)&=\exp \left(-i \frac{\theta}{2} Z\right)=\left(\begin{array}{cc} e^{-i \frac{\theta}{2}} & 0 \\ 0 & e^{i \frac{\theta}{2}} \end{array}\right) \end{align*}\]

Rx・Ry・Rzゲートでは三角関数の角度は\(\theta\)ではなく\(\displaystyle\frac{\theta}{2}\)で定義することが多いです。量子コンピューターではこのようなゲート操作も可能です。

CNOTゲート

 前項目では指数の肩にパウリ行列が1つ乗っているものを考えましたが、パウリ行列の積の場合はどうでしょうか? 例えば、\(\exp (-i \theta Z_1 \otimes Z_0)\)について考えてみましょう。これを計算すると、以下のようになります。確かめてみてください。

\[\exp (-i \theta Z_1 \otimes Z_0) = (I_1 \otimes I_0) \cos \theta-i (Z_1 \otimes Z_0) \sin \theta\]

このような操作を行う量子回路を実装する際、複数の量子ビットにまたがって操作を行う必要があるためこれまで紹介してきたような単一の量子ビットにのみ作用するゲートだけでは構成できません。そこで、以下のCNOTゲートと呼ばれるゲートを用います。

CNOTゲート
図3. CNOTゲート

CNOTとはControlled NOTの略称で、上の図で制御(Control)を表す小さな丸に0のビットが入力されたら何もせず、1のビットが入力されたらプラスのアイコンが書かれた大きな丸を通る標的(Target)ビットを反転させるというゲートです。数式で表現すると、以下のようになります。

\[CNOT = I_1 \otimes \ketbra{0}{0} + X_1 \otimes \ketbra{1}{1}= \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \end{pmatrix}\]

先程のRzゲートをCNOTゲートで以下のように挟んでみましょう。

Rzzゲート
図4. Rzzゲート

すると、

\[\begin{aligned} &CNOT (I_1 \otimes e^{-i\theta Z_0 / 2}) CNOT \\ &=I_1 \otimes \ketbra{0}{0}e^{-i \theta Z_0 / 2}\ketbra{0}{0}+X_1 \otimes \ketbra{1}{1}e^{-i \theta Z_0 / 2}\ketbra{0}{0} \\ &\qquad\qquad + X_1 \otimes\ketbra{0}{0}e^{-i \theta Z_0 / 2}\ketbra{1}{1}+I_1 \otimes \ketbra{1}{1}e^{-i \theta Z_0 / 2}\ketbra{1}{1} \\ &=e^{-i \theta / 2}(I_1 \otimes\ketbra{0}{0})+e^{i \theta / 2} (I_1 \otimes \ketbra{1}{1}) \\ &=e^{-i \theta / 2}\left(I_1 \otimes \frac{I_0+Z_0}{2}\right)+e^{i \theta / 2}\left(I_1 \otimes \frac{I_0-Z_0}{2}\right) \\ &=\left(I_1 \otimes I_0\right) \cos \theta+i \left(Z_1 \otimes Z_0\right) \sin \theta \end{aligned}\]

となり、目的の操作が可能です。

 過程は省きますが、\(\exp (-i \theta Z_1 \otimes X_0)\)のような場合では\(HZH=X\)の関係を利用して以下のように実装できます。

Rzxゲート
図5. Rzxゲート

アウトロ

 ここまで読んでいただきありがとうございました。まだまだ書きたいことがありましたが、疲れてしまったのでこのあたりで終わりとします。冒頭で量子コンピューターは古典コンピューターより優れている、みたいな口調で書きましたが実のところそうでもありません。VQEアルゴリズムの項目を読むとわかるように、エネルギーを求めること自体は量子コンピューターで行いますが、パラメータを少しづつずらして最適化する作業は古典コンピューターの役目となっています。それぞれ得意不得意があり、古典コンピューターでできる計算がそのまま量子コンピューターでできるわけではないのです。

 また、本文では述べませんでしたが量子コンピューターでの計算を行う際に出てくるノイズの影響がまだまだ大きいです。量子ビットは繊細であり、環境からの影響を受けて容易に状態が変化してしまいます。ノイズの程度は量子ビットそのものの操作のしやすさとトレードオフの関係にあり、ノイズが小さいということは周辺環境の影響を受けづらいということですから、当然外部から量子ビットを操作することが困難になるわけです。現在の量子コンピューターはNISQデバイス(Noisy Intermediate-Scale Quantum device)と呼ばれ、量子コンピューターが現実の問題を解くには計算中に発生したノイズをその都度修正する誤り訂正技術の発達を待つ必要があります。

参考文献

[1] 西久保靖彦(2021)『図解入門 よくわかる最新半導体の基本と仕組み』秀和システム.
[2] Michael A. Nielsen, Isaac L. Chuang. Quantum computation and quantum information 10th Anniversary Edition. 2010.
[3] 杉﨑研司(2020)『量子コンピュータによる量子化学計算入門』KS化学専門書.
[4] 原田義也(2007)『量子化学 下巻』裳華房.