Development of Accurate Electron Correlation Calculation Method Using Annealing Machine and Machine Learning

Abstract

The Xia–Bian–Kais (XBK) transformation, which maps the second quantized Hamiltonian into the Ising Hamiltonian, is a method for optimizing the configuration interaction (CI) wave function using an annealing computer. In this method, r replicas of qubit sequences, each of which corresponds to a Slater determinant, are used for mapping. To express a CI wave function with determinants having different magnitudes of coefficients, r should be increased. We developed a low-cost and highly-accurate ABCI (Annealing + Bayesian-optimization CI) method by introducing weight bits into the XBK method, which are determined by the Bayesian optimization.

Translated Abstract

The Xia–Bian–Kais (XBK) transformation, which maps the second quantized Hamiltonian into the Ising Hamiltonian, is a method for optimizing the configuration interaction (CI) wave function using an annealing computer. In this method, r replicas of qubit sequences, each of which corresponds to a Slater determinant, are used for mapping. To express a CI wave function with determinants having different magnitudes of coefficients, r should be increased. We developed a low-cost and highly-accurate ABCI (Annealing + Bayesian-optimization CI) method by introducing weight bits into the XBK method, which are determined by the Bayesian optimization.

シュレディンガー方程式を厳密に解きその波動関数を得ることは理論計算化学の大きな目標のひとつである． この厳密解は，与えられた基底の範囲では完全配置間相互作用(full CI)波動関数により表される．CI波動関数は，所与のM個の軌道の占有状態を表す電子配置$|n_1{n}_2\cdots n_M⟩$とそれに対するCI係数$c_{n_1}{}_{n_2\cdots n_M}$を用いて，以下のように表現できる．   

$|\Psi ⟩=\sum c_{n_1{n}_2\cdots n_M}|n_1{n}_2\cdots n_M⟩$

ここで電子配置は，電子が軌道n_iを占有していれば1，非占有であれば0のようにビットで表現する．N電子系の波動関数であれば，$\sum _{i=1}^M{n}_i=N$である．一般にfull CI波動関数を求める計算に要する時間は，系の電子数の階乗に比例して増大するため，古典計算機を用いた計算は十数電子系が限界である．そこで，解ける問題は限定的でありながらも，古典計算機よりも速く計算を行える可能性がある量子コンピュータをはじめとするnon-von Neumann型計算機を利用した量子化学計算に注目が集まっている [1]．

Isingハミルトニアンの最小化問題を解くアニーリング計算機を用いてCI計算を行う方法のひとつに，Xia–Bian–Kais (XBK)変換法がある [2, 3]．XBK変換法ではまず，第二量子化されたハミルトニアンをJordan–Wigner変換などによりスピンハミルトニアンに置き換える．さらに，スピンiのPauli演算子(${\sigma }_x^i$, ${\sigma }_y^i$, ${\sigma }_z^i$)をIsingハミルトニアンで用いられるz型のPauli演算子のみに置き換える以下のマッピングを行う．   

$\begin{array}{c}{\sigma }_x^i\to \frac{1-{\sigma }_z^{i_j}{\sigma }_z^{i_k}}{2},{\sigma }_y^i\to i\frac{{\sigma }_z^{i_k}-{\sigma }_z^{i_j}}{2},\\ {\sigma }_z^i\to \frac{{\sigma }_z^{i_j}+{\sigma }_z^{i_k}}{2},I^i\to \frac{1+{\sigma }_z^{i_j}{\sigma }_z^{i_k}}{2}\end{array}$

すなわち，軌道の占有状態に対応する元のスピン${\sigma }^i$を，zのみで表されるr組のスピン列(配置ビットと呼ぶ)のレプリカを用意して表現する．これは，各配置に対応するレプリカの数でCI係数の絶対値を表現していることに他ならない．CI係数の符号については，スピン${\sigma }_z^{i_j}$ごとに符号を適当に割り当てて全探索を行い，エネルギーが最も低くなるものを採用する．しかし，多くの場合，CI波動関数は係数の大きいいくつかの主配置と係数の小さい多数の配置の和で表されるため，このような波動関数を表現するにはrを大きくしなければならなかった．

そこで本研究では，XBK法に各電子配置の重みを表現する重みビットを導入したABCI (Annealing + Bayesian-optimization CI)法を開発した．Figure 1にABCI法で使われる重みビット・配置ビットと波動関数の関係を模式的に示す．配置ビットは1行が一つの電子配置に対応し，r個のレプリカを用いてCI波動関数を表す．XBK法では，各行に+1または−1の符号を先に与え，エネルギーが最小になるように配置ビットがアニーリング計算により決定される．これをすべての符号の可能性に対して実施し，最もエネルギーが低かったものを最適波動関数として採用する．ABCI法では，配置ビットを表す各行に対し，より柔軟な係数を与える．係数の絶対値は2ⁿに対応する重みビットのワンホットバイナリ表現で表す．例えば，重みビット数3の場合で，Figure 1のように重みビットと配置ビットが決まった場合，波動関数としては5Φ₁₂−Φ₃₄を表していることに相当するが，同じレプリカ数を用いたXBK法ではこのような波動関数は表せず，ABCI法により柔軟な表現が可能となっていることが確認できる．しかし，ABCI法では重みビットの可能性が多数に上るため，XBK法と同じように重みビットを全探索すると大きな計算コストがかかる．そのため，これを古典計算機によるベイズ最適化で決定することとする．このようにABCI法では計算に必要なqubit数の増加を抑えつつ効率的に重みを決定できるハイブリッドなアルゴリズムを採用した．レプリカ数rが同じ場合，アニーリング計算で決定する配置ビット数はABCI法とXBK法で同じため，必要なqubit数は変化しない．また，アニーリング計算による配置ビットの最適化は固定された重みビットに対し実行されるのに対し，重みビットのベイズ最適化はこれまでに得られた最適配置ビットに対するエネルギーに基づき行われる．すなわち，固定された配置ビットに対して重みビットを最適化している訳ではないことに注意されたい．

Figure 1.

 Schematic representation of the ABCI method.

実装は，参考文献 [4]で公開されているXBK法のプログラムを改変して行った．具体的には，積分計算はpySCF [5]を，スピンハミルトニアンへの変換はOpenFermion [6]を，ベイズ最適化はCOMBO [7]を用いた．重みビットを決定するベイズ最適化では，まず10回のランダムサーチを実行後，Thompsonサンプリングで50回の探索を行い，最小エネルギーの結果を採用した.アニーリング計算は，D-Wave Ocean SDK [8]のsimulated annealingを用いて行った．

まずは，同じqubit数を利用するABCI法でXBK法よりも柔軟な波動関数を表せることを，H₂分子のポテンシャルエネルギー曲線の計算で確認した(Figure 2)．表せる波動関数がどのように変わるかが明確にわかるよう，レプリカ数r = 3，ABCI法の重みビット数は2と小さく設定した．基底関数はSTO-3Gであり，CASCIの結果が本手法の目標となる．CASCIの活性空間は2電子2軌道，すなわちfull CIである．R_H-H ≤ 0.7 Åでは，XBK法もABCI法も単一配置のHartree–Fock (HF)波動関数が得られるが， ABCI法ではR_H-H = 0.8 Åから4Φ_gg−Φ_uuが，R_H-H = 1.4 Åから3Φ_gg−Φ_uuが得られる．この波動関数は，r = 3のXBK法では描けない．XBK法ではR_H-H = 1.2 Åから2Φ_gg−Φ_uuが単一配置よりも低エネルギーな解として得られるが，ABCI法でこの解になるのはR_H-H = 1.6–1.8 Åの間であり，この部分だけABCI法とXBK法の結果が一致する．このことからも，ABCI法により波動関数が柔軟に表現できていることが確認できる．

Figure 2.

 Potential energy curve for H₂ dissociation by HF, CASCI, XBK (r = 3), and ABCI (2 bit, r = 3) methods.

Figure 3.

 Potential energy curve for C=C bond dissociation of ethylene molecule by HF, CASCI, XBK, and ABCI methods.

次に，ABCI法の重みビット数を増加させたときの収束の様子を調べた．エチレンのC=C二重結合の解離ポテンシャルエネルギー曲線をCASCI (2, 4)，HF，XBK (r = 5)，ABCI (r = 5，重みビット数2∼10)法で計算した結果をFigure 3に示す．基底関数は6-31Gを用いた．平衡構造付近までは，XBK法ではHF解が得られているが，ABCI法では重みビット数2でもCASCIに近いエネルギーが得られている．ABCI法は，どの結合長においてもXBK法の結果よりもCASCIに近いエネルギーが得られた．ただし，重みビット数が多くなるほどにCASCIに近い結果が得られる傾向は見られるものの，必ずしもそうなっていない場合があることに注意されたい．これは主に，ベイズ最適化で最小エネルギーを与える重みビットの組み合わせが探索できていないことに由来すると考えている．本手法をより系統的な改善を得られるものにするためには，ベイズ最適化を停止する条件やアニーリング計算のパラメータなど，まだ多くの検討の余地があることもわかった．

本稿では，アニーリング計算機を用いてCI波動関数を得るXBK変換法に重みビットを導入したABCI法を提案した．実質的にXBK法の計算に必要なqubit数を削減することに成功したものである．導入した重みビット数の決定は，今回は古典計算機を用いたベイズ最適化により行ったが，最低エネルギー解が得られない場合もあり，改善の余地があることもわかった．今後は，重みビットについてもアニーリング計算により決定するなどの改良を検討している．ただし，配置ビットと重みビットを同時にアニーリング計算で最適化すると，qubit数が増加してしまうため，交互最適化の方法を模索している．また，汎用的に計算ができるようにするためには，二次化や占有電子数の制約に伴うペナルティなどのアニーリング計算のパラメータを系に合わせて自動的に最適化する手法の開発なども必要であり，検討を進めていきたい．

謝辞

本研究は，北海道大学機能強化プロジェクト「フォトエキサイトニクス研究拠点~光励起状態制御の予測と高度利用~」の助成を受けて実施された．文部科学省世界トップレベル研究拠点プログラム(WPI)により設置された北海道大学化学反応創成研究拠点(ICReDD)からも支援を受けた．計算の一部は岡崎の計算科学研究センターと九州大学情報基盤研究開発センターの計算機を利用して行った．この場を借りて謝意を表したい．

参考文献

• [1]   S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, X. Yuan, Rev. Mod. Phys., 92, 015003 (2020). doi:10.1103/RevModPhys.92.015003
• [2]   R. Xia, T. Bian, S. Kais, J. Phys. Chem. B, 122, 3384 (2018). doi:10.1021/acs.jpcb.7b10371 PMID:29099600
• [3]   M. Streif, F. Neukart, M. Leib, Lect. Notes Comput. Sci., 111, 11413 (2019).
• [4]   J. Copenhaver, A. Wasserman, B. Wehefritz-Kaufmann, J. Chem. Phys., 154, 034105 (2021). doi:10.1063/5.0030397 PMID:33499641
• [5]   Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, et al., Wiley Interdiscip. Rev. Comput. Mol. Sci., 8, e1340 (2018). https://pyscf.org/doi:10.1002/wcms.1340
• [6]   J. R. McClean, N. C. Rubin, K. J. Sung, I. D. Kivlichan, X. Bonet-Monroig, Y. Cao, et al., Quantum Sci. Technol., 5, 034014 (2020). https://github.com/quantumlib/OpenFermiondoi:10.1088/2058-9565/ab8ebc
• [7]   T. Ueno, T. D. Rhone, Z. Hou, T. Mizoguchi, K. Tsuda, Mater. Discov., 4, 18 (2016). https://github.com/tsudalab/combo3/doi:10.1016/j.md.2016.04.001
• [8]   https://github.com/dwavesystems/dwave-ocean-sdk