＊ここに掲載するのは、ユニバーシティカレッジ・ロンドンのLianheng Tong博士によるHECToRレポート「dCSE Project on Improving Automatic Load Balancing and Molecular Dynamics in Conquest, Lianheng Tong, University College London, 2013/04/17」を要約したものです。
概要
物質の構造と特製の研究において、量子力学的手法の利用は現代科学の基礎となっています。CASTEPやVASPといった一般的に用いられるDFTコード[5]は、HECToRにおいて利用頻度が高く、その計算負荷は原子数の3乗に比例し、メモリー容量は2乗に比例します。この性質により計算可能な系のサイズは多くて数千原子に制限されます。生物科学においては、この問題を回避するのに、古典的な力場に囲まれた小規模の量子力学系を用いるのが一般的です。
この15年来、線形にスケールするDFTコードの開発が研究されています。ConquestはHECToRで稼働する線形スケールコードの一つで、1万原子を超えて数百コアを効率的に用いる事が実証されています[2]。
このプロジェクトの目的は、Conquestの負荷バランスを改善し、10万原子へMDシミュレーション規模を拡大することです。対処する問題は、プロセッサー間に計算を最初に分散することと、シミュレーション中に動的に負荷バランスを行うことです。
初期負荷バランスには現在3種の方法が用いられています：空間充填曲線[3]を用いた自動手法、ユーザパラメータに従ってシミュレーションセルを分割してデータをプロセッサへ割り当てる外部Pythonライブラリによる方法（これは自動手法が不十分な場合よく用いられます）、および（あまりよく利用されない）焼き鈍し法です。自動手法はヒルベルト曲線[7]を用いますが、これは空間的にコンパクトな利点があります。つまり、原子やグリッドデータが空間的に近く、同じプロセッサに割り当てられる可能性が高く、通信量が少なくなります。標準的なヒルベルト曲線と空間充填曲線は、一般に負荷バランサに固定の制約を課します。これは、全てのシミュレーションセル分割が3方向で同じで、2のべき乗個でなくてはならないことです。この制約は、セルが立方体とみなされない場合は適しません。このプロジェクトが行う自動手法の改善は、曲線のコンパクト性を保ちつつヒルベルト曲線上の制約の緩和することです。
Conquestでは、行列保存形式とインデックスはシミュレーションセル内の原子の位置に依存しており、分子動力学の実装をより困難なものにしています。現状のConquestはボルンオッペンハイマーMDアルゴリズムを実装しており、シミュレーション中に同じパーティションに原子は配置されていることが仮定されています。これは小さな原子の移動ならばよいですが、原子が遠く離れた位置へ移動すると、隣接リストが不正になり、動作に適しません。さらに、原子の移動は隣接リストの変更を生じますが、行列要素のプロセッサ分担は原子の位置に依存しているため、スパース行列形式とプロセッサ上のサイズはシミュレーションの進行に従って変化します。Conquestはこの問題を、行列をリセットして各時間ステップで最初から計算を開始することで回避しています。これは常にボルンオッペンハイマー的に正しい結果を生成しますが、分子動力学では不十分です。このプロジェクトは、各MDステップでインデックスを動的に割り当てて、行列要素を再配列する実装を行います。こうすることで、どのような原子の移動にも対応した安定的なMDシミュレーションと、次のステップで計算済みの行列要素を再利用することを可能にします。
自動領域分割
Conquetはシミュレーションセルを均等なパーティションに分割します。これは原子データのプロセッサへの割り付けの基本ブロックです。ユーザの指定がない場合は、実空間グリッドポイントはパーティションに従ってプロセッサへ分配されます。つまり、特定のパーティションを担当するプロセッサは、そのパーティション内のグリッドポイントブロック群にも対応します。よって良好な領域分割が負荷バランスに必須です。
Conquestは、一つのパーティションが許容する原子の最大数をベースにしてパーティションセルサイズを推定します。この最大数はユーザ定義が可能です。コードは、均一な原子分布を仮定した場合の初期値から初めて、各方向にパーティション数を徐々に増やし、パーティション内の最大原子数が指定の数を超えない状態になるまで続けます。パーティションが得られたら、標準的なヒルベルト曲線法を用いて、パーティション座標(ix,iy,iz)をヒルベルト曲線上の一つの1次元整数へ変換します。この整数インデックスを用いて、パーティションをプロセッサへ割り当てます。この時、各パーティションの原子数で重みをつけられるため、プロセッサはほぼ同じような原子数を担当することになります。
ヒルベルト曲線を用いてパーティションを1次元化する事には、ヒルベルト曲線が空間的にコンパクトであるという明確な利点があります。つまり、分割数にかかわらず、空間的に近いパーティションはヒルベルト曲線上で近くなることを意味します。一旦パーティションがプロセッサへ割り当てられれば、そのプロセッサに属する原子群は空間的に近くに存在し、ヒルベルト曲線上で離れており、つまり異なるプロセッサに割り当てられているパーティションに属する原子群は、空間的に離れている可能性が高くなります。こうしてMPI通信量は削減され、有限距離内の相互作用のみが考慮されることになります。
しかしながら標準的な3Dヒルベルト曲線法にも不利な点があります。ヒルベルト曲線の再帰的な生成手法では、3次元方向の分割は同数、かつ2のべき乗でなければなりません。シミュレーションセルが立方体でない場合は、最短次元方向に不必要に多くの分割が生じます。最良の場合でも、空のパーティションが多数作成されて計算が遅くなります。 最悪の場合、これは重大な負荷の不均衡を招き、最終的にコード実行はす破たんします。この事例がこのレポートに示されています。
ここでは、異なる空間次元で異なる再帰レベルを可能にする3Dヒルベルト曲線法：非立方体ヒルベルト曲線法を構築しました。
·負荷バランス
既存のConquestは、パーティションのプロセッサへの割り当てにおいては、負荷バランスのために、パーティション数か原子数で重みをつけて行います。通常は、原子がシミュレーションセル内に均一に分散されていたとしても、パーティション数の重みは用いません。このような場合でも、各パーティション内の原子数は様々なので負荷分散が生じます。空のパーティションを担当する場合も可能性としてあり、この場合は実行に失敗します。
生体分子のような多くの原子種がある場合、原子数による重みは負荷分散を招きます。これはサポート関数数が原子種により大きく異なるためです。
既存の領域分割コードはConquest内に深く埋め込まれ、修正するには労力がかかりすぎるため、新しい領域分割と負荷バランスコードを作成して組み込みました。さらに下記の新規機能も同時に作成しました：
・系の形状を考慮した自動領域分割の組込み
・サポート関数による重み付けの組込み
・各方向の分割数の指定機能
テスト結果
新しい非立方体ヒルベルト曲線領域分割のテストのために、以下のデータを用意して既存の分割と比較しました：
- BulkSi512C：512原子バルクシリコン(41.04 a.u.×41.04 a.u.×41.04 a.u.)
- BulkSi512F：512原子バルクシリコン(82.09 a.u.×82.09 a.u.×10.26 a.u.)
- BulkSi512L：512原子バルクシリコン(656.72 a.u.×82.09 a.u.×10.26 a.u.)
- SlabSi2048M ：セル中心に配置された2048原子のxy平面スラブシリコン (82.09 a.u. 82.09 a.u. 92.35 a.u.)
- SlabSi2048W ：セル下端に配置された2048原子のxy平面スラブシリコン、原子の半分はセル上端へ周期境界条件が適用される(82.09 a.u. 82.09 a.u. 92.35 a.u.)
- GeSi5333：シリコンスラブ基盤上のゲルマニウムhut cluster、5333原子、スラブはセル下端のxy平面に配置 (164.18 a.u. 205.22 a.u. 51.31 a.u.)
- GeSi22746：シリコンスラブ基盤上のゲルマニウムhut cluster、22746原子、スラブはセル上端のxy平面に配置されセル下端へ周期境界条件が適用される(328.36 a.u. 328.36 a.u. 89.79 a.u.)
- GeSi39130：シリコンスラブ基盤上のゲルマニウムhut cluster、39130原子、スラブはセル上端のxy平面に配置されセル下端へ周期境界条件が適用される(328.36 a.u. 328.36 a.u. 97.48 a.u.)
- DNAH2O3439：水中のDNA、3439原子、サポート関数重みによる負荷バランスが、原子群重みよりどのくらい 性能向上するかを検証するのに用いる
テストは全て、補助行列トレランス10-4Haとした、非セルフコンシステントエネルギー最小化1点計算をHECToRフェーズ3で行いました。コンパイラとライブラリは、Cray compiler suit PrgEnv-cray/4.0.46とxt-libsci/12.0.00です。コンパイルフラグは-O3を用いました。
·セル形状別の性能
512原子バルクシリコンに対して、自動領域分割を行いました。原子種が1種類しかないのでサポート関数重みの効果はありません。
既存アルゴリズムは、BulkSi512Lでは並進方向の分割数が多すぎて実行が中断されました。また、BulkSi512Cから平面シェルのBulkSi512Fへ大きな性能劣化がみられました。これはからのパーティションの生成によるものです。新しい領域分割は全てのケースで最適な負荷バランスを示し、丸め誤差に対する補正の点からも正しい原子配置を示しました。
より大きな系に対しては、さらに大きな性能改善を見せました。既存コードは、SlabSi2048Mに対してメモリーオーバーで実行中断されました。これは原子数が多いのではなくか空のパーティション生成による結果です。あまりに多くの積分ブロックポイントが割り当てられたことが原因です。
·Ge Hut Clusterケースの性能
シリコン基板上のGe Hut Clusterについては、既存コードは全てのケースでメモリーオーバーで実行できませんでした。既存コードは全ての原子を全てのプロセッサに割り当てることが出来ずに失敗しています。これに対して新しいコードは領域分割と負荷バランスの改善により実行可能です。
しかしながら、既存コードでもユーザ定義の外部ライブラリによるパーティションファイルを用いれば実行可能です。
結果として新しい領域分割のほうが若干高速でした。これはプロセッサ間の通信量が減ったことによるものです。
·水中のDNAケースの性能
これは生体系に特有な、原子に関する負荷バランスが最適な性能を生じないケースです。サポート関数の重み付けにより既存コードに比べて約1割の性能向上がみられました。
動的再割り付けによる分子動力学
既存のConquestでは、最初にパーティションへ割り当てられた原子は、シミュレーションが終わるまでそのパーティションに固定されます。各原子はシミュレーションセル内のグローバルインデックスを持ちますが、実際にプロセッサが度の原子を所有するかや、離れた場所にある隣接原子に関連するデータをどこでフェッチすべきかは、パーティションとハロに基づくインデックスシステムに基づいて行われます。これはプライマリセット内の原子の指定された相互作用範囲内の少なくとも1つの原子を含むパーティションの集合です。
MPIプロセスのプライマリーセットはプロセスに所有された全ての原子のセットです。もしMD中に、原子がそのパーティションから別のパーティションへ移動、あるいはシミュレーションセル境界を跨いでも、設定は更新されません。よって、パーティション移動が発生しても感知されず、元のパーティションはその原子が依然として所有しているように扱います。これは隣接リストでのエラーとグリッド上の値の配列送受信でのサイズ矛盾を生じます。この回避のために、通常は時間ステップ数回に一度、MD計算をリスタートすることが必要です。これがConquestの大規模MDシミュレーションを困難にしています。
さらに、Conquestは各MDステップの最初に、補助行列に対するMcWeenyの初期化を行います。補助行列はエネルギー最小化処理内で用いられる変数です。もし原子がそれほど遠方に移動しなければ、次のMDステップで補助行列を再利用して、McWeeny初期化もスキップすることが可能です。
よってConquestのMDの改善点は2つあります。最初に、パーティションとそのデータを動的に更新する機能を追加します。これにより原子の移動の度にリスタートする必要がなくなります。次に、原子の移動に際して旧所有者から新所有者へ行列データを移動して、補助関数を前のMDステップから計算しなおす方法を実装しました。
この作業は、日本のNIMSのMichiaki Arita博士、Tsuyoshi Miyazaki博士と協力して行いました。
·テスト結果
テストには、様々な大きさの氷の立方体を用いました。MDの設定は以下の通りです：
・初期イオン温度：300 K
・時間ステップサイズ：0.5 fs
・DFT汎関数：PW92(LDA)
・セルフコンシステンシー：エネルギー最小化と電荷セルフコンシステンシーの両方
768原子の水のMDシミュレーションを、HECToRフェーズ3の1ノード32コアで実行しました。ステップ毎にMcWeeny初期化を行いました。このやり方と既存コードの違いはメンバーデータが更新される点です。よって以前のように原子がパーティション境界を跨いでも実行は中断されません。このとき、エネルギー収束条件として補助行列トレランスを10-3としました。
McWeeny初期化に時間が掛かるため、比較は12時間かかって136MDステップのみで行いました。32コアという数は、768原子系にとって既に大きなもので、負荷バランスが悪くなります。既存コードのままでは高速化のためにコア数を増やすことは現実的でありません。
結果は、シミュレーション時間に対してエネルギーの振動が共に示されました。MacWeeny初期化版の平均エネルギーは一定値ですが、補助関数再利用を行うコードは振動は素早く減少してより低いエネルギーへ収束していきます。
また、補助行列再利用手法はMcWeenyよりも40％高速でした。これは主に、McWeeny初期化コストが削減されたためです。
拡張ラグランジアン・ボルンオッペンハイマーMD
補助行列再利用は低いトレランスでは劇的に高速化しましたが、より厳しいトレランスの場合は効果は減少してしまいます。
この問題の解決策の一つに、現状のMDを拡張ラグランジアン形式[6]へ拡張することです。これは、ボルンオッペンハイマーMDへ電子的な自由度を組み込むことです。同様に補助行列再利用も可能です。
拡張ラグランジュ形式を用いたConquest を8原子のシリコンセルで実行しました。結果は、トレランス10-3を用いた場合、エネルギーの振動は生じませんでした。
8原子シリコンセルの結果は、拡張ラグランジアン形式は補助行列再利用版よりも時間が掛かり、繰り返し数も増加しますが、安定性は優れていました。
謝辞
このプロジェクトは、nAG Ltd.が運営するHECToRの分散計算科学および工学(CSE)サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびnAG Ltd.のCSEサポートサービスにより管理運営されています。
文献
| [1] | D. R. Bowler, T. Miyazaki, and M. J. Gillan. Parallel sparse matrix multiplication for linear scaling electronic structure calculations. Comput. Phys. Commun., 137(2):255-273, June 2001. |
| [2] | D. R. Bowler, T. Miyazaki, and M. J. Gillan. Recent progress in linear scaling ab initio electronic structure techniques. J. Phys.: Condens. Matter, 14(11):2781, 2002. |
| [3] | V. Brázdová and D. R. Bowler. Automatic data distribution and load balancing with space-filling curves: implementation in CONQUEST. J. Phys.: Condens. Matter, 20(27):275223, 2008. |
| [4] | F. Gray. Pulse code communication, 1953. |
| [5] | W. Kohn and L. J. Sham. Self-consistent equations including exchange and correlation effects. Phys. Rev., 140(4A):A1133, 1965. |
| [6] | A. M. N. Niklasson. Extended born-oppenheimer molecular dynamics. Phys. Rev. Lett., 100(12):123004-, Mar. 2008. |
| [7] | H. Sagan. Space Filling Curves. Springer, 1994. |
| [8] | J. Skilling. In G. Erickson and Y. Zhai, editors, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 23rd Int. Workshop, number CP707, page 381-7, New York, 2004. American Institute of Physics. |