$$
\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で、ハミルトニアンの期待値が計算できるようになったので、量子化学計算に初挑戦です。Quantum Native Dojoの以下の章に、VQE(Variational Quantum Eigensolver)を使って、水素化ヘリウムイオン$H-He^+$の基底エネルギーを求める例が掲載されています。これと同じことをqlazyでやってみたいと思います。
VQE (Variational Quantum Eigensolver)
時間依存を除いたシュレーディンガー方程式は、
\hat{H} \ket{\psi} = E \ket{\psi}
と記述されます。$\hat{H}$と$\ket{\psi}$は、ヒルベルト空間上で定義された演算子(行列)とベクトルなので、これを解くということは、固有値問題を解くということに等しいです。なんだこんな問題、と思われるかもしれませんが、粒子数が多くなると次元数は指数的に増大しますので、簡単に解けなくなります。そこで編み出された近似法の一つが、VQE(Variational Quantum Eigensolver)です。基本的な考え方は、とても簡単です。ある量子状態に対するハミルトニアンの期待値は、基底エネルギーよりも高くなる、つまり、
\bra{\psi} H \ket{\psi} \geq E_{0}
という性質が成り立ちますが、これを利用します。いろんな状態{$\ket{\psi_{i}}$}をランダムにもってきて、左辺を計算してみて、一番小さい値が基底エネルギーであろう、という考え方です。実際には、あるパラメータセット{$\phi_{i}$}で量子状態を$\ket{\psi(\phi_{i})}$という具合に表現しておいて、
\bra{\psi(\phi_{i})} H \ket{\psi(\phi_{i})}
を最小にするような{$\phi_{i}$}を、ベイズ最適化などの賢い方法で探索します。この探索は古典コンピュータ側でやりますので、VQEは量子と古典のハイブリッドアルゴリズムということになります。
計算実行
では、水素化ヘリウムイオンの基底エネルギーを計算します。
ハミルトニアンは、
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に従うことにします。
さて、コードは、以下の通りです。
【2021.9.5追記】qlazy最新版でのソースコードはここに置いてありま
す。
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部分を見てください。
M_str = set_hamiltonian_str()
M = Observable(M_str)
で、ハミルトニアンを文字列として記述しておいて、それを引数にして、Observableクラスのインスタンスを生成します。
init = np.random.rand(6)
callback(init)
で、$\phi_{i}$の初期値をランダムに生成して、それに対するエネルギー期待値を計算して、表示します。callback関数の中で、期待値を計算するcost関数を呼び出しています。
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によると、正解は"-2.8626207640766816"ということなので、それと比べると厳密には一致していないですが、小数点第3位まで同じなので、「はじめての量子化学計算」としては無事成功！ということで良いですよね。
おわりに
この程度の問題であれば、トータルで実質数十行くらいのコードでできてしまうことがわかり、qlazyの動作確認もできました(オブザーバブルの文字列指定が、もう少しスマートな仕様にできないかなぁ、という気がしますが、気が向いたら検討します)。
しかし、量子化学の理論や量子状態のパラメータ化について、どういう発想でこんなことをしているのかなど、よくわかっていないので、そのあたりの理解が今後の課題です。
以上