2022年度 第1回計算科学フォーラム
2022.10.4

大規模量子化学計算法
divide-and-conquer法の現在と展望
小林 正人
北海道大学大学院理学研究院
北海道大学WPI-ICReDD

本日の講演内容
 イントロダクション 量子化学計算とは
 なぜ大規模系を計算したいのか？

 分割統治(Divide-and-Conquer)法
 DC法とその他の大規模量子化学計算法
 DC法の理論とバッファ領域の自動決定 （DC法の応用）

 量子計算・新概念計算による量子化学計算
 Annealing + Bayesian-optimization CI (ABCI)法 まとめ2

3

量子化学計算とは 化学（物質科学）の対象 ＝ 分子（集団）・固体
 構成要素原子核電子(原子核と)電子の量子論的な振舞い

物質のあらゆる性質数値的に求める = 量子化学計算

 典型的量子化学計算(ab initio)
 入力: 原子の座標、電子数、計算手法・基底関数
 出力: 分子のエネルギー、電子分布（分子軌道）、プロパティ、エネルギー微分（場合によって）

量子化学計算のジレンマと近似計算法
 量子化学計算のジレンマ
ディラック(1929)[1]『化学を支配する法則は完全にわかった。 しかし、複雑すぎて解けない方程式に行き着くのが問題だ』

 解くべき方程式: Full CI （完全にわかっている）|Ψ⟩𝑐⋯|𝑛 𝑛 ⋯ 𝑛 ⟩• |𝑛 𝑛 ⋯ 𝑛 ⟩: 電子配置【ni=1（占有） or 0（非占有）】∑ 𝑛𝑁

 古典計算機を利用した計算: O(Ne!)
様々な近似計算法の開発により超克= 量子化学計算の歴史[1] P.A.M. Dirac, Proc. Roy. Soc. A 123, 714 (1929).

45

様々な量子化学計算手法と計算時間・精度
 量子化学計算にかかる計算時間と精度

分子の大きさの3乗に比例して計算時間増大
方法

Hartree-Fock(HF)法

密度汎関数理論(DFT)

MP2法
(摂動法)

CCSD法

CCSD(T)法

計算時間O(N3)

O(N3)

O(N5)

O(N6)

O(N7)

計算精度

定性的

計算時間は
 精度の低い理論でもO(N3) 精度が上がるにつれて莫大に正確

7

量子化学計算にかかる時間と精度
 様々な量子化学計算手法と精度分子の大きさの3乗に比例して計算時間増大
方法

Hartree-Fock(HF)法

密度汎関数理論(DFT)

MP2法
(摂動法)

CCSD法

CCSD(T)法

計算時間

O(N3)

O(N3)

O(N5)

O(N6)

O(N7)

計算精度

定性的

1000倍性能の計算機で×10

正確

×10

計算時間は
 精度の低い理論でもO(N3) 精度が上がるにつれて莫大に×4.0

×3.2×2.7『富岳』をただ使うだけでは大きな分子を扱えない最先端の分子理論も不可欠量子化学計算の種類
 波動関数理論: N電子波動関数
 Schrödinger方程式に依拠
 Hartree-Fock (HF)法による1電子の波動関数から構築• 平均的なポテンシャル中の1電子描像

 HF波動関数を系統的に改善する電子相関計算法• 摂動論: MP2, MP3, MP4, …
• 結合クラスター（CC）法: CCSD, CCSD(T), CCSDT, …

 密度汎関数理論: 電子密度の空間分布
 Hohenberg-Kohnの定理に依拠• 厳密なエネルギーは電子密度分布関数𝜌 𝒓 の汎関数
• HF法とほぼ同じ手続きで電子相関を考慮可能• 厳密なエネルギー汎関数は今のところ不明
 様々な近似汎関数: e.g.) B3LYP, M06, ωB97X-D, …

8

9

量子化学計算の手続き
Atomic coordinates, # electrons,
method, basis set

N: system size (# atoms)

increases withcubic of size

Guess electron density DConstruct Hamiltonian F from D・・・ O(N2-4) time

Diagonalize F to get MOs {ψp}・・・ O(N3) time 𝜓Construct density D from MOs
NoD converged?
Yes
HF or DFT?
NoAccurate post-HF calculation
End

𝐷𝐶 𝐶

Self-consistent field (SCF)
Yes

End・・・ O(N>5) time𝐶 𝜒

大規模量子化学計算手法
 並列化とプログラムの工夫で頑張る ProteinDF （タンパク質密度汎関数プログラム）
 RSDFT （実空間密度汎関数理論）

 密度行列を近似計算
 CONQUEST 数学・アルゴリズムにより高速化
 積分のRI計算, Laplace変換MP2法[1] D.G. Fedorov and K. Kitaura eds., The Fragment Molecular Orbital Method(CRC Press, 2009);[2] A. Imamura, Y. Aoki, and K. Maekawa, J. Chem. Phys. 95, 5419 (1991).

10

『方法序説』 に学べ ルネ・デカルトの名著 (1637)
 「近代哲学の父」
 Cogito ergo sum （我思う、故に我あり）
 デカルト座標系の提案

 物事を考える4つの原則のひとつ
 「検討する困難（大きな問題）を問題をよりよく解くために必要な小部分（小さな問題）に分割せよ。」

世界で最初の「分割統治」という考え方

11

12

大規模量子化学計算手法
 並列化とプログラムの工夫で頑張る ProteinDF （タンパク質密度汎関数プログラム）
 RSDFT （実空間密度汎関数理論）

 密度行列を近似計算
 CONQUEST 数学・アルゴリズムにより高速化
 積分のRI計算, Laplace変換MP2法

 系を分割して計算を簡略化

O(N)へ

 フラグメント分子軌道（FMO）法[1] エロンゲーション法[2] 分割統治（DC）法[1] D.G. Fedorov and K. Kitaura eds., The Fragment Molecular Orbital Method(CRC Press, 2009);[2] A. Imamura, Y. Aoki, and K. Maekawa, J. Chem. Phys. 95, 5419 (1991).

13フラグメント分子軌道(FMO)法
 分子を（単結合で）切断
 結合に使われる電子は混成軌道を使って一方に寄せる（HOP）

 エネルギー等のプロパティ
は多体展開で求める
EFMO1   EI (FMO1)
I

EFMO2   EI   ( EIJ  EI  EJ ) (FMO2)I

I J

EI: モノマーのエネルギー
EIJ: ダイマーのエネルギー
EIJK: トリマーのエネルギー

・・・

EFMO3  EFMO2   ( EIJK  EIJ  EJK  EIK  EI  EJ  EK ) (FMO3)
I J K 多体展開を打ち切り計算の大幅な高速化[1] D.G. Fedorov and K. Kitaura eds., The Fragment Molecular Orbital Method(CRC Press, 2009).

分割統治(DC)量子化学計算

14 部分系に分けて計算するのはFMO法と同じ
 部分系の周辺をバッファ領域として計算に含め、環境とのあらわな相互作用を考慮
 部分系の電子数はフェルミ準位を導入して自動決定(予め指定する必要なし)
SCFと電子相関計算で取り扱いが大きく異なる
 分割統治SCF法 全系の密度行列を部分系の寄与の和で表現
 分割行列を導入

 分割統治post-HF電子相関計算
 部分系の相関エネルギーを部分系の分子軌道で表現

15

分割統治量子化学計算のスキーム[1,2]
中央領域
(部分系)

局在化領域

バッファ領域

divide

部分系ごとの方程式を解く

部分系の分子軌道共通のフェルミ準位を設定

全電子数を保存conquer

バッファ領域が精度に大きく影響

SCF
全系の電子密度・エネルギー[1] MK and H. Nakai, in Linear-Scaling Techniques in Computational Chemistry
and Physics (2011), pp.97-127.
[2] W. Yang and T.-S. Lee, J. Chem. Phys. 103, 5674 (1995).

16

DC-HF計算: CPU時間
 最初のSCF cycleにかかるCPU時間[1]基底関数: 6-31G**
中央領域: C2H2(3) (1 unit)
バッファ: 6 units
Pentium4 / 3.20 GHz

 ポリエン鎖 CnHn+2
3000

Fock⾏列対⾓化

600

O(n )

400
200

O(n1.3)003060n90120

Fock⾏列⽣成(w/ FMM)

2500

Conventional HF
DC-HF
3.2

CPU time [sec]CPU time [sec]

800

150

Conventional HFDC-HF

2000

O(n1.6)

1500
1000

O(n1.1)500
003060n90

DC法によりほぼリニアスケーリングを達成[1] MK and H. Nakai, in Linear-Scaling Techniques in Computational Chemistryand Physics (2011), pp.97-127.

120

150

DC-SCF計算に対するバッファ領域自動決定法

17 2 layerのバッファ領域を利用して、DC法の精度のカギとなるバッファ領域を自動決定する方法を開発
i.

外側バッファ領域に属する原子が与えているエネルギー
誤差Δ𝐸 を摂動的に見積もり
Δ𝐸

2Δ𝐷 𝐹
∈

∈𝑺

外側バッファ
内側バッファ

中央[1] MK, T. Fujimori, and T. Taketsugu, J. Comput. Chem. 39, 909 (2018).

DC-SCF計算に対するバッファ領域自動決定法

18

 2 layerのバッファ領域を利用して、DC法の精度のカギとなるバッファ領域を自動決定する方法を開発
i.

外側バッファ領域に属する原子が与えているエネルギー
誤差Δ𝐸 を摂動的に見積もり
ii. 外側バッファ領域に属する原子を内側バッファ領域へ
iii. Δ𝐸 の大きい原子を基準に外側バッファ領域を再構築
外側バッファ内側バッファ

中央

中央[1] MK, T. Fujimori, and T. Taketsugu, J. Comput. Chem. 39, 909 (2018).

中央・・・

19

バッファ領域自動制御DC法の精度と計算時間
 誤差と計算時間の初期バッファサイズ依存性
 系: 水分子1000個の箱形モデル系
 計算レベル: PM3
3.5
4.5
4.0
5.0
4.5
5.5
5.0
6.0
5.5
6.5
Standard PM3

閾値: 0.1 μEh

Energy/Eh(Diff./μEh•atom−1)

Time/s−11945.190938−11945.190942−11945.190837−11945.190719−11945.190414−11945.192376(+0.48)(+0.48)(+0.51)(+0.55)(+0.65)

250
246
233
209
2092443

1原子当たりのエネルギー誤差を約0.5 μEh程度に制御
 計算時間を大幅に低減
[1] MK, T. Fujimori, and T. Taketsugu, J. Comput. Chem. 39, 909 (2018).

20

分割統治post-HF電子相関計算[1] DC-HF計算 → 部分系MO
 Local canonical MO 占有・非占有軌道ともバッファ領域に広がっている

中央領域

 部分系の相関エネルギー 部分系のMOから相関エネルギーを計算
 中央領域の相関エネルギーのみが必要
中央領域

エネルギー密度解析(EDA)[2]の利用[1] MK, Y. Imamura, and H. Nakai, J. Chem. Phys. 127, 074103 (2007).
[2] H. Nakai, Chem. Phys. Lett. 363, 73 (2002).

バッファ領域

局在化領域

21分割統治post-HF電子相関計算
 相関エネルギー（Nesbetの定理）
occ vir

局在化領域

E   (ia | jb)[2tia , jb  tib , ja ]i , j a ,b

 部分系のΔEを部分系の軌道から計算
occ( ) vir( )E  


i, j中央領域

 a ,b
 


C(a|jb)[2tt iia , jb
ib , ja ]S ( )EDAS() : 部分系αの中央領域のAO  (ia|jb )t MP2の場合: ia , jb           abij 全相関エネルギー = 部分系の相関エネルギーの和E subsystemE DC-MP2法

22DC-MP2法の計算時間
 MP2計算にかかる時間[1]rb ＝ 6.0 Å
部分系: 1原子MP2/6-31G

 βストランドグリシンペプチド (Gly)nXeon (Paxville) 2.8 GHz (1CPU)

1800

O(n5.1)1600

CPU time [min]

1400
1200
1000
800

Canonical MP2
DC-MP2600
400O(n1.5)200005

10152025

計算時間: 系の大きさに対してほぼ線形 必要メモリ量: 系の大きさにほぼ非依存30n[1] MK, Y. Imamura, and H. Nakai, J. Chem. Phys. 127, 074103 (2007).

(最大部分系の大きさに依存)

DC-MP2計算に対するバッファ領域自動決定法

23

 MP2で取り込まれる動的電子相関エネルギーは、HFエネルギーよりも局所的 SCFよりも小さなバッファ領域を用いてDC電子相関計算を効率化[1] DC-SCF計算で自動決定したバッファ領域の各原子に対して、MP2相関エネルギーへの寄与を算定
 AO-Laplace MP2法[2]の考え方を利用
Δ𝐸𝑋∈𝑺𝑌𝐴2𝐴∈max 𝐀 /𝜔𝐴𝜇𝜈 𝜇𝜈

⁄X, Y: エネルギー重み付き密度行列[1] MK and H. Nakai, Int. J. Quantum Chem. 109, 2227 (2009).
[2] M. Häser, Theor. Chim. Acta 87, 147 (1993);
MK and H. Nakai, Chem. Phys. Lett. 420, 250 (2006).

DC-MP2計算に対するバッファ領域自動決定法
 自動制御型DC-MP2計算の手続き[1]i. バッファ領域自動制御型DC-HF計算を実行
ii. バッファ領域に属する各原子の電子相関エネルギーのおよその寄与Δ𝐸を算定
iii. Δ𝐸が小さい原子をバッファ領域から削除

中央

中央[1] T. Fujimori, MK, and T. Taketsugu, J. Comput. Chem. 42, 620 (2021).

24

25

バッファ領域拡大の様子
 自動制御型DC法によるバッファ領域の変化[1] Trp-cageのDC-HF/DC-MP2計算の場合
Initial

DC-HF

DC-MP2

異方性を持って拡大
部分系が適切なバッファ領域を選択[1] MK, T. Fujimori, and T. Taketsugu, J. Comput. Chem. Jpn. 20, 48 (2021).

26

自動制御型DC-MP2計算の精度
 自動制御型DC-MP2計算の初期半径依存性[1] 系: シニョリン
 計算レベル: MP2/6-31G(d)
𝑟/Å

/Å

𝑟

MP2 Energy/Eh(Diff.)/μEh atom−1<𝑙,/Å

><𝑙,> /Å

3.54.5-11.194860(+0.43)

8.248

7.606

4.05.0-11.194810(+0.79)

8.121

7.606

4.55.5-11.194825(+0.68)

8.174

7.582

5.06.0-11.194820(+0.71)

8.304

7.606

5.5

6.5-11.194814(+0.76)

8.151

7.606

Standard-11.194919

初期条件によらず、エネルギー誤差を一定に制御
 HF計算よりも小さい主軸半径を選択
[1] T. Fujimori, MK, and T. Taketsugu, J. Comput. Chem. 42, 620 (2021).

27

Density functional tight binding (DFTB)法[1] DFTB法: DFTをベースとした半経験的計算法
 DFTエネルギー
Tight binding近似： ρ = ρ0 + δρ
交換相関項に対して2次のTaylor展開: Exc[ρ0 + δρ]

密度に依存する項に対してモノポール近似︓ δρA ≈ ΔqAF00Y00
DFT計算から決定するパラメータ DFTBエネルギー
Mulliken charge

Charge independent

Charge dependent• Hamiltonian行列 𝐻𝐻

𝑆Repulsive energy

1𝛾2• Mulliken電荷を自己無撞着に決定 (SCC)[1] D. Porezag, Th. Frauenheim, Th. Köhler, G. Seifert, and R. Kaschner,
Phys. Rev. B 51, 12947 (1995).

𝛾2電子積分不要!

Δ𝑄𝜇 ∈ 𝐴, 𝜈 ∈ 𝐵

28

DC-DFTB計算のパフォーマンス
 計算時間のサイズ依存性[1] 水 Nwater 分子系
 計算機: Intel Xeon 1ノード (8コア)1ノードでも十分に高速（水2000分子が3分）
 高並列化も達成
[1] H. Nishizawa, Y. Nishimura, MK, S. Irle, and H. Nakai, J. Comput. Chem. 37,
1983 (2016).

量子化学計算と解析法
 量子化学計算で得られる結果
 エネルギー: 系全体で単一の量• 応答量（微分）もコストをかければ計算可双極子モーメント、勾配、分極率･･･

 波動関数: 大きな系になるほど複雑化 理解・予測・記述子開拓のためには解析法が重要 特にローカル（原子・残基・結合・成分）な量への射影• 電子密度解析• エネルギー解析

29

大規模量子化学計算と解析法

30

 フラグメント分割型大規模量子化学計算法
フラグメントに基づくローカルな量の算定が容易
 フラグメント分子軌道(FMO)法• IFIE: フラグメント間相互作用エネルギー解析• PIEDA: IFIE + 諸熊エネルギー分解

 DC法• 部分系エネルギー解析
 部分系の密度行列から部分系のエネルギーを算出• 電子相関エネルギー: 部分系に対してadditiveな表式

31HIV-1逆転写酵素(RT) ウィルスの一本鎖RNAを鋳型として二本鎖のプロウィルスDNAを合成する酵素（HIVの増殖に不可欠） 非核酸系逆転写酵素阻害剤(NNRTIs)
 低濃度で高選択的にHIV-1と結合
 MK-4965: 野生(WT)および主要変異ウィルスに対して高い阻害活性を示す新奇NNRTI
MK-4965

WTおよび変異HIV-1 RTに対する阻害活性[1]CompoundEfavirenz

IC50 (nM)
WT

K103N

Y181C

MK-4965

0.23 ± 0.09

0.40 ± 0.15

0.39 ± 0.16

Nevirapine

137 ± 5.89300 ± 500

8900 ± 1800

Efavirenz

0.42 ± 0.06

10.2 ± 6.6

0.31 ± 0.17

Etravirine

0.60 ± 0.12

1.3 ± 0.17

1.96 ± 0.9

Nevirapine[1] M.-T. Lai et al., Antimicrob. Agents Chemother. 53, 2424 (2009).

Etravirine

32HIV-1 RTとMK-4965の複合体
MK4965Tyrosine (Y)

野生型 (WT)MK4965

Cysteine(C)Y181C

DC計算によるMK-4965の相互作用解析

33

 HIV-1 RT (WT/Y181C)とMK-4965の間にはたらく相互作用エネルギー
 相互作用エネルギー: Eint = Ecomplex − EAA − EMK-4965
Tyr/Cys181との相互作用は小 Lys102及びLys103とSL3の間に大きな相互作用


SL3
SL2
SL1

データを高次元化して解析[1] P. Saparpakorn, MK, S. Hannongbua, and H. Nakai, Int. J. Quantum Chem.
113, 510 (2013).

本日の講演内容
 イントロダクション 量子化学計算とは
 なぜ大規模系を計算したいのか？

 分割統治(Divide-and-Conquer)法
 DC法とその他の大規模量子化学計算法
 DC法の理論とバッファ領域の自動決定 （DC法の応用）

 量子計算・新概念計算による量子化学計算
 Annealing + Bayesian-optimization CI (ABCI)法 まとめ34

量子化学計算のジレンマと近似計算法
 量子化学計算のジレンマ
ディラック(1929)[1]『化学を支配する法則は完全にわかった。 しかし、複雑すぎて解けない方程式に行き着くのが問題だ』

 解くべき方程式: Full CI （完全にわかっている）|Ψ⟩𝑐⋯|𝑛 𝑛 ⋯ 𝑛 ⟩• |𝑛 𝑛 ⋯ 𝑛 ⟩: 電子配置【ni=1（占有） or 0（非占有）】∑ 𝑛𝑁

 古典計算機を利用した計算: O(Ne!)
量子コンピュータや新概念計算機を利用して解く！
[1] P.A.M. Dirac, Proc. Roy. Soc. A 123, 714 (1929).

35

37

XBK変換法[1,2] CIハミルトニアンをQUBOハミルトニアンに変換
 第二量子化ハミルトニアン:
𝐻ℎ 𝑎 𝑎Γ𝑎 𝑎 𝑎 𝑎 スピンハミルトニアンに変換（Jordan-Wigner or Bravyi-Kitaev）
 r個のレプリカを使うことで、σzだけのハミルトニアンに変換
𝐻ℎ 𝜎ℎ,,𝜎 𝜎ℎ, ,, , 2次化 → QUBO

[1] R. Xia, T. Bian, and S. Kais, J. Phys. Chem. B 122, 3384 (2018).
[2] M. Streif, F. Neukart, and M. Leib, Lect. Notes Comput. Sci. 111, 11413 (2019).

𝜎 𝜎 𝜎

⋯38

XBK変換法[1,2]38

 CIハミルトニアン→QUBOハミルトニアン
軌道の数 M

符号+++−
全通りで試行1234

単一電子配置|𝑛 𝑛 ⋯ 𝑛 ⟩3Φ12−Φ34
レプリカ数rqubit

アニーリング計算で決定

 CI係数が大きい電子配置と小さい電子配置が混在
レプリカ数r → 大[1] R. Xia, T. Bian, and S. Kais, J. Phys. Chem. B 122, 3384 (2018).
[2] M. Streif, F. Neukart, and M. Leib, Lect. Notes Comput. Sci. 111, 11413 (2019).

39

ABCI（Annealing + Bayesian-opt. CI）法 各電子配置に重みビットを付加
軌道の数 M

重みビット−12

係数= +21234= +2= +15Φ12−Φ34
レプリカ数r= −1波動関数の自由度が拡大
重みビットの組み合わせをベイズ最適化で決定

40H2のポテンシャルエネルギー曲線
 重みビット数を変化させたときのPECの変化
 レプリカ数: r=4（固定=qubit数変化なし）
 重みビット数: 1（XBK法と同じ）, 2, 3, 4

重みビット数増加
→精度向上
6-31G
Active space: (2,3)

量子計算技術と大規模計算
 量子計算・新概念計算
 複雑な電子状態を計算するポテンシャル
 当面はqubit数に制限
困難は分割せよ

分割統治（DC）法のような分割型計算との融合が必須

41

まとめ：量子化学計算の大規模系への展開

42

 大規模系の量子化学計算を実現する
分割統治（DC）法を開発
 バッファ領域を用いて全系を重なりながら分割• バッファ領域を自動決定する方法も開発 電子相関計算（MP2, CCSD, CCSD(T)）にも対応
 DFTB法やxTB法と組み合わせた超大規模計算

 データ科学nativeの時代・量子計算の時代に向けて
 大規模計算からのデータ抽出が重要• 構造分類の自動化、相互作用解析

 量子計算機による量子化学計算にもDC法が必要

43

共同研究者と謝辞
 北海道大学
 武次 徹也 先生 赤間 知子 先生
 藤森 俊和 博士
 田代 啓介 君
 能條 小夜子 さん 米山 亮 君
 藤原 孝太郎 君

 研究費
 JSTさきがけ
 京都大学ESICB
 科研費
 ポスト京 重点課題5

 計算機
 早稲田大学
 京コンピュータ（理研）
 中井 浩巳 先生
 計算科学研究センター（岡崎） P. Saparpakorn 先生（Kasetsart大）
 九州大学 吉川 武司 先生（東邦大）
 ほか、両研究室のみなさま44

ご清聴ありがとうございました