# Fileset

[59_11_771.pdf](https://mdr.nims.go.jp/filesets/6d5c0f3f-1f96-46a8-a45e-bec7ccf1a844/download)

## Creator

[中野 晃佑](https://orcid.org/0000-0001-7756-4355)

## Rights

[In Copyright](http://rightsstatements.org/vocab/InC/1.0/)

## Other metadata

[第一原理量子モンテカルロ計算の現状と材料科学分野での応用展望](https://mdr.nims.go.jp/datasets/61cd340b-62cc-42cd-8aa4-6512d9f7c7b1)

## Fulltext

771セラミックス　59（2024）No. 111. セラミックス分野における計算科学　マテリアルズ・インフォマティクスと呼ばれる分野が立ち上がって以降，セラミックス分野においても，計算科学が果たす役割はより重要になっている．広く計算科学では，現象の時間・空間的スケールによって，異なる手法が用いられる．本寄稿が扱う計算は，その中でも，ミクロな電子状態を扱う「第一原理計算」である．セラミックス分野で最も使われている第一原理計算は，密度汎関数理論（Density Functional Theory: DFT）に基づく計算であろう．DFT では，多体問題を一体問題に落とし込んだコーン・シャム（KS）方程式を解くことによって，KS軌道とそのエネルギー（バンド分散やフェルミ面）を得る．得られた情報を元に，フォノン計算などを行い，物性計算を行うこともできる．　DFT はこれまで大成功を収めているものの，問題点の 1 つとして，結果の交換相関汎関数（eXchange-Correlation functional: XC）依存性が挙げられる．DFT においては，不明な電子相関の効果を XC に押し込めたことで計算を可能としている．XC の具体系に関しては，多くの種類が開発されており，それぞれ一長一短である．つまり，すべての物質に対して普遍的に良い結果を与える XC はない．この問題に対して，別の角度からアプローチを行う方法の一つが，筆者が専門とする第一原理量子モンテカルロ法（Quantum Monte Carlo: QMC）である 1）．2. 第一原理量子モンテカルロ法の概要　QMC は，多体シュレーディンガー方程式を，モンテカルロ法を使って解く手法であり，原理的には，XC のようなパラメタが必要ない理論の形式になっている．QMC の実装のうち，最も利用されている実空間のQMCとして，変分量子モンテカルロ法（Variational Monte Carlo: VMC） と拡散量子モンテカルロ法（Diffusion Monte Carlo: DMC）が挙げられる 1）．　VMC が依拠する原理は，変分原理である．すなわち，所与の多体波動関数とそのエネルギーをΨ� ��r rN1( ) ，および E ，真の基底状態の多体波動関数とそのエネルギーを，Ψ0 1� ��r rN( )，および E0 ，系のハミルトニアンをHと表した時，Edr dr r r H r rdr dr r rEN N NN N=∫ ( ) ( )∫ ( )≥∗1 1 11 12Ψ ΨΨ00が成り立つという原理である（等号は = 0Ψ Ψ の時成り立つ）．ここで， N は系に含まれる電子数，ri は i番目の電子の座標である．多体波動関数には変分パラメタ（e.g.，分子軌道係数）が含まれている．VMCでは，エネルギーが下がるように，そのパラメタを最適化することで，真の基底状態の波動関数とエネルギーに迫る．VMC においては，最適化の出発点となる波動関数はユーザーが与える必要がある．最もよく使われる波動関数系は，ジャストロー・スレーター型の関数であり，Ψ r r e Det F r rNJ r rNN1 11( ) = ⋅ ( )( )( )とジャストロー部分 eJ r rN� ��1( ) とスレーター行列式部分Det F r rN� ��1( )( )に分解される．ジャストロー部分は，スレーター行列を超えた電子相関を取り込むために考慮され，スレーター行列式は多体波動関数が粒子の交換に対し反対称でなければならないという要請から利用される．ジャストロー，行列式部分は，別々に複数のパラメタを持ち，変分原理にしたがってそれぞれ最適化される．実務的には，最適化の初期波動関数をDFT で作成するため，厳密に言えば結果が XC に依存する．系によっては，いくつかの初期条件を試す必要がある．　VMC において， E の評価には，多体波動関数の3N 次元における積分評価が必要になる． E の評価をどのように行うか，という点に関しては，式変形を行い，E dr dr r r e r rN N L N= ∫ ( ) ( )� � � � � �� � �1 1 1π© 日本セラミックス協会第一原理量子 モンテカルロ法の 現状と物質科学分野での応用展望─原子に働く力の計算─Current Status and Perspectives of ab initio Quantum Monte Carlo Methods in Materials ScienceKey-words：ab initio calculation, Quantum Monte Carlo, Atomic force中 野　晃 佑Kosuke NAKANO（National Institute for Materials Science（NIMS））特集　計算材料設計最前線2024772 セラミックス　59（2024）No. 11についても，現時点では VMC に限り有効である 5）．　ところで，物質科学系の研究者からよく聞かれる質問として，DFT と QMC の算出できる値の違いがある．両者とも，同じくシュレーディンガー方程式を解こうとするものであるから，基本的には同じ物性値を算出可能である．全エネルギーを基盤とする量，例えば，結合エネルギー，反応エネルギーなどは QMC でも問題なく計算可能であり，多くの QMC の物質科学分野での応用は，エネルギーに基づいた研究である．しかし，バンド分散やフェルミ面などの，（準粒子）1 体描像により得られるものは，QMC では計算が難しい．これは，QMC の結果得られる波動関数が 3N 次元の多体波動関数であって，DFT とは異なり，1 体近似された波動関数（KS 軌道）ではないという理由による．多体波動関数を準粒子描像に落とし込むことができれば，バンド分散等の描画も可能かと思われるが，今後の研究進展が待たれる．フォノン計算や分子動力学計算は，原子に働く力が確立しつつあり，DFT と同様に計算可能となっている（c.f.，フォノンは準粒子描像）．　近年では，計算パッケージも徐々に整備されてきている．現在，アクティブに開発・保守が進められている代表的な実空間 QMC パッケージを表 1 に示した．DFT と同じく，基底関数，実装されている計算可能な物性値，ライセンス形態などに違いがある．3. 物質科学分野における応用例　QMC の中でも，特に DMC は，固体周期系に対するゴールドスタンダードの計算と呼ばれ，DFT が取り扱いに難渋する物質群に対する「参照結果」として用いられてきた．これまでの応用は，例えば，高圧下における液体／固体水素の状態図解明 11），12），層状物質を含む結晶多形間の相対的安定性の議論 13），結晶表面の分子吸着に関する研究 14），二次元層状物質の剥離エネルギーに関する研究 15），分子結晶の格子エネルギーの推定 16），などが挙げられる．巻末に関連する参考文献を記載した．参考にされたい．　最近では，DFT 分野と同様に，QMC 分野でもハイスループット計算技術が確立し始めているとすることから始める．π� �� �� � � ���� �r rr rdr dr r rNNN N1121 12( ) = ( )∫ ( )ΨΨ，e r rH r rr rL NNN111( ) = ( )( )ΨΨ である．すると，この 3N 次元の積分は，マルコフ連鎖モンテカルロ法により，π� ��r rN1( ) に従う配置 R r ri N i= ( )1 を M 個作成し，その配置を使って，EMe RiML i≈ ( )=∑11と評価できる．このように，多体波動関数の積分にモンテカルロ積分を使うことから，変分量子モンテカルロ法と呼ばれている．　一方，DMCは，VMCとは異なり，多体シュレーディンガー方程式の基底状態を射影演算子によって抽出する手法である．紙面の制約もあるため，詳しいアルゴリズムの詳細に興味がある読者は総説論文 1）などを参照されたい．端的に言えば，DMC は，多体波動関数の値が 0 になる「節」と呼ばれる多次元空間上の座標を固定した上で，振幅のみを緩和させる手法である（c.f.，電子がフェルミオンであるため．節固定近似と呼ばれる）．多体波動関数の節は，通常，DFT で計算されて固定される．その意味で，DMC の結果は，厳密に言えば DFT の計算結果に依存するが，DFT の結果は節のみに影響するため，異なる XC を用いても同様の結果が得られることが多い 2）．　さて，一般的に，その厳密解への近さでは，VMCよりも DMC に軍配が上がる．VMC では，前記の通り，変分原理にしたがって多体波動関数に含まれる変分パラメタを最適化する．多体波動関数に含まれる変分「パラメタ」ということから分かるように，最終的に得られるエネルギーは，与える関数系（とパラメタ数）に制約される．無限にフレキシブルな関数系（e.g.，深層ニューラルネットワーク）を設定することで厳密解に迫ろうという試みも最近なされているが 3），4），系が大きくなると最適化の難しさとのトレードオフが顕在化し，一筋縄では行かないと筆者は考えている．一方，DMC においては，固定するのは節のみであり，より厳密解に近い状態が得られることが通常である．　しかし，実装や計算結果解釈の容易さという点では，VMC に軍配が上がる．VMC では，基底状態として得られる波動関数が陽に与えられるのに対し，DMCでは，射影演算が数値的に行われるため，基底状態の波動関数は，重みを持った粒子（ウォーカー）の分布としてのみ得られる．そのため，例えば，物性計算はVMC に限られていることが多い．本寄稿のサブタイトルである，筆者が最近確立した原子に働く力の計算表 1  現在開発・保守が進められている QMC コードのうち，周期境界条件に対応しているパッケージ．名称 主要開発国 ライセンスQMCPACK6） 米国 NCSACASINO7） 英国，ドイツ 独自TurboRVB8） 日本，フランス GPLv3CHAMP（EU）9） オランダ GPLv3PyQMC10） 米国 MIT773セラミックス　59（2024）No. 11に働く力を，原子位置の変位に対するエネルギーの全微分，� �F dE dRα αdef − / ，として定義する．VMC では，項を次のように分けて評価する 20），21），� �− −( )��F eRe EREcdcdRLLiMiiαα α α= − ∂∂∂∂− ∂∂=∑log21Ψここで，Rα は原子α の位置，Ψは多体波動関数， eLは局所エネルギーと呼ばれる量（e HL /Ψ Ψdef ），Hは多体ハミルトニアン， ci は波動関数に含まれるパラメタ（e.g.，分子軌道係数）， M はパラメタの総数である．〈⋯〉 はこれらの項がモンテカルロ法によって評価されることを示す．QMC においては，これまで，第 1 項（Hellmann–Feynman 項）と第 2 項（Pulay 項，局在基底を利用する場合）のみで力や圧力が評価され，第 3 項は無視されてきた 22）．しかしながら，2021 年頃，米国のグループ 20），および，筆者らのグループ 21）が，扱う電子数が大きくなると第 3 項の影響が無視できないことを発見した．例えば，c-BN では圧力の評価値に 5 GPa 程度の誤差が生じる 5）．第 3 項を無視したことによる誤差を self-consistency 誤差と呼ぶ 20）．　もちろん，多体波動関数に含まれる変分パラメタをすべて最適化できれば，第 3 項は無視できる．しかし，実務上，数百万を超える変分パラメタをすべて最適化することは難しい．そこで，通常，VMC においては，ジャストロー関数の部分のみを変分最適化して，スレーター行列式の部分は DFT で得られたものに固定しておく，という方策が採用される．この場合，ジャストロー部分の変分パラメタ ciJ は変分最小（ ∂ ∂ =E ciJ/ 0）（e.g.，QMCPACK に付随する Nexus17））．筆者が最近開 発 し た ソ フ ト ウ ェ ア は，TurboRVB に 対 す るTurboGenius18）というソフトウェアである．Pythonで実装されており，TurboRVB の諸機能をラップする形になっている．典型的な計算は Python 上で自動化され，ジョブの投入や監視なども行うことができる．このようなソフトウェアの登場により，複雑であったQMC の計算手順が（隠蔽によって）簡素化され，確立した計算手法であれば，ルーチンワーク的に実施可能になった．実務的観点では大きな進歩であろう．4. 原子に働く力の計算に関する最近の進展　QMC の応用は，セラミックス分野を含めた物質科学分野では，DFT ほどは広まっていないことは間違いない．計算コストが大きいという計算資源的な問題はあるものの，それよりも本質的な問題の 1 つとして，原子に働く力や圧力が計算できなかったことが大きいと筆者は考えている．物質科学分野においては，エネルギーのみならず，その応答の理解が非常に重要である．例えば，エネルギーの原子変位に対する応答，すなわち，「原子に働く力」が計算可能であると，フォノン計算や分子動力学計算を通じた，物質の安定性の議論が可能である．原子に働く力の計算は，DFT では確立された計算であるが 19），QMC では計算が困難であることが知られていた．その理由の 1 つが，self-consistency 誤差と呼ばれる誤差である 20），21）．　以下では，VMC に焦点を絞って話を進める．原子図 1　立方晶 BN の圧力を，ポテンシャルエネルギー（実線）の数値微分（波線），および微分の直接評価（菱点）から計算した結果．波線が真の値．（a）Quantum Espresso（DFT）を用いた計算結果：波線と菱点が一致している，つまり，微分の直接評価によりバイアスのない圧力が計算可能．（b）TurboRVB（VMC）を用いた計算結果：波線と菱点が一致しない，つまり，微分の直接評価によりバイアスのない圧力が計算できない．四角点は筆者が考案した方法により求めた圧力であり，正しい圧力と一致している．結晶構造はVESTA26）で描画した．（b）は文献5）より転載［Creative Commons Attribution 4.0 International license］．（a） Quantum Espresso（DFT）の計算結果 （b） TurboRVB（VMC）の計算結果774 セラミックス　59（2024）No. 11　謝　辞　本研究を進めるにあたってお世話になった方々に深く感謝致します．まずは，（故）Sandro Sorella 教授（SISSA）に感謝と， 深く哀悼の意を示したいと思います．Michele Casula 研究員（CNRS），および Giacomo Tenti 博士課程学生（SISSA）には，VMC による原子に働く力の計算について，共同研究を行えたことを感謝します．前園 涼 教授（JAIST），本郷 研太 准教授（JAIST）には，QMC の基礎をご教授頂くと共に，QMC の方法論開発に携わるチャンスを与えて頂いたことに大変感謝しております．東後 篤史グループリーダー（NIMS）には，日々の議論に感謝しております．本寄稿に関連する研究は，科研 費（JP21K17752，JP21K03400），JSPS 海 外 特 別 研 究 員，MEXT 卓越研究員（JPMXS0320220025），NIMS 材料数値シミュレーター，HPCI 富岳利用研究課題（hp210038，hp220060，hp230030）の支援を受けて実施されました．文　　献1）  W. M. C. Foulkes et al., Rev. Mod. Phys., 73, 33（2001）.2）  Y. S. Al-Hamdani et al., Nat. Commun., 12, 3927（2021）.3）  X. Li et al., Phys. Rev. Res., 2, 033429（2020）.4）  J. Hermann et al., Nat. Chem., 12, 891-897（2020）.5）  K. Nakano et al., Phys. Rev. B, 109, 205151（2024）.6）  P. R. C. Kent et al., J. Chem. Phys., 152, 174105（2020）.7）  R. J. Needs et al., J. Chem. Phys., 152, 154106（2020）.8）  K. Nakano et al., J. Chem. Phys., 152, 204121（2020）.9）  https://trex-coe.github.io/champ-user-manual/10）  W. A. Wheeler et al., J. Chem. Phys., 158, 114801（2023）.11）  H. Niu et al., Phys. Rev. Lett., 130, 076102（2023）.12）  L. Monacelli et al., Nat. Phys, 19, 845-850（2023）.13）  Y. Nikaido et al., J. Phys. Chem. C, 126, 6000-6007（2022）.14）  B. X. Shi et al., J. Am. Chem. Soc., 145, 25372-25381（2023）.15）  T. Ichibha et al., Phys. Rev. Materials, 5, 064006（2021）.16）  A. Zen, Proc. Natl. Acad. Sci. U.S.A., 115, 1724-1729（2018）.17）  J. T. Krogel, Comput. Phys. Commun., 198, 154-168（2016）.18）  K. Nakano et al., J. Chem. Phys., 159, 224801（2023）.19）  R. Martin. Electronic structure: basic theory and practical methods（2020）.20）  J. Tiihonen et al., J. Chem. Phys., 154, 204111（2021）.21）  K. Nakano et al., J. Chem. Phys., 156, 034101（2022）.22）  S. Sorella et al., J. Chem. Phys., 133, 234111（2010）.23）  F. Jensen, Introduction to computational chemistry（2017）.24）  K. Nakano et al., Phys. Rev. B, 103, L121110（2021）.25）  A Tirelli et al., Phys. Rev. B, 106, L041105（2022）.26）  K. Momma and F. Izumi, J. Appl. Crystallogr., 44, 1272（2011）.であっても，行列式部分のパラメタ ciDet は変分最小ではないため（∂ ∂ ≠E ciDet/ 0），第3項を無視すると，誤った 力 の 計 算 が な さ れ る こ と に な る． な お，self-consistency 誤差は，DFT には存在しないことを強調したい．なぜなら，KS 方程式は，エネルギーに対する KS 軌道の変分が 0 になるように導出された方程式であるから 19），第 3 項は常に 0 である（図 1（a））．平面波基底では，上記第 1 項に相当する項のみで力が計算できる（i.e.，Hellmann–Feynman 力）19）．　この問題に対して，筆者は最近，第 3 項を量子モンテカルロ法で評価する方法を定式化し，QMC においても「正しい」力や圧力が計算できることを示した 5）．そ の 手 法 を， 自 身 が 開 発 す る TurboRVB お よ びTurboGenius に実装し，H2，Cl2 の原子に働く力の計算，および，c-BN の圧力計算に応用した 5）．なお，第一原理計算の別の実装である量子化学計算においても，変分最小化されていない波動関数を利用する際に同様の事態が生じることが知られており，同時に，そのバイアスを補正する手法も確立されている 23）．今回，筆者が考案した方法は，その量子化学計算で行われている手法と類似する手法である．詳細は原著論文を参照 5）されたい．図 1（b）に示したように，筆者の考案した方法により，VMC においても，正しい圧力が計算できることが示された．6. 今後の展望　筆者の提案した手法により，これまでは計算できなかった原子に働く力，および，圧力の計算が VMC で可能になったことは，今後の物質科学分野への応用に関して大きな進歩であると言える．筆者は，VMC による原子に働く力を利用して，フォノン分散を計算する研究 24）や，機械学習力場 25）を構築するような研究にも継続して取り組んでいる．しかしながら，現状の実装は，有限差分法を利用しており， 3NA 回の DFT計算が必要（NA は系に含まれる原子数）であることなど，実装上の課題は多く残っている．応用の観点ではこの点を改善していく必要があると考えている．　長期的には，VMC を超えて，DMC による力，圧力の計算手法確立が必須であろう．例えば，真にQMC がその威力を発揮するファンデルワールス力が支配的な物質などでは，VMC はその結合状態を描けないこともある．物質科学分野で興味が持たれるような複雑な物質に対する QMC の応用を進めていくためには，VMC を超えて DMC による力の計算を確立する必要があり，現在，計算手法の開発を進めている．筆 者 紹 介中野　晃佑（なかの　こうすけ）　2008 年， 大阪府立三国丘高等学校を卒業．2012 年，京都大学 工学部 工業化学科を卒業．学士（工学）．2014 年，京都大学大学院 工学研究科 物質エネルギー化学専攻を修了．修士（工学）．2017 年，北陸先端科学技術大学院大学 先端科学技術研究科 先端科学技術専攻（社会人）を修了．博士（情報科学）．旭硝子株式会社（現：AGC）常勤研究員，Scuola Internazionale Superiore di Studi Avanzati（SISSA）博士研究員，および，JSPS 海外特別研究員，北陸先端科学技術大学院大学 助教を経て，2023 年より，物質・材料研究機構（NIMS）に独立研究者（MEXT 卓越研究員）として勤務．専門： 第一原理量子モンテカルロ法，ソフトウェア開発．［連絡先］　〒 305-0047　茨城県つくば市千現1-2-1　物質・材料研究機構（NIMS）E-mail：NAKANO.Kosuke@nims.go.jp