プログラムの高速化・並列化サービスの事例

＊ここに掲載するのは、ケンブリッジ大学のMike Towler博士らによるHECToRレポート「Adding the molecular dynamics functionality to the quantum Monte Carlo code CASINO, Mike Towler1), Norbert Nemec2), 1)Cambridge University, 2)UCL, November 15, 2011」を要約したものです。

概要

拡散モンテカルロ（DMC）ウォーカー数を僅かに変化させた試行関数へ更新するGrossman and Mitas [1]のアイデアに従い、分子動力学（MD）シミュレーションに対してDMCを実行することを行います。現段階$2010$の形式においてはCASINOは、保存された前回のウォーカー数を用いてDMCを継続実行します。本来はシステムの変更無しに継続計算する様設計されましたが、試行波動関数ファイルを単純に置き換えることで問題なく動作し、平衡までの時間を削減できました。

ウォーカー群の更新の最も単純な実装のやり方は、DMC継続実行時の最初にディスクからウォーカー数を読み込むようにルーチンを修正することです。ユーザ観点ではこの処理は以下の様になります：

1. ・MD時間ステップで新たなCASINOプロセスを起動する。
2. ・次のMD時間ステップはCASINOを引き継いで実行される。前回のMDステップが書出したconfig.outファイルは、config.inと改名されて次のMDステップの入力として読み込まれる。
3. ・新しいフラグdmc_reweight_configsによりconfig.in入力後に再重み付けステップが実行される。
4. ・新しいフラグdmc_spacewarp_configsによりconfig.in入力後にスペースワーピングス$Space warping$テップが実行される。

実装は以下のステージで行います：

1. I. 再重み付け[2]：各ウォーカーの統計的重みは、新旧の波動関数の二乗の絶対値間の比により調節されます。必要なデータのほとんどはconfig.outに保存されており、実装はtwist average計算部分のみです。
2. II.  Space warping[3]：継続計算では新旧両方のイオン座標が必要なため、イオン座標はconfig.inへも保存する必要があります。こうすれば、space warpingを再重み付けループに実装できます。

ケーススタディ：シリコン結晶

ここでは、Car and Parinelloにより計算されたMDシミュレーション[4]に従います。対象とした系の条件は以下のものです：

・ダイヤモンド構造のシリコン結晶
・ユニットセルは8原子を含む
・Trail-Needs擬ポテンシャル
・LDA
・カットオフエネルギー：200Ry
・格子定数：10.26Bohr
・光学フォノンモード（Γ点フォノン）

初期調査から以下のことが判りました：

1. ・DMC継続計算は平衡化計算の除去には極めて効果的である。
2. ・この系では再重み付けは全く必要性が無い。DMC平衡化計算は正確な配位に鈍感である。一つの波動関数で平衡化されたウォーカー数は、似た波動関数に対してもほぼ平衡化されている。これは一般に一様な結晶に対しても同様と考えられる。
3. ・space warpingの利得はただ一つ、再重み付けの性能改善だが、この系に対してはspace warpingを用いても改善は見られなかった。これは、電子が特定の核に束縛されず全結晶に広がっているためである。
4. ・この時系列MD計算では、DFT、VMC、DMCのエネルギー差は定数に支配されている。DMCはより低コストの手法の結果と比較して、如何なる異なる物理的な特徴を示さない。この場合は再重み付けとspace warping は効果的だが、DMCを利用した場合に系に対してこれがどの様に働くかは不明である。

シリコン結晶では、DFT、VMC、DMCのエネルギー差が一定である事しか特徴がありません。このような場合はDMCはVMC以上の利点を発揮できません。
本当に知りたいのは、VMCとDMCのエネルギー差が一定でない場合に何が起こるのか、です。DMCの真の目的は、エネルギー絶対値の改善ではなく、新しい物理の探索あるいは物理的に関連するエネルギーの差をより正確に求めることです。DMCの計算コストは、DMCによるエネルギー曲面がVMCに比べて一定ではない差を見せる場合にのみ正当化されます。

こうして、DMCの虚数時間発展形式の開発が必要なことが明らかになりました。試行波動関数の比を前提にした再重み付け処理では、試行波動関数には無い特性を正確に出現させることを期待できません。
よって、DMCエネルギーがVMCとは質的に異なる特性を示す系の研究に大きな期待が寄せられます。ファンデルワールス錯体はそうした系で、その研究に必要な相関がVMCには無くDMCに存在します。このためにはDMCが必要です。

CASINOへのDMC-MD機能の組込み

パブリックバージョンのCASINOやDFTコードへのDMC-MDの組込みにおいては、$1$単純なコマンド実行可能なユーザフレンドリーな形式が望ましく、$2$CASINOがサポートする全てのアーキテクチャへの対応（シングルコアPCからペタスケール・スーパーコンピュータのバッチ処理まで）が必要です。

最終的にDMC-MDはCASINOへ組み込まれ、最初に2011年11月にバージョン2.10として公開されました。また、quantum espressoプロジェクト（DFTコードpwscfが含まれます）のPaolo Giannozzi博士らと相談して、CASINOのDMC-MD計算機能がバージョン4.3に組み込まれました。
PwscfへのCASINOサポート追加、擬ポテンシャルとコンバーター追加、runpwscfおよびrunqmcmdスクリプトの変更、追加を行いました。

並列効率ベンチマーク

一般的には量子モンテカルロ法は本質的に並列処理であり、そのため、新世代および次世代の大規模並列コンピュータに理想的には適しています。これは言うまでも無く、VMCとその波動関数最適化に用いられる様々な技法に当てはまります。純粋なVMCの場合は、シミュレーション中に一切プロセッサー間で通信は発生しません。各プロセッサーは、個別の波動関数と乱数系列を用いた独立したランダムウォークを実行します。得られた結果は最後に全プロセッサーに渡って平均化されます。平衡化の計算時間を無視すれば、N_p個のプロセッサーで経過時間Tで計算される場合、同じデータを用いるとシングルプロセッサーではT N_pの時間が掛かります（無論、ランダムウォークは2つのケースで異なるため、結果は統計誤差内で一致します）。よってVMCは任意のプロセッサー数に対してスケールすべきです（このため今後は波動関数最適化手法には言及しません）。

よって問題があるとすればDMCアルゴリズムに存在することになり、その問題の大半は負荷バランスです。DMCはVMCと似たやり方、つまりウォーカーを異なるプロセッサーへ割当てることで並列化されますが、対照的にDMCでは通信が発生ます。負荷バランス処理を行うということは、ウォーカー数を動的に変更することを意味します。つまり、波動関数の形状を変えてウォーカーは中断されたり複製されたりすることで、実行中にその数が揺らぎます。ここでプロセッサー間通信が発生する理由は、シミュレーション中に、ある種のフィードバック機構を通じて、この数が初期ターゲットへ適合するように調節されるためです。これは各時点の全エネルギーに基づいているため、各時間ステップ後に全プロセッサーに渡って計算と平均化が必要になります。

しかしながら最大の問題点は、コア間でウォーカーを転送する必要がある事です（ここではウォーカーは、配位に対する現状の電子位置リストであり、エネルギーと波動関数に関係する様々な値を伴います）。この転送は純粋に効率のための操作です。時間ステップのコストは最大集団を持つプロセッサーにより決まるため、負荷バランスのためには各コアに同程度数のウォーカーを持たせることが重要です。プロセッサー間で通信するウォーカーの総数はコア数と共に増加し、大規模並列マシンでの実行において唯一かつ最大の非効率性を生じます。

ここでは、強スケール性をベンチマークします。波動関数のサンプル回数をNに固定してコア数を変化させた場合の実行時間を計測します。Nはウォーカー数×イオン移動数です。
最初に、ウォーカー数もイオン移動数も固定した場合を考えます。理想的な場合は、プロセッサー数を半分にすれば時間も半分になります。これはQMCを並列マシンで利用する最良の方法とは言えませんが、幾つかの重要な特性を見ることが出来ます。

50炭素原子を含む2D周期セルで表現されたグラファイトシート状に一つの水分子が吸着する、という典型的なDMCシミュレーションを行いました。ウォーカー数は486,000個に固定し、JaguarPFマシン（米国オークリッジ国立研究所の224,256コアのCray XT5）でベンチマークを行いました。このマシンの特性からコア数は12の倍数としました。このウォーカーに対して10回の統計的なイオン移動を行った場合、CASINO2.6は約2000コアまでは程々にスケールしますが、それ以降は劣化して5000コアを超えると遅くなります。

並列効率改善

この性能劣化問題に対して、ウォーカーの再配置コストを完全に消去しました。この性能改善は以下の様に行いました：

1. ・非同期・ノンブロッキングMPI通信の実装：既存実装で用いられているブロッキングMPI_SENDとMPI_RECVをノンブロッキングタイプに変更し、同時に計算をオーバーラップさせました。
2. ・プロセッサー間で転送されるウォーカーの決定処理の並列化：マスタープロセスは、各コアの現状のウォーカー数を受信後に、ウォーカーの転送を決めます。このアルゴリズムは、送信と受信側のプロセッサー要請を満足するように、かつ転送数を最小にするように、転送を決定します。残念ながらこの処理コストはプロセッサー数に比例して増加し、最も効率的な転送の決定は実際の転送処理よりも高コストになります。この処理は通常は目立ちませんが、大規模プロセッサー数ではCASINOの実行に大きな制約となります。そこで、約500コアの小さなグループで転送を決定するように修正したところ、前述の例では1，2秒で終了しました。

この結果82944コアでは以前に比べて3倍以上高速化しました。

大規模並列マシンではプロセッサー数を増やすとプロセッサー当たりのウォーカー数が減少してしまうため、スケール性が制限されます。よってウォーカー総数固定ではなくて、プロセッサー当たりのウォーカー数で考えなくてはなりません。
ここで採った方法は、N個のウォーカーから始める場合、このウォーカー数とプロセッサー数を2倍にして、信頼区間を維持するためにイオン移動数を半分にすることです。この方法では、コア当たりの実行時間は半分となり、余分な通信は発生しません。

結果として82944コアにおいて既存コードに対して3割高速化し、線形にスケールしました。またウォーカーの再配置コストは412秒から1秒へ短縮化されました。

謝辞

このプロジェクトは、nAG Ltd.が運営するHECToRの分散計算科学および工学$CSE$サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびnAG Ltd.のCSEサポートサービスにより管理運営されています。

文献

[1]	J. C. Grossman and L. Mitas, Phys. Rev. Lett. 94, 056403 $2005$.
[2]	N. Nemec, Phys. Rev. B 81, 035119 $2010$.
[3]	C. Filippi and C. J. Umrigar, Phys. Rev. B 61, R16291 $2000$.
[4]	R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 $1985$.