第一原理分子動力学コード FEMTECK における 計算技法の紹介

1. 第一原理分子動力学コード FEMTECK における計算技法の紹介 産業技術総合研究所　計算科学研究部門 土田 英二

2. 構成 • 第一原理計算とは　　　　　　　～１０分 • FEMTECK の数値計算技法　 ～２５分 • 有限要素法 • 一般化固有値問題 • 多重格子法 • 最近の展開　　　　　　　　　　　～１０分 • T2K 上でのベンチマーク • 混合精度計算

3. Part I:第一原理計算とは

4. 第一原理計算の紹介 • 量子力学に基き、原子・分子から凝縮系まで、様々な系を同じように扱うことができる • ミクロな平均構造や振動スペクトルを計算することができる • 実験と相補的な役割を果たす • 古典的な力場で扱い難い系で威力を発揮 • 分極などの多体効果も自動的に考慮 • 電解質系のシミュレーションに最適の手法

5. フラーレン（ C60）の構造

6. フラーレンの電子密度

7. FEMTECK とは • Finite Element Method based Total Energy Calculation Kit の略称 • １０年程前から開発を続けている • 現在マニュアルを準備中 • 内部公開を目指している段階 • 解説記事：東大情報基盤センターニュース（2009年2月）

8. FEMTECK の概略 • ソースは 800 Kb （約３万行） • 前後処理はＣ言語 （上のサイズには含まず） • カーネル部分は全てFORTRAN • 必要なライブラリ • BLAS/LAPACK/ScaLAPACK • MPI • フラットＭＰＩで並列化されている • 最大 4096 並列まで動作確認済（T2K＠東大）

9. 応用計算 • 大規模な不規則系が中心 • 溶液系、界面など • オーダーＮ法にも対応している • 燃料電池やリチウム電池が重要なテーマ • プロトンジャンプが起きる • 古典ＭＤでは扱うことが難しい

10. リン酸 （Phosphoric Acid） • H3PO4： • 肥料や食品添加物、燃料電池の電解質等として広く使われている • E.Tsuchida, J. Phys. Soc. Jpn. 75, 054801 (2006).

11. シミュレーションの条件 • リン酸 ４８分子 • 原子数：384、軌道数：768 • 要素数：64000 (=40^3) • 総変数～４億個 • 周期境界条件 • 日立 SR8000 を６４ノード使って計算した • 理論値:500 GFlops、実測値:300 GFlops • 分子動力学～8000 ステップ（１０ピコ秒）

12. 分子動力学シミュレーション

13. 液相の動径分布関数

14. 振動スペクトル

15. NAFION (I) • 燃料電池の電解質として広く使われている • 水和のレベルを変えてシミュレーションを行い、構造やプロトンのダイナミクスを比較 • Y-K.Choe, E.Tsuchida, T.Ikeshoji, S.Yamakawa, and S.Hyodo, J.Phys.Chem.B 112, 11586 (2008). • 産総研と豊田中研の共同研究（CREST）

16. NAFION (II) λ＝４．１ λ＝１２．７

17. 計算の詳細 • ナフィオン２本＋水４８分子 • 原子数：424、波動関数の数：1022 • 基底関数の数～５０万個 • Opteron クラスタ×１２８プロセッサを使用 • 理論ピーク性能 = ５１２ Gflops • 約９０℃でプロダクション・ラン • １ステップ毎に約５億個の未知数について最適化 • １ステップ当たり５分程度かかる • 20000 ステップ程度→二ヶ月以上かかる

18. 液体エタノール（C2H5OH） • オーダーＮ法の実証計算 • 1125 原子、１０ピコ秒 • ５倍弱の高速化を実現 • J.Phys.:Condens. Mat. 20, 294212 (2008).

19. その他の応用例 • 水素結合性液体 • 水・メタノール・ホルムアミド・硫酸水溶液等 • リチウム伝導体 • 炭化水素系高分子水溶液 • 半導体の表面構造 • 水・シリコン界面の解離吸着反応 • 大雑把に言って 15～20年前に古典分子動力学法で扱っていた程度の大きさの系を第一原理計算で調べることが可能

20. Part II:FEMTECK における数値計算技法の概略

21. Numerical Recipes • 全１７章から成る数値計算の標準的教科書 • FEMTECK では、前後処理まで含めると、ほぼ全ての章を参考にしている • 第一原理計算は非常に幅広い数値計算の知識が必要とされる

22. 主要な計算手法 • 有限要素法 • 基底関数、座標変換 • 一般化固有値問題の解法 • レイリー商の最適化問題に変換 • 準ニュートン法（省メモリ版） • ポアソン方程式の解法 • マルチグリッド法、共役勾配法

23. Kohn-Sham 方程式 (I) ハミルトニアン 電子密度

24. Kohn-Sham 方程式 (II) • 基本的には、一電子に対するシュレーディンガー方程式と同形 • 一種の固有値問題に相当 • 以下では全て「周期境界条件」を仮定 • 原子の位置や他の電子の影響をポテンシャルとして取り込んでいる　　　　　　　　　　　→　ハミルトニアンはψにも依存する

25. 計算の手順 １：原子配置を決める ２：ハミルトニアンを計算 ３：ハミルトニアンに対応する波動関数を計算 ４：波動関数を使って原子に働く力を計算する ５：力に応じて原子を動かす→１へ

26. 基底関数 (I) • Kohn-Sham 方程式を解くためには、波動関数ψを基底関数で展開する必要がある • ハミルトニアンは「行列」になる • 従来の基底関数： • 原子基底（スレーター/ガウス関数） • 平面波 （≒スペクトル法） • いずれも物理的な意味を持つ

27. 基底関数 (II) • 原子基底 • 主に孤立境界条件 • 少数の原子からなる「分子」の計算から発展 • 平面波 • 主に周期境界条件 • 周期性が本質であるような「結晶」から発展 • 現在、大規模計算のターゲットは「乱雑系」 • 原子基底や平面波が必ずしも最適ではない

28. 基底関数 (III) • ミクロ系の場合も、三次元で偏微分方程式を解くことに変わりはない • マクロ系の標準解法である差分法・有限要素法が使えるはず • 近年、これらの手法が急速に適用され始めた • J. Chelikowski, 1994～（差分法） • E. Tsuchida and M. Tsukada, 1995～（有限要素法） • T. L. Beck, Reviews of Modern Physics 72, 1041 (2000).

29. 有限要素法の特長 (I) • 解像度を局所的に変えられる • 局所分割や座標変換を利用 • 基底関数の数を大幅に減らせる • 精度を体系的に改善できる • 「機械的に」精度を上げることができる • 基底関数同士が疎である • 計算量は基底関数の数に比例（⇔ 原子基底） • 基底関数は常に線形独立

30. 有限要素法の特長 (II) • オーダーＮ法に向いている • 本来は原子数の三乗で計算量が増える（メモリは二乗） • 多少の近似の下で、両方ともほぼ一乗にすることができる • 高い並列性を持つ • 最新のスパコンは全て大規模並列 • 通信の大半は局所的で済む

31. 一次の有限要素展開 • 各節点での値のみで決まる

32. 三次の有限要素展開 • 各節点での値と微分で決まる

33. 三次の有限要素基底 各節点においた基底関数の線形結合で波動関数を近似する

34. 行列の特徴 • ハミルトニアン： • 三次元の場合、基底関数は最近傍の基底関数のみと重なりを持つ→H0 は 8*33 =216個の非零成分を持つ疎行列 • 実際には密行列のブロックに分解する→他の行列との演算はほぼ BLAS3 を使える

35. メッシュの取り方 (I) • 各原子の近傍では波動関数が激しく振動する • 全空間を一様に分割するのは効率が悪い • 何らかの方法で原子付近だけ解像度を上げたい • 局所分割 • 座標変換 • ...

36. メッシュの取り方 (II):局所分割 • 原子が動いた場合の扱いが難しい • 並列計算機上での振舞いはあまり芳しくない • 解像度を極端に上げる時には有力

37. メッシュの取り方 (III):曲線座標系 • 原子の近くでメッシュが密になるよう座標変換 • 原子が複数あるときは線形に足合わせる • 原子が移動すると　座標系も追随する • 並列化し易い

38. メッシュの取り方 (IV): 座標変換 • 座標変換の式： • x→ x は直接計算できないので、Newton-Raphson 法を使って反復して求める • 極端に adaptation が強いと変換不能になる • ２倍強程度ならば特に問題はない

39. 電子状態計算の方法 (I) • Kohn-Sham 方程式を有限要素法で展開： • H, Fは巨大な疎行列 • バンド幅は２００程度 • 有限要素法の場合は基底関数が非直交なため、一般化固有値問題になる

40. 電子状態計算の方法 (II) • ハミルトニアンの次元をNとすると、N/100～N/1000 程度の固有値を（下から）求める必要がある • N は１００万を超えることもある • 求めたい固有値の上限と、その次の固有値の間には通常ギャップがある • 「半導体」「絶縁体」に相当 • 「金属」ではギャップがない→ここでは考えない

41. 電子状態計算の方法 (III) • 直接ハミルトニアンを変形するのはほぼ不可能なため、反復法で計算する • 固有値問題を最適化問題に変換して解く • 求めたい固有値が一つの場合： 　をψについて最適化すると、最小の固有値＋固有ベクトルが求まる レイリー商

42. 電子状態計算の方法 (IV) • 複数の固有値を同時に求めたい場合： 　をψについて最適化すれば良い。ただし

43. 最適化問題に変換するメリット • 「一般化」固有値問題であってもほとんど手間が変わらない • 特に求めたい固有値群と、その次の固有値の間にギャップがある場合には非常にうまく行く • 全ての固有ベクトルを同時に求めるので、ブロック化が容易→ Flops が高い • Hが非線形であることを考慮し易い • 第一原理計算特有の事情

44. 最適化の方法 • 準ニュートン法の省メモリ版（Limited memory BFGS method）を使用 • Liu and Nocedal, Math. Prog. 45, 503 (1989). • 国内ではマイナーだが、６００回以上引用 • 通常の準ニュートン法の収束性を保ちつつ、メモリの使用量を大幅に削れる • 非線形問題に関しては共役勾配法よりも高速 • １反復毎に目的関数を一回、微分を一回、前処理一回計算すれば良い（一次元の最適化がほぼ不要）

45. 準ニュートン法の収束

46. 波動関数の初期値 • 最初のＭＤステップ • 粗いメッシュで乱数から解いた結果を利用（1-way multigrid） • ２ステップ以降 • 直近２ステップ程度の波動関数から「外挿」 • 簡単な割に効果が大きい

47. ポアソン方程式の高速解法 (I) • 電子状態が更新される度に、電子密度が作るポテンシャルを再計算する必要がある • 周期境界条件の下でポアソン方程式を解くことで計算できる

48. ポアソン方程式の高速解法 (II) • 曲線座標系を使っているため、ＦＦＴでは済まない • 大規模な線形方程式に相当する • Ax=b の形 • A は疎行列 / x, b は通常のベクトル • 直接解法は現実的ではない • 共役勾配法のような反復法を使う必要がある

49. ポアソン方程式の高速解法 (III) • 共役勾配法は系が大きくなると収束率が急激に悪化する • 多重格子法（マルチグリッド）を用いた解法 • 系のサイズによらず、一定の収束率を持つ • 単独で使うよりも、共役勾配法の前処理として利用した方が効率が良い

50. マルチグリッド法 (I) • 標準的な「Ｖサイクル」を使用 • S=Smoother, E=Exact Solver • 色々な波長の残差を同じ比率で減らせる Ｓ Ｓ Ｓ Ｓ Ｅ