量子化学における 超大規模半正定値計画問題と 並列計算による高速求解

1. 量子化学における超大規模半正定値計画問題と並列計算による高速求解 神奈川大学工学部情報システム創成学科 山下真 Makoto.Yamashita@is.kanagawa-u.ac.jp 東京工業大学 Mituhiro Fukuda 東京大学 　　　　　　 中田真秀 RAMP 18th 2006/10/13

2. Outline • 量子化学の基本的で重要な問題に • 基底状態のエネルギー(Ground State Energy) • そのときの電子構造 を求めることがある • この問題は、半正定値計画問題(SemiDefinite Programming, SDP)に緩和できる • 並列計算を用いて実用時間内に解くことができる • 量子化学＋数理最適化＋並列計算の３つの分野

3. 目次 • 量子論とSDP • 基底状態におけるエネルギーとは • SDPへの帰着 • 並列計算による高速化 • 数値実験

4. 量子論とSDP(その１：座標回転） • まずは、大まかな話から • 波動関数と物理量のオペレータの行列表現をそれぞれ　　　　とすると、物理量は　　　　　　　の２次形式で計算できる • 適切な回転行列　　を用いると対角行列　　で • 本質的な量（固有値）と座標回転を同時に考えることが必要

5. 量子論とSDP（その2：対称行列） • 本質的な量（固有値）と座標の回転を同時に考えることが必要 • 対称行列で変数の表現 • 電子構造を対称行列で表現すると、（直観的には）固有値が電子の個数となる • 変数が半正定値行列に制約される • 半正定値行列を変数に持つ線形計画問題が半正定値計画問題

6. 半正定値計画問題（SemiDefinite Programming, SDP)

7. 基底状態におけるエネルギーとは • 直観的には、原子・分子系の最も安定した状態が基底状態 • 化学反応で発生するエネルギーなどを計算するときに必要となるエネルギー • これを求めることは、量子化学における基本的なこと

8. 量子力学における定式化 • Heisenberg の行列力学 • Schrödinger の波動方程式 • これらは、Hilbert 空間の内積が一致するため、物理学的には等価 時間座標と空間座標を変数分離

9. 時間に依存しないSchrödinger方程式 • . • 最小固有値がＧＳＥ • 水素原子のときには • 微分方程式を解かずに最小固有値のみを求める

10. 行列表現 • 水素原子の場合、全空間での電子の数は 1 であるため、波動関数は Hilbert 空間 　　　に属している • Hilbert 空間での基底をとすると、波動関数を展開できる • このとき波動関数　　とベクトル　　を同一視できる

11. 行列の最小固有値 • Schrödinger 方程式を行列の固有方程式にただし、 • したがって、(Ritz の変分原理ともよばれている）

12. 多電子ハミルトニアン • ハミルトニアンに電子-電子間などの項が必要 • 原子核の位置を固定して近似(Born-Oppenheimer approximation) • . constant

13. 多電子波動関数の基底 • 多電子波動関数はPauliの排他原理を満たす必要がある • 一電子波動関数の基底をとすると、 電子波動関数の基底全体はである。ただし、

14. 既存手法 • 一般の多電子波動関数はと展開される • GSE計算の既存手法 • Hartree-Fock • Full CI • CCSD,CCSD(T), etc

15. Hartree-Fock • Self Consistent Field 法で得られた一電子波動関数 • エネルギーの低い順に • 基底状態の波動関数を　　　　　　　　　　　　　　　　とする • 使用しているのは、のひとつの項のみ • GSEと誤差１％程度とされている

16. Full CI • 一次結合に現れるすべての項を用いて最小固有値計算 • 与えられた基底については、厳密解を与える • 計算量が膨大　　　　　　　　　の行列の固有値分解が必要 • 大きな分子に適用できない

17. CCSD(T) • Coupled Cluster expansion using Single and Double excitations with perturbation treatment of Triples • 一次結合　　　　　　　　　　　　　　の中で効果の高い項のみで計算 • 量子化学の分野で信頼されている手法 • これと同等の精度を出すことがひとつの目安

18. SDPによるGSE計算 • Full CI 同様に最小固有値の計算をベース • Full CI では、von Neumann 密度行列と同等 • 必要な情報量は、２次の縮約密度行列のみ • ２次の縮約密度行列を変数とすることで計算量を抑える

19. SDPによるGSE計算の歴史（１） • 1960年代：Coleman, Garrodら • 半無限線形計画問題への帰着 • ＳＤＰ緩和 • Erdahl による独自のアルゴリズムでの数値実験 • SDP緩和が不十分なため良い結果を得られず • 1990年代：数理最適化の分野で主双対内点法の研究が活発に

20. SDPによるGSE計算の歴史（２） • 2001年：中田ら • 主双対内点法のソフトウェアSDPA [SemiDefinite Programming Alogorithm]を用いて計算 • P,Q,G条件などで良好な精度 • 2004年:Zhao ら • T1,T2条件を導入して精度向上 • 並列ソフトウェアSDPARA [SDPA paRAllel version]を使用

21. SDPへの帰着 • 主問題側　　　へと帰着させる方法と双対問題側　　　へと帰着させる方法がある • 計算効率は　　　の方がいいが、まずは分かりやすい へと帰着させる • ハミルトニアンも電子間に相互作用がないところから帰着

22. ＧＳＥ計算と１次縮約密度行列 • . • ここで、行列 をで定義すると 新しい変数 離散化１次縮約密度行列

23. N-representability • 密度行列はN-representabilityとよばれる条件を満たす必要がある。１次なら • 直観的には、それぞれの波動関数上の電子の数を 0 から 1 に制約していることになる • まとめると、 SDPのひとつだが標準形ではない

24. 標準形への変形とSDPとしてのサイズ • 対角行列で標準形の SDP へと変化させる • このとき標準形は 2Kを基底関数の数とする は、それぞれが の２つのブロックに分解

25. 双対問題側への定式化による効果 Primal Dual 少ないメモリで計算できる

26. ２次縮約密度行列 • 電子間相互作用に関するハミルトニアンを組み込む • この項を計算するために 2-RDM Γ を導入 • Γが波動関数に対応するための必要十分条件もN-representabilityである • 組み合わせ的な要素を含んでいるため、全部の条件を列挙して計算することは不可能

27. P,Q,G条件の導入 • P,Q,G条件はN-representabilityの中でも効率的で効果的にΓの領域を表現できる • それぞれ、次の行列の半正定値条件に対応 • P,Q,G行列はΓからの線形変換で得られる • これらをP,Q,Gの半正定値条件を SDP の制約条件に追加する（SDP緩和が得られる）

28. GSE計算のためのSDP • 2-RDMは物理量を計算するのに十分な情報を持っている 3-RDMの知識から得られる2-RDMの制約を追加 PQG T1T2 計算時間と数値精度に影響する

29. 2-RDMがあるときのSDPのサイズ • 2K=28のときには、単一CPU上のSDPA で8.6GB のメモリを使用する • 並列計算へと移行する

30. SDPARA http://homepage.mac.com/klabtitech/sdpa-homepage/ • SDPA parallel version(generic SDP solver) • 主双対内点法 • PCクラスタ上の分散メモリで超大規模なSDPを格納 • 並列計算による計算時間の短縮 • 実装に使用しているライブラリ • MPICH　　　　 （通信用ライブラリ） • ScaLAPACK　（並列版 LAPACK） • マルチプロセッサ化にともなって、普通のPCでも並列計算が標準の時代に

31. SDPの標準形

32. KKT 条件 • 最適解は KKT 条件で特徴付けられる • 主双対内点法では、行列の正定値性を維持しながら反復的に最適解を求める

33. 主双対内点法のフレームワーク

34. 探索方向の計算 2.CHOLESKY Schur complement matrix ⇒ Cholesky Factorizaiton Schur complement equation 1.ELEMENTS

35. 単一プロセッサ上でのボトルネック 計算時間の単位は秒Opteron 246 (2.0GHz) SDPARAでボトルネックを並列化

36. SDPAにおける疎性を活用した計算 • 行方向に計算式を変えて計算 • これを維持したまま並列計算へ F1 F2 F3

37. Row-wise distribution for evaluation of the Schur complement matrix • CPUが4 個の場合を表示 • それぞれのCPUは　行方向を担当 • . • CPU間の通信は一切なし • 効率的なメモリ管理

38. 並列 Cholesky 分解 • Schur complement matrix の Cholesky 分解には ScaLAPACK のルーチンを使用 • 行方向のメモリ管理から２次元巡回配置に分散メモリ上で再配置して効率的に 再配置

39. 分子 LiOH に関する計算時間 AIST super clusterOpteron 246 (2.0GHz) 6GB memory/node

40. 分子 HF の計算時間のスケーラビリティー ELEMENTS 63 times CHOLESKY 39 times Total 29 times ELEMENTSの行方向分割が非常に強力

41. 計算時間 （単位は秒） • 一番大きな分子 (NH) にも並列計算は効果 • 単一プロセッサの場合には、メモリスペースの関係で解くことができない • 22 日 (単一プロセッサ上の計算時間の予測値) → 7 時間に短縮

42. 数値精度比較 • 得られたGSEの精度を Full CIからの差で表現する • 既存手法 • HF (Hatree-Fock) • CCSD(T) (Coupled Cluster expansion using Single and Double excitations with perturbation treatment of Triples) • SDP 緩和 (PQG, T1T2)

43. 数値精度 • PQG は HF よりも精度がよい • T1T2は CCSD(T) とほぼ同等の精度 • T1T2は化学的精度の 1kcal/mol より高い精度 in Hartree (1 Hatree = 4.36x10^(-18) Joule)

44. まとめと今後の展開 • ＧＳＥ計算はSDP緩和に帰着できる • SDP緩和はGSEを高い精度で計算できる • メモリ消費や計算時間から、並列計算をすることが必要 • Hubbard モデルへの適用 • 双極子モーメントの計算などなど • 数値的安定性の向上 • メモリ消費の抑制

45. ご清聴ありがとうございました Ｔａｎｔｅ　Ｇｒａｚｉｅ