Transcript ppt

第一原理分子動力学コード
FEMTECK における
計算技法の紹介
産業技術総合研究所 計算科学研究部門
土田 英二
構成
第一原理計算とは
～１０分
FEMTECK の数値計算技法 ～２５分
有限要素法
一般化固有値問題
多重格子法
最近の展開
T2K 上でのベンチマーク
混合精度計算
～１０分
Part I:
第一原理計算とは
第一原理計算の紹介
量子力学に基き、原子・分子から凝縮系ま
で、様々な系を同じように扱うことができる
ミクロな平均構造や振動スペクトルを計算
することができる
実験と相補的な役割を果たす
古典的な力場で扱い難い系で威力を発揮
分極などの多体効果も自動的に考慮
電解質系のシミュレーションに最適の手法
フラーレン（ C60 ）の構造
Bond length[A]
Bond angle[degrees]
Ｔｈｉｓ ｗｏｒｋ
1.39±0.003
1.45±0.002
108.0±0.17
120.0±0.22
Other
1.39
1.44
108
120
Exp.
1.40±0.015
1.45±0.015
フラーレンの電子密度
FEMTECK とは
Finite Element Method based Total
Energy Calculation Kit の略称
１０年程前から開発を続けている
現在マニュアルを準備中
内部公開を目指している段階
解説記事：東大情報基盤センターニュース
（2009年2月）
FEMTECK の概略
ソースは 800 Kb （約３万行）
前後処理はＣ言語 （上のサイズには含まず）
カーネル部分は全てFORTRAN
必要なライブラリ
BLAS/LAPACK/ScaLAPACK
MPI
フラットＭＰＩで並列化されている
最大 4096 並列まで動作確認済（T2K＠東大）
応用計算
大規模な不規則系が中心
溶液系、界面など
オーダーＮ法にも対応している
燃料電池やリチウム電池が重要なテーマ
プロトンジャンプが起きる
古典ＭＤでは扱うことが難しい
リン酸 （Phosphoric Acid）
H3PO4：
肥料や食品添加物、燃料電池の電解質等とし
て広く使われている
E.Tsuchida, J. Phys. Soc. Jpn. 75, 054801
(2006).
シミュレーションの条件
リン酸 ４８分子
原子数：384、軌道数：768
要素数：64000 (=40^3)
総変数～４億個
周期境界条件
日立 SR8000 を６４ノード使って計算した
理論値:500 GFlops、実測値:300 GFlops
分子動力学～8000 ステップ（１０ピコ秒）
分子動力学シミュレーション
液相の動径分布関数
振動スペクトル
NAFION (I)
燃料電池の電解質として広く使われている
水和のレベルを変えてシミュレーションを
行い、構造やプロトンのダイナミクスを比較
Y-K.Choe, E.Tsuchida, T.Ikeshoji,
S.Yamakawa, and S.Hyodo,
J.Phys.Chem.B 112, 11586 (2008).
産総研と豊田中研の共同研究（CREST）
NAFION (II)
λ＝４．１
λ＝１２．７
計算の詳細
ナフィオン２本＋水４８分子
原子数：424、波動関数の数：1022
基底関数の数～５０万個
Opteron クラスタ×１２８プロセッサを使用
理論ピーク性能 = ５１２ Gflops
約９０℃でプロダクション・ラン
１ステップ毎に約５億個の未知数について最適化
１ステップ当たり５分程度かかる
20000 ステップ程度→二ヶ月以上かかる
液体エタノール（C2H5OH）
オーダーＮ法の実証計算
1125 原子、１０ピコ秒
５倍弱の高速化を実現
J.Phys.:Condens. Mat. 20, 294212 (2008).
その他の応用例
水素結合性液体
水・メタノール・ホルムアミド・硫酸水溶液等
リチウム伝導体
炭化水素系高分子水溶液
半導体の表面構造
水・シリコン界面の解離吸着反応
大雑把に言って 15～20年前に古典分子動力
学法で扱っていた程度の大きさの系を第一原理
計算で調べることが可能
Part II:
FEMTECK における
数値計算技法の概略
Numerical Recipes
全１７章から成る数値
計算の標準的教科書
FEMTECK では、前後
処理まで含めると、ほ
ぼ全ての章を参考に
している
第一原理計算は非常
に幅広い数値計算の
知識が必要とされる
主要な計算手法
有限要素法
基底関数、座標変換
一般化固有値問題の解法
レイリー商の最適化問題に変換
準ニュートン法（省メモリ版）
ポアソン方程式の解法
マルチグリッド法、共役勾配法
Kohn-Sham 方程式 (I)
ハミルトニアン
電子密度
Kohn-Sham 方程式 (II)
基本的には、一電子に対するシュレーディ
ンガー方程式と同形
一種の固有値問題に相当
以下では全て「周期境界条件」を仮定
原子の位置や他の電子の影響をポテン
シャルとして取り込んでいる
→ ハミルトニアンはψにも依存する
計算の手順
１：原子配置を決める
２：ハミルトニアンを計算
３：ハミルトニアンに対応する波動関数を計算
４：波動関数を使って原子に働く力を計算する
５：力に応じて原子を動かす→１へ
基底関数 (I)
Kohn-Sham 方程式を解くためには、波動関
数ψを基底関数で展開する必要がある
ハミルトニアンは「行列」になる
従来の基底関数：
原子基底（スレーター/ガウス関数）
平面波 （≒スペクトル法）
いずれも物理的な意味を持つ
基底関数 (II)
原子基底
主に孤立境界条件
少数の原子からなる「分子」の計算から発展
平面波
主に周期境界条件
周期性が本質であるような「結晶」から発展
現在、大規模計算のターゲットは「乱雑系」
原子基底や平面波が必ずしも最適ではない
基底関数 (III)
ミクロ系の場合も、三次元で偏微分方程式を解くこ
とに変わりはない
マクロ系の標準解法である差分法・有限要素法が使える
はず
近年、これらの手法が急速に適用され始めた
J. Chelikowski, 1994～（差分法）
E. Tsuchida and M. Tsukada, 1995～（有限要素法）
T. L. Beck, Reviews of Modern Physics 72, 1041
(2000).
有限要素法の特長 (I)
解像度を局所的に変えられる
局所分割や座標変換を利用
基底関数の数を大幅に減らせる
精度を体系的に改善できる
「機械的に」精度を上げることができる
基底関数同士が疎である
計算量は基底関数の数に比例（⇔ 原子基
底）
基底関数は常に線形独立
有限要素法の特長 (II)
オーダーＮ法に向いている
本来は原子数の三乗で計算量が増える（メ
モリは二乗）
多少の近似の下で、両方ともほぼ一乗にす
ることができる
高い並列性を持つ
最新のスパコンは全て大規模並列
通信の大半は局所的で済む
一次の有限要素展開
各節点での値のみで決まる
三次の有限要素展開
各節点での値と微分で決まる
三次の有限要素基底
各節点においた基底関数の線形
結合で波動関数を近似する
行列の特徴
ハミルトニアン：
三次元の場合、基底関数は最近傍の基底
関数のみと重なりを持つ→H0 は 8*33
=216個の非零成分を持つ疎行列
実際には密行列のブロックに分解する→他の
行列との演算はほぼ BLAS3 を使える
メッシュの取り方 (I)
各原子の近傍では波動関数が激しく振動
する
全空間を一様に分割するのは効率が悪い
何らかの方法で原子付近だけ解像度を上
げたい
局所分割
座標変換
...
メッシュの取り方 (II):局所分割
原子が動いた場合
の扱いが難しい
並列計算機上での
振舞いはあまり芳し
くない
解像度を極端に上
げる時には有力
メッシュの取り方 (III):曲線座標系
原子の近くでメッシュ
が密になるよう座標
変換
原子が複数あるとき
は線形に足合わせる
原子が移動すると
座標系も追随する
並列化し易い
メッシュの取り方 (IV): 座標変換
座標変換の式：
x → x は直接計算できないので、NewtonRaphson 法を使って反復して求める
極端に adaptation が強いと変換不能になる
２倍強程度ならば特に問題はない
電子状態計算の方法 (I)
Kohn-Sham 方程式を有限要素法で展開：
H, F は巨大な疎行列
バンド幅は２００程度
有限要素法の場合は基底関数が非直交
なため、一般化固有値問題になる
電子状態計算の方法 (II)
ハミルトニアンの次元をNとすると、N/100
～N/1000 程度の固有値を（下から）求め
る必要がある
N は１００万を超えることもある
求めたい固有値の上限と、その次の固有値の
間には通常ギャップがある
• 「半導体」「絶縁体」に相当
• 「金属」ではギャップがない→ここでは考えない
電子状態計算の方法 (III)
直接ハミルトニアンを変形するのはほぼ
不可能なため、反復法で計算する
固有値問題を最適化問題に変換して解く
求めたい固有値が一つの場合：
レイリー商
をψについて最適化すると、最小の固有値
＋固有ベクトルが求まる
電子状態計算の方法 (IV)
複数の固有値を同時に求めたい場合：
をψについて最適化すれば良い。ただし
最適化問題に変換するメリット
「一般化」固有値問題であってもほとんど
手間が変わらない
特に求めたい固有値群と、その次の固有値の
間にギャップがある場合には非常にうまく行く
全ての固有ベクトルを同時に求めるので、
ブロック化が容易→ Flops が高い
H が非線形であることを考慮し易い
第一原理計算特有の事情
最適化の方法
準ニュートン法の省メモリ版（Limited memory
BFGS method）を使用
Liu and Nocedal, Math. Prog. 45, 503 (1989).
国内ではマイナーだが、６００回以上引用
通常の準ニュートン法の収束性を保ちつつ、メモリ
の使用量を大幅に削れる
非線形問題に関しては共役勾配法よりも高速
• １反復毎に目的関数を一回、微分を一回、前処理一回
計算すれば良い（一次元の最適化がほぼ不要）
準ニュートン法の収束
波動関数の初期値
最初のＭＤステップ
粗いメッシュで乱数から解いた結果を利用
（1-way multigrid）
２ステップ以降
直近２ステップ程度の波動関数から「外挿」
簡単な割に効果が大きい
ポアソン方程式の高速解法 (I)
電子状態が更新される度に、電子密度が
作るポテンシャルを再計算する必要がある
周期境界条件の下でポアソン方程式を解く
ことで計算できる
ポアソン方程式の高速解法 (II)
曲線座標系を使っているため、ＦＦＴでは済
まない
大規模な線形方程式に相当する
Ax=b の形
A は疎行列 / x, b は通常のベクトル
直接解法は現実的ではない
共役勾配法のような反復法を使う必要がある
ポアソン方程式の高速解法
(III)
共役勾配法は系が大きくなると収束率が
急激に悪化する
多重格子法（マルチグリッド）を用いた解法
系のサイズによらず、一定の収束率を持つ
単独で使うよりも、共役勾配法の前処理として
利用した方が効率が良い
マルチグリッド法 (I)
標準的な「Ｖサイクル」を使用
S=Smoother, E=Exact Solver
色々な波長の残差を同じ比率で減らせる
Ｓ
Ｓ
Ｓ
Ｓ
Ｅ
マルチグリッド法 (II)
Smoother
収束率が良い方法は大体逐次性が強い
• ガウス・ザイデルなど
並列性やメモリアクセスの観点からチェビシェフ
多項式を使用している
ブロック対角型前処理
Exact solver
教科書的には厳密に解くべきだが、実際には
Smoother を何回か施せば十分
レイテンシーを抑えられる
ポアソン方程式の収束
Part III:
最近の展開
Ｔ２Ｋ上でのパフォーマンス
東大センターのT2Kシステム上でベンチ
マークを取った
理論ピーク＝9.2 GFlops/core, Quad
ASC-P32 は 4.0 GFlops/core, Dual
理想的には2.3倍速くなるべき
ベンチマークの詳細は東大情報基盤セン
ターニュース（2009年2月）
単体性能（小規模系）
ダイヤモンド結晶のＭＤ（６４原子）
P32:112.9秒（PGI+ACML）
T2K:112.7秒（日立+netlib）
T2K:87.8秒（IFC+MKL）
T2K:87.1秒（IFC+ACML）
インテルコンパイラ＋ACMLが最有力
P32 から 1.3 倍程度の高速化
• 最新バージョンだともう少し出る
BLAS/LAPACK 部分は理論性能近く出ている
• 系が大きければもう少し改善するはず
並列化の方法
実空間を「曲線座標上で」分割し、各領
域をプロセッサ毎に割り当てている
各プロセッサは同数の要素を持つ
負荷分散はほぼ完全
必要な通信は主に隣接プロセッサ間での、
波動関数等の通信
グローバルな通信は最小限で済む
コア単位でＭＰＩ並列
曲線座標系
並列性能（中規模系）
水１２５分子系のＦＰＭＤ計算
NUMA は指定した方が良い（特に高並列時）
FEMTECK の profile 分析（Ｉ）
100%
80%
その他
通信
O(N)
O(N^2)
O(N^3)
60%
40%
20%
0%
水１２５分子
水６４０分子
FEMTECK の profile 分析（ＩＩ）
O(N) の部分はスカラー演算が多い
行列を生成する部分に相当
規則性が低いので、ベクトル・GPU には不向き
O(N2), O(N3) の部分はほぼ行列演算
極力 BLAS3 を使うようにしている
アクセラレータの恩恵を受けるはず
特に大規模計算では大きな効果が期待できる
混合精度計算の検討 (I)
単精度／倍精度の性能の差
Intel, AMD: ２倍
GPU, CELL: ～１０倍
次世代スパコン：？？
性能差が大きい程、工夫の余地が広がる
メモリ使用量は半分にできる
通信コストも（ほぼ）半分になる
混合精度計算の検討 (II)
第一原理計算は「反復法」
最終的な波動関数の精度は５～６桁くらい
全てを単精度で計算するのは困難
• 特に収束領域で飽和しやすい
ある程度収束するまで単精度で計算するのは
有力
収束領域でうまく行く方法を検討
グラム・シュミットの直交化 (I)
波動関数の正規直交化
一種のQR分解
「DTRMM」 に相当
グラム・シュミットの直交化 (II)
そのまま単精度化すると大幅に精度が落
ちるので、
と分割して（L は三角行列）
と変換すれば良い
反復部分の効率化 (I)
残差の計算：
本当の式はもっと複雑
桁落ちし易い
収束領域で単精度化するのは困難
反復部分の効率化 (II)
未知変数の更新：
「前処理」部分は収束領域でも単精度化できる
前処理部分に手間をかけて反復回数を減らし、
「単精度化しにくい」残差計算を最小限にしたい
単純な対角型前処理は損と言える
混合精度化による高速化
単精度が２倍速いシステムであれば、
FEMTECK の場合には全体で 20-30%
程度の高速化が期待できる
１０倍の差があれば、２～３倍？
単精度と倍精度の性能比は最終的にど
の位に収まる？
まとめ
有限要素法に基く第一原理電子状態計算に
ついて紹介した
大規模・高精度向けの手法
条件が整えば数千台規模の分散メモリ型並列
計算機上でも十分な性能を発揮できる
次世代スパコンに向けてチューニングを進める
予定
濃リン酸の性質
リン以外全ての原子が水素結合に寄与しうる
プロトン移動が頻繁に起こる
→ 液体の構造は複雑
実験を中心に数多くの研究が行われてきた：
古典ＭＤ：定性的には良いが、定量的にはかなり
誤差がある
第一原理計算から信頼できる情報を出したい
特に動的な性質やプロトン移動を記述したい
リン酸のまとめ
動径分布関数など、実験と直接比較可能な
物理量は定量的に一致している
リン酸分子の約30%はイオン化している
各プロトンは平均してほぼ１本の共有結合と
１本弱の水素結合を持つ
プロトン移動が頻繁に起きている事を確認
この際、共有結合と水素結合が切り替わる
同時に P-O と P=O も入れ替わる
基底関数の基準
波動関数は通常 singularity を持たない
ある程度高次の基底関数を使った方が、
基底の数を減らせる
ただし、基底関数同士の重なりが増えるので、
ハミルトニアンが密行列に近づく
三次元問題では、三次か五次程度の基底
関数が最適と考えられる