はじめての量子化学計算：水素化ヘリウムイオンの基底エネルギーを計算してみた

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

## はじめに

自作の量子計算シミュレータ[qlazy](https://github.com/samn33/qlazy)で、ハミルトニアンの期待値が計算できるようになったので、量子化学計算に初挑戦です。[Quantum Native Dojo](http://dojo.qulacs.org/ja/latest/index.html)の以下の章に、VQE(Variational Quantum Eigensolver)を使って、水素化ヘリウムイオン$H-He^+$の基底エネルギーを求める例が掲載されています。これと同じことを[qlazy](https://github.com/samn33/qlazy)でやってみたいと思います。

* [Quantum Native Dojo: 5-1.Variational QuantumEigensolver（VQE）アルゴリズム](http://dojo.qulacs.org/ja/latest/notebooks/5.1_variational_quantum_eigensolver.html)

## VQE (Variational Quantum Eigensolver)

時間依存を除いたシュレーディンガー方程式は、

```math
\hat{H} \ket{\psi} = E \ket{\psi}
```

と記述されます。$\hat{H}$と$\ket{\psi}$は、ヒルベルト空間上で定義された演算子(行列)とベクトルなので、これを解くということは、固有値問題を解くということに等しいです。なんだこんな問題、と思われるかもしれませんが、粒子数が多くなると次元数は指数的に増大しますので、簡単に解けなくなります。そこで編み出された近似法の一つが、VQE(Variational Quantum Eigensolver)です。基本的な考え方は、とても簡単です。ある量子状態に対するハミルトニアンの期待値は、基底エネルギーよりも高くなる、つまり、

```math
\bra{\psi} H \ket{\psi} \geq E_{0}
```

という性質が成り立ちますが、これを利用します。いろんな状態{$\ket{\psi_{i}}$}をランダムにもってきて、左辺を計算してみて、一番小さい値が基底エネルギーであろう、という考え方です。実際には、あるパラメータセット{$\phi_{i}$}で量子状態を$\ket{\psi(\phi_{i})}$という具合に表現しておいて、

```math
\bra{\psi(\phi_{i})} H \ket{\psi(\phi_{i})}
```

を最小にするような{$\phi_{i}$}を、ベイズ最適化などの賢い方法で探索します。この探索は古典コンピュータ側でやりますので、VQEは量子と古典のハイブリッドアルゴリズムということになります。


## 計算実行

では、水素化ヘリウムイオンの基底エネルギーを計算します。

ハミルトニアンは、

```math
H = (-3.8503 - 0.2288 X_{1} - 1.0466 Z_{1} - 0.2288 X_{0} + 0.2613 X_{0} X_{1} + 0.2288 X_{0} Z_{1} - 1.0466 Z_{0} + 0.2288 Z_{0} X_{1} + 0.2356 Z_{0} Z_{1})/2
```

です。

また、探索の際に使う量子状態は、６つのパラメータ{$\phi_{i}$}を以下の量子回路にセットし、$\ket{00}$をこの回路に通して作成するようにします。

```
0 --RX(phi[0])--RZ(phi[1])--CX--------------------------
1 --RX(phi[2])--RZ(phi[3])--*---RZ(phi[4])--RX(phi[5])--
```

なぜ、ハミルトニアンがこうなって、量子状態をこうするのかは、一旦置いておきます。すみません、まだ理解できていません(汗)。全面的に[Quantum Native Dojo](http://dojo.qulacs.org/ja/latest/index.html)に従うことにします。

さて、コードは、以下の通りです。

【2021.9.5追記】[qlazy](https://github.com/samn33/qlazy)最新版でのソースコードは[ここ](https://github.com/samn33/qlazy/blob/master/example/py/VQE)に置いてありま\
す。

```python
import numpy as np
import scipy.optimize
from qlazypy import QState,Observable

#------------------------------------
#  functions
#------------------------------------

def set_hamiltonian_str():

    s = "{}".format(-3.8503/2)
    s += "-{}*x_1".format(0.2288/2)
    s += "-{}*z_1".format(1.0466/2)
    s += "-{}*x_0".format(0.2288/2)
    s += "+{}*x_0*x_1".format(0.2613/2)
    s += "+{}*x_0*z_1".format(0.2288/2)
    s += "-{}*z_0".format(1.0466/2)
    s += "+{}*z_0*x_1".format(0.2288/2)
    s += "+{}*z_0*z_1".format(0.2356/2)

    return s

def ExpectVal(hm,qs):

    return qs.expect(observable=hm)

def cost(phi):

    qs = QState(2)

    qs.rx(0,phase=phi[0])
    qs.rz(0,phase=phi[1])
    qs.rx(1,phase=phi[2])
    qs.rz(1,phase=phi[3])
    qs.cx(1,0)
    qs.rz(1,phase=phi[4])
    qs.rx(1,phase=phi[5])

    exp = ExpectVal(M, qs)

    del qs

    return exp
    
def callback(phi):

    print("energy = ", cost(phi))

#------------------------------------
#  main
#------------------------------------

# Hamiltonian

M_str = set_hamiltonian_str()
M = Observable(M_str)

# VQE

init = np.random.rand(6)
callback(init)
res = scipy.optimize.minimize(cost, init,
                              method='Powell', callback=callback)
del M
```

何をやっているか、少し説明します。下の方のmain部分を見てください。

```python
M_str = set_hamiltonian_str()
M = Observable(M_str)
```

で、ハミルトニアンを文字列として記述しておいて、それを引数にして、Observableクラスのインスタンスを生成します。

```python
init = np.random.rand(6)
callback(init)
```

で、$\phi_{i}$の初期値をランダムに生成して、それに対するエネルギー期待値を計算して、表示します。callback関数の中で、期待値を計算するcost関数を呼び出しています。

```python
res = scipy.optimize.minimize(cost, init,
                              method='Powell', callback=callback)
```

で、エネルギー期待値を最小にする{$\phi_{i}$}をPowell法を使って探索しています。scipyの機能を使って一発でできます(超簡単！)。

結果は、以下の通りです。

```
energy =  -2.4316103068324915
energy =  -2.7072897381116556
energy =  -2.862206867330209
energy =  -2.8625199957532996
energy =  -2.8625207446389243
```

繰り返しごとに、そのときのエネルギー期待値を順に表示して、収束の様子をテキストで可視化しています。このときは5回で収束完了しました(何度か実行しましたが、だいたい5回くらいでした)。最終結果は、"-2.8625207446389243"です。

[Quantum Native Dojo](http://dojo.qulacs.org/ja/latest/index.html)によると、正解は"-2.8626207640766816"ということなので、それと比べると厳密には一致していないですが、小数点第3位まで同じなので、「はじめての量子化学計算」としては無事成功！ということで良いですよね。


## おわりに

この程度の問題であれば、トータルで実質数十行くらいのコードでできてしまうことがわかり、qlazyの動作確認もできました(オブザーバブルの文字列指定が、もう少しスマートな仕様にできないかなぁ、という気がしますが、気が向いたら検討します)。

しかし、量子化学の理論や量子状態のパラメータ化について、どういう発想でこんなことをしているのかなど、よくわかっていないので、そのあたりの理解が今後の課題です。


以上