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

目次

1  イントロダクション

2  作業1：ハイブリッド並列DL_POLY_4の実装

2.1  イントロダクション

2.2  実装

     リンクセルへの原子の割当てとVerletリストの計算

     力の計算

       長距離力

       拘束力

       通信

2.3  結果

2.4  現状

2.5  結論

3  作業2：百万原子計算

3.1  イントロダクション

3.2  64bitコードテスト

4  結論

5  謝辞

1  イントロダクション

DL_POLY_4[1]は、STFCデアーズベリー研究所[2]のI.T.TodorovおよびW.Smithにより開発された、一般用途古典分子動力学(MD)シミュレーション・ソフトウェアパッケージです。このソフトウェアはモデルの原子レベルの時間発展モデリングのために用いられ、物性科学や、固体化学、生体シミュレーションおよびソフトマター研究において広く利用されています。また、CCP5[3]の後援の元に開発され、英国および世界的に広く利用されています。このソフトウェアは完全な分散データコードであり、空間領域分割、リンクセル、Verletリスト(Verlet neighbour　list)[4]、および3Dデアーズベリーフーリエ変換(DaFT)[5]といった手法を用いています。コードの性能は極めて高く、大規模なプロセッサー数に対してスケールします。

本プロジェクトの目的は、DL_POLY_4で効率よく利用可能なコア数の増大のみならず、研究に用いる系のサイズの拡大にあります。最初の目的については、OpenMPを用いた第2レベル並列化の導入により、リンクセルの制約を回避します。これはセクション2で詳しく述べます。系のサイズ問題はDL_POLY_4が全体を通して32ビット整数を用いていることが要因のため、本プロジェクトの第2の目的としては、DL_POLY_4を64ビット整数化し、HPCシステム上で現在研究可能なサイズを超えて計算対象系のサイズを最大限に拡大することです。これはセクション3で記述されています。

2   作業1：ハイブリッド並列DL_POLY_4の実装

2.1  イントロダクション

上で述べたようにリンクセル・アルゴリズム[6]はDL_POLY_4が並列化を行う基盤であり、均等空間領域分割を用いています。シミュレーションにおけるセルは、多数の同体積かつ同形の並行四面体に分割され、その各領域はMPIプロセスへ割当てられます。多くの原子間相互作用は有限距離(カットオフ)を超えると無視でき、各領域のエッジ近傍の原子座標データだけは原子間相互作用計算のために隣接コアと通信が必要であり、カットオフは領域エッジ長よりかなり小さな値とされます。リンクセル法は、グリッド上の微分方程式を解く際の袖領域交換(halo exchange)の手法と多くの点で似ていますが、その違いは空間内に不規則に配置され、時間ステップで移動する一般化されたグリッド点、即ち原子を用いる事です。

しかしながらMPIプロセス数が増えてくると、リンクセル・アルゴリズムは限界を示し始めます。多くのMPIプロセスは小さな領域で実行され、かつその領域の辺がカットオフより大きいものでなくてはなりません。実際に特に大きなカットオフを持つ原子間相互作用ポテンシャルによるシミュレーションでは、リンクセル・アルゴリズムの性能劣化は利用可能なコア数にさらに大きな制約になりますが、その時点までは良好なスケール性能を示します。ここで、もしこの性能劣化がP個のプロセスにおいて生じ、さらにT個のスレッドが用意できてハイブリッド並列が可能ならば、T*P個の総コア数を利用できることになります。よって十分な数のスレッドベースの第2レベル並列化を導入すれば、このような問題を研究する際にさらに多くのコアを利用できて、比較的短いカットオフを用いた長期間の現象を研究するのに役立つと考えられます。

2.2  実装

DL_POLY_4への並列レベルの追加の実装は、原理的には直截的なものです。リンクセル法は問題に対して領域分割を用いる事から、少なくとも0次においてはMPIフレームワークは関係なく、そこで主になすべきことは以下の2点です

1. 各原子をリンクセルへ割当ててVerletリストを計算する
2. 様々な力の計算を実行する

MPIフレームワーク問題に行く前に、このオーダーでの実装と、スレッドを用いて性能が改善するかどうかについて、まずは議論します。

     リンクセルへの原子の割当てとVerletリストの計算

リンクセル法において単位セルは、出来る限りカットオフに近い辺の長さを持つ平行六面体に分割されます。この平行六面体はリンクセルと呼ばれます。これは、各原子はそれが属するリンクセル内の原子あるいは隣のセル内原子とのみ相互作用することを意味します。複数プロセッサーでこのモデルを実行するには、同数のリンクセルを各プロセッサーへ割り当てます。通信量を最小化するためには、空間領域のセットを各プロセッサーへ割当てるのが最良な方法です。その性能評価は通常良くある並列PDEソルバーと似ており：力の計算コストは領域の体積、つまりそこに含まれる原子数に比例し、通信量は表面積に比例し、これら2つの量の比が並列スケール性能を決定します。

従って、リンクセル法は相互作用を計算する際の、原子を検索すべき空間を劇的に減らすことにより、原子間力計算を加速します。これに対してDL_POLY_4はさらにVerletリスト(Verlet Neighbour List[7])を計算します。これは各原子に対する、そのカットオフ内に存在してそれらと相互作用する他の原子のリストで、力の計算を加速するものです。このようなリストがリンクセル法をさらに劇的に加速することは驚くべきことのように思えますが、立方体に内接する球の体積を考えた場合その体積比が約1.9であることから、この方法によって力の計算毎にほぼ2倍の効果があると期待できます。

リンクセルとVerletリストの活用は極めて効果的で、原理的には力の計算を$O(N^2)$から$O(N)$へ削減しますが、都合の悪い現象も生じます：力の計算は極めて高速になりますが、原子のリンクセルへの割当てと近傍原子リスト生成コストが実行時間において無視できない量となります。ある場合にはそれは実行時間の30%まで達します。こうして本プロジェクトは、

1. もし力の計算をOpenMPで加速することを考える場合には、この部分も含めて解決しなくてはなりません。
2. この部分へアルゴリズムの改良を組み込む第2のdCSEプロジェクトとなりました。

近傍原子リスト生成の並列化は直截的で、各近傍原子リストは他のリストと独立なため、各スレッドへ互いに異なる原子セットの近傍原子リストを生成させることとします。しかしながらリンクセルへの原子の割当ては少し込み入っています。シリアル処理では以下のような流れです

　For 全原子；
　　　原子がどのセルに有るかを評価
　　　該当セルのリストの最後に原子を追加

この処理の最後では、もし共有メモリー並列を用いる場合は注意が必要です。何故なら複数のスレッドが一つのリストに書き込む際に競合が生じる可能性があるからです。これに関して様々な方法を試しましたが、最終的に最も有効な方法は、各スレッドに独立なリストを作成させて、同期ポイントを設定するよりもシリアル実行させてリストをマージすることでした。これはかなり解きやすい問題ですが、プロジェクトを通して何度も生じた一つの問題を提示しています：リンクセル・アルゴリズムとその後の領域分割が生じる構造は、グリッド上でPDEを解くことを思い出させるものですが、原子が移動して密度が微妙に揺らぐという事実は、静的な通常のグリッドではないことを示しており、このコードは、与えられた時刻に隣にいる原子がどれか、言い換えれば与えられた領域に何個原子があるかを知ることが出来ません。これの意味することは、コード中の至る所に同じようなリスト生成を入れなければならず、これはシリアル実行と比較してスレッド化したコードにオーバーヘッドを引き入れることを意味します。

     力の計算

DL_POLYに実装された力場の一般系は以下の通りです

系のポテンシャルエネルギーは系を構成する全ての原子種の座標の関数として与えられています。これを見てわかる通り、力の計算の並列化は単純に展開された項数分あります。

多くの項は独立した変数、多くは常にそうとは限りませんが原子数のループで、それらは簡単に並列化可能ですが、注意深く述べるとすれば、エネルギー、力、テンソル力およびその他の量はループの最後にスレッド間で縮約(reduction)が必要です。これらの項は全く直截的ですが、コード内に静的変数を使用しており、OpenMP 2は並列領域でのエラーハンドリングの仕組みがないため、並列領域で必要となる変数の絞り込みには十分な注意が必要です。

しかしながらルーチンの数は絞り込むのでなく、別々に扱うことが望ましいです。特にこれらは下記の項を含みます

1. Ewald法における長距離成分
2. 原子の拘束による力

ここでは、各時間ステップの最後においてハロの再構築を行うルーチンを検討することも有用でしょう。

       長距離力

DL_POLY_4は原子間のクーロン力を計算するEwald法用にSPME[8]を用いています。この相互作用の1成分は長距離力であり、上述のリンクセル法では収まらず、その代わりにフーリエ変換の基づく方法で計算されます。この方法は以下のステップを踏みます

1. 実空間グリッド上の原子の電荷密度（に関連した量）を計算する
2. グリッドから波数空間へフーリエ変換する
3. 変換結果をEwaldポテンシャルの波数空間表示とコンボリューションする
4. 結果を実空間へ逆フーリエ変換する
5. エネルギー、テンソル力および原子に掛かる力を上の結果を使って計算する

3番目と5番目のステップはそれぞれ、グリッドおよび原子という独立した量に関するループであり、先述のやり方に適合します。1番目のステップに関して様々な手法を調査した結果、DL_POLY_4でFFT計算に用いられるDaFT[9]はOpenMPによる並列化が必要となりました。

グリッド上の電荷密度の計算は原理的には単純で、特定の領域の電荷密度へ寄与する原子のループであり、各原子は適切なグリッドへその寄与を付加します。よってOpenMPによる並列化を行うには、異なる原子のサブセットに対して電荷密度を計算するスレッドを生成し、異なる原子が同じグリッドへ寄与する際には競合を避けて適切に同期をとるようにします。この同期方法として3つの方法を検討しました：ループ後の各グリッド上でのreduction、および加算時のatomicとcriticalの使用です。reductionが最も性能が良いものと期待されましたが、グリッドはどの手法が有効かを調べるのにはメモリー的にコスト高であり、reductionは各スレッドにグリッドのコピーを持たせることが必要でした。実際にはatomicとcriticalはオーバーヘッドがあり全く有効ではなく、メモリーの問題はありますが、最終的にはreductionを選択しました。

DaFTに実装されたFFTは以下の2つのステップを踏みます

1. ローカルデータに対して1D FFTをシリアル実行する
2. 2つのMPIプロセスが持つデータを結合する。この時点で通信が発生し、複雑な回転因子が必要とされます。

最初のステップは、1D FFTルーチン内のClive Tempertonにより実装されて一般に利用されているFFTルーチンGPFA[10]を、ACMLルーチン[11]へ変更して並列化されました。これらが選ばれた理由には、それがスレッド対応されているだけでなく、幅広いHPCアーキテクチャーで高い性能を示すためであり、かつDL_POLY4での利用に関してライセンスの問題が発生しないためです。

次のステップには2つの手法を検討しました。それらは、スレッド並列がx,y,z方向(3D変換ステージに依存します)のベクトルに対して行われて同期通信が一回発生しすることと密接に関係しており、一方においては、1スレッドが通信している間にその他は計算することにより、通信と計算がオーバーラップすると仮定されます。実際には2番目の手法に性能利得は示されませんでした。最終的に、コードが簡潔であるという理由から最初の手法が選択されました。

       拘束力

DL_POLY_4の化学結合を記述やり方は拘束によるものです。2つの原子は常時一定の距離を保つように拘束されます。これはSHAKE[12](およびRATTLE[13])アルゴリズムで実装されています。この時、原子は最初に拘束条件無で動きますが、その後正しい配置に戻すために拘束を記述する非線形方程式がガウス・ザイデル法で繰り返し法によって解かれます。つまり拘束条件は2つの原子座標に影響を与えます。さらに各原子は複数の拘束条件を持ちます：例えば簡単な例として水分子を考えると、酸素分子は2つの水素分子との拘束条件を持ちます。

こうしたDL_POLY_4の構造から、最も単純なSHAKEルーチンの並列化単位は拘束条件であると言えます。しかしながらこれまでの議論から、原子の座標は複数の拘束条件から決定されるため、その更新時には競合が発生する可能性があります。この回避策として様々な手法を検討した結果、スレッド間の同期が生じますがatomic指示行が最も効果的であることが解りました。

実際にはSHAKEとRATTLEを、この同期が不要な原子で並列化することも可能です。しかしこれを行うためにはDL_POLY_4のかなり大きな修正が必要であり、その実施作業の時間がありませんでした。

       通信

時間ステップ毎のグリッド上のPDE計算の場合と同様に、DL_POLY_4にもプロセッサー上の情報を正しく更新するためのデータ通信が必要となります。DL_POLY_4の場合はその通信には、プロセッサー担当領域のハロ部の更新のみならず、原子が動くため、プログラム全体としての領域それ自身の更新が含まれます。

PDEの場合は、OpenMP/MPIハイブリッド並列プログラムにおいて、スレッドに対して通信と計算をオーバーラップさせることによりかなりの性能向上が示されています。例を図1に示します。これは2Dポアソン方程式の高速化を示しており、もし通信と計算のオーバーラップ(CCO)可能ならば、ほぼ線形にスケールします。

同様な効果がDL_POLY_4にもたらせられればよいが、これは現状のフレームワークの範囲内では困難です。結局これは、並列スレッドを用いてリンクセルを構築する場合に似通った状況です。つまりグリッドとは違い、通信すべき原子がどれなのかが事前にわからないので通信と計算のオーバーラップが困難となります。そのため原子のリストとデータは通信前に生成する必要が生じ、以前と同様に並列オーバーヘッドが発生してしまいます。事実リンクセルの場合は、領域内の全ての原子がリンクセルに含まれ、通信が必要なサブセットのみを扱う場合には、同期オーバーヘッドは悪化します。

通信段階の効率的な並列化に関して様々な手法を試しました。例えば、通信バッファーが送信可能になるサイズになったら時点で、他のスレッドが次のバッファーの構築を開始するやり方です。何れにしても元のアルゴリズムがより効果的であることが判明し、最終的にはユーザーの選択で並列通信を可能とし、デフォルトではシリアルコードが設定されました。しかしながらここで開発された手法は別のプロジェクトで組み込まれ、より一般化された方法でコード内の通信が改善されました。

図1　：　2Dラプラス方程式ソルバーの性能

2.3  結果

このセクションではOpenMP/MPIハイブリッド並列コードの性能についてレビューします。これら全てにおいて、横軸の「Cores」はMPIプロセス数でなく、利用したコア総数を意味します。よって256コアは、1スレッドを持つ256MPIプロセスか、2スレッドの128MPIプロセスのどちらかの場合を意味します。全ての場合においてノードは飽和させています。これらは全てGNUコンパイラー4.4以降のバージョンで作成された実行モジュールで、HECToRのCrayコンピューターおよび、制限がありますがintelおよびSUNコンパイラーによる他のマシンにおいても実行しました。

図2は、HECToRで実行した、最も単純な系の一つである液体アルゴンのシミュレーションです。力場はファン・デル・ワールス項のみで、12-6レナード-ジョーンズ・ポテンシャルでモデリングしました。系は256,000個の原子で構成され、ユニットセルは1辺が210.36Åの立方体、カットオフは9Åとしました。

図2　：　カットオフ9Åでのアルゴン系の性能

シングルスレッドのコードが全体でほぼ優位な性能を示しています。これはある意味当然の結果です。並列リンクセル法はユニットセルの辺より十分小さいカットオフを持つ短距離力に対する解法で、かつOpenMPはこの過程が破たんする場合のスケーリング性能の改善策でしたが、今回はそのケースではありません。
図3は、同じ系をカットオフ15Åで計算した結果です。

図3　：　カットオフ15Åでのアルゴン系の性能

ここでは、明らかにスレッディングされたコードがそのスケーラビリティを改善しています。実際に、青色のシングルスレッドの線は512コアで中断されていますが、これはリンクセル法がこれ以上計算不能となっており、長いカットオフにより各領域内にリンクセルが1つしか存在しないためです。しかしながらOpenMP化によるスレッディングにより、8スレッドを持つ256MPIプロセスを使用すれば2048コアまで効果的な利用が可能になりました。複数スレッドのスケーラビリティ改善の理由として、MPIプロセスが少なくなった結果より大きなメッセージを用いていることがあります。この結果、効率的なスケーリングを妨げる遅延領域で実行がトラップされることが避けられています。

図4　：　ナトリウム/カリウム・二ケイ酸塩ガラス系の性能

アルゴン系での力場は単純ですこし特殊なケースでしたが、図4はより現実的なケースの実例です。これはナトリウム/カリウム・二ケイ酸塩ガラスの場合です。シミュレーションは69,120個の粒子から成り、カットオフ12.03Åで、一辺が96.72Åの立方体セルを用いました。このとき力場は、塩化ナトリウムの場合でのVan Der WaalsやEwald項のみでなく3体項を含め、ほぼすべての重要な非結合項を用いています。ここでもスレッディングが利用可能なコア数を増加させています。以前と同様にシングルスレッドの線はMPIプロセスで利用可能なコア数により制限されています。今回はより現実的な問題であり、意味もなくカットオフを大きくしているわけではありません。

図5は最後の例として、水中の生体分子グラミシジン-Aの結果を示しています。これは抗生物質的な振る舞いを示すポリペプチド分子です。シミュレーションは総数99,120個の原子を含み、カットオフは8Åで、ユニットセルはa=94.6Åおよびc=112.7Åの正方晶形です。ここでの力場は先ほどのものよりも複雑で、上述の力場に加えてボンド結合を含みます。これは水分子に対しては分離拘束として、グラミシジン分子に対しては、バネ結合、結合角および二面角としてモデリングされます。

図5　：　水中の生体分子グラミシジン-A系の性能

このケースでは、性能は2スレッドまでに抑えられており、それより多重のスレッドは性能が減少しています。この理由として拘束条件が主たる原因であることが判明しました。表1は64MPＩプロセスでのスレッド数に対する性能を示しています。

Threads	Time/s
1	6.64
2	4.81
4	4.10
8	3.67

表1：拘束条件ルーチンのスケーリング性能

atomic指示行による同期の影響が明らかで、より良好なスレッド並列性能が必要であればSHAKEルーチンの改良をして、二ケイ酸塩ガラスに適したスレッド数でスケールするように直すべきでしょう。

よってこうしたスレッド化は一般には有利な結果となりますが、スレッドとプロセスを如何にうまく組み合わせるかは、力場および対象の系の性質に依存します。実際に最も効率の良いコードを決めるには、このような組み合わせの実験的な試行計算を幾つか繰り返して決める必要がありますが、経験上ではハイブリッド並列が有利です。

2.4  現状

上述の修正を含めたコードは現在、CCPForgeのDL_POLYリポジトリーのexperimental branchに有り、近い将来にこれらをmain branchへマージしてリリースする予定です。

2.5  結論

上述の例から明らかなように、スレッドの利用はDL_POLYで利用可能なコア数を増やして性能を改善します。先述のガラスのようないくつかの例においては、コア数と性能共に大きく増加させることが出来ます。達成可能なことはシミュレーションで用いる力場の形に依存していますが、スレッドの利用により、特にMPIプロセスによる制限に関して、DL_POLYの能力は明確に向上しました。

3  百万原子計算

3.1  イントロダクション

DL_POLY_4を用いた大規模シミュレーションは現在500万原子のオーダーです。これらは放射線核崩壊時の反跳による物質中のダメージを議論する際に利用されます。核崩壊で生じる超高エネルギーのために、これは物質中の広範囲に影響を及ぼします。よってこのシミュレーションには極端に大きなセルと膨大な原子数が要求されます。しかしながら現状のシミュレーションで用いている原子数では、その体積は実際の高エネルギー反跳を扱うにはまだ十分に大きくなく、現実的なダメージ現象のシミュレーションに単独の多重およびoverlapping反跳しか扱うことが出来ません。このためより大きな系を用いる事が望まれていますが、こうするにはコード中の単純かつ多数の制限が露わになります。それはコード全体を通して用いられているデフォルトの32bit整数で、これらは正の最大値が$2^{31}=2,147,483,648$であり、現状実行されているシミュレーションよりも高々4倍程度の大きさでしかありません。

よってコード全体で利用される整数インデックスは、より大きな整数種別へ変更しなくてはなりません。最新のFortran2008標準においては、$+/-10^{18}$の整数サポート[14]、およびこれらの整数のFortranコンパイラーでのサポート[15]が規定されており、コード内でグローバルインデックスとして用いる全ての変数に対し、この整数種別を用いる事が推奨されています。

実行中に特定の原子を識別するためにこのようなグローバルインデックスが用いられています。簡単な例としては出力内での、原子の運動をフォローするための各原子に関する様々な値があります。内部的には原子種の結合を定義するのに用いられる様々なテーブルもあります。つまり水分子に関する力を計算する際には、酸素原子が2つの水素原子と結合していることだけでなく、それがどの水素原子なのかを知らなければなりません。グローバルインデックスはコード全体に散りばめられているため、どれを変更し、どれを32bitのまま残すかを注意深く正確に識別していかなくてはなりません。DL_POLY_4の開発者であるIlian Todorov博士には不明瞭な点に関し助言をいただきました。事実こうした議論により、最終的なコードは、特にどの変数がコード内の異なる場所でローカルおよびグローバルインデックスが用いられているかについて、元のコードより明確になりました。

この仕事は、実行時の整数オーバーフローをチェック可能なnAGコンパイラー5.3.1以降のバージョンを援用しつつ、単純ですが骨の折れる仕事でした。さらに困ったことには、全ての整数を更新するのは得策でないことが判明しました。この主たる要因は、コード内の最大のデータ構造であるVerlet近傍リストで、これはローカルインデックスであり、32bit種別整数で十分なものです。64bitへ変更してしまうと、メモリーがほぼ倍になり、コア当たりのメモリーが少なくなるトレンドにおいては好ましくありません。

3.2  64bitコードテスト

64bitコードのテストが可能な最大のマシンはHECToRで、少なくともマシンリソースの半分が必要でした。実際にはいくつかのケースでは、あえてカットオフを小さくしなければならず、精密な実行が出来ませんでした。しかしながらそれでもコード内の分子動力学計算部分のテストには十分でした。また、ディスクスペースを考慮して多くのI/Oを止めなければならず、その部分については希望通りのテストが出来ませんでした。

また、ある部分については、HECToR上にあるMPIライブラリーの64bitの取り扱いにバグがあったため、順調には進みませんでした。例えば、ルーチンMPI_TYPE_CREATE_F90_INTEGERは、その後の通信で用いられる要求された範囲のMPI整数型へハンドルを返します。実際に、ハンドルが返された後の通信は、MPIライブラリーが要求する方が存在しないために中断されてしまいました。したがってMPI標準では不要ですが、MPI_INTEGER8を用いて修正を行いました。

基本的なコードチェックとしてアルゴン(かなり小さなカットオフを用いて)と、非結合項のチェックとして塩化ナトリウム(多少小さなカットオフで)および結合項のチェックとして水分子(多少小さなEwald許容差で)の、3つのケースを実行しました。各ケースでの検証は、結果の正確性と、スケール性能および原子当たりのエネルギーが所望の誤差範囲内で一定であるかという点を見ました。これまでと同じくGNUコンパイラーを用いました。
アルゴンの結果を表2に示します。

Number Of Particles	Number/2^31	Total Energy	Total Energy/Particle
2048000	0.0009536743	379597287015.312	185350.237800445
16384000	0.0076293945	3036778296122.86	185350.237800468
131072000	0.0610351563	24294226368983.3	185350.237800471
1048576000	0.48828125	194353810951866	185350.23780047
4294967296	2	796073209658845	185350.237800471

表2：アルゴン系

原子当たりのエネルギーは極めて良く保存されており、コードの正確性を証明しています。ちなみに最大のケースではシミュレーションセルの一辺は539nmであり、可視光波長と同等の長さです。

Number Of Particles	Number/2^31	Total Energy	Total Energy/Particle
4096000	0.0019073486	-159282467501.217	-38887.3211672893
32768000	0.0152587891	-1274643765558.87	-38899.040696987
262144000	0.1220703125	-10197290041763	-38899.5744390984
2097152000	0.9765625	-81566710329911.9	-38894.0383576927
4298942376	2.0018510409	-167173602476120	-38887.1466175987

表3：塩化ナトリウム

原子当たりのエネルギー誤差はアルゴン系に比べて増加しますが、ほぼ矛盾のない値と言えます。事実、Ewald和の許容誤差はほぼ5桁目に現れますが、これが示されており、コードの正確性がここでも確かめられました。
表4が最後の水分子の結果です。

Number Of Particles	Number/2^31	Total Energy	Total Energy/Particle
768	3.576278686*10^-7	-2484.9971734628	-3.2356734029
52931328	0.0246480703	-171254197.138575	-3.2354033728
578742528	0.2694979906	-1872529167.86996	-3.2355133367
2152873728	1.0025099516	-6965827487.98538	-3.2355950084

表4：水分子

ここでもほぼ5桁目の違いが示され、許容誤差のオーダーとなっています。
これらの結果は全てコードの正確性を支持しています。これらテストは十分ではないですが、力場、特に64bit化が影響を与える重要な部分のほとんどをカバーしています。さらに多くのテストが必用ですが、それを実施するにはHECToR上の32,000以上のコアの利用許可を必要とし、そのスケールでのデバックは極めて多くの時間を必要とするでしょう。

4  結論

作業は成功裏に完了しました。スレッドの利用は、DL_POLY_4のスケーラビリティを、5-10倍の利用可能コア数へ拡張し、40億原子までに制限されたシミュレーションでしたが64bit整数の導入を実施しました。こうした改善が取り込まれたDL_POLY_4の正式パッケージが近い将来にリリースされます。64bit化は現在はHECToR上でのデモンストレーション実行のみですが、次期ARCHERや他の高性能マシン上で、他の一般目的MDシミュレーションパッケージでは不可能なスケールのシミュレーションを可能にするでしょう。

5  謝辞

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

注記

[1]	http://www.ccp5.ac.uk/DL_POLY/
[2]	http://www.sci-techdaresbury.com/properties/daresbury-laboratory/
[3]	http://www.ccp5.ac.uk/
[4]	I.T. Todorov, W. Smith, K. Trachenko & M.T.Dove, J. Mater. Chem., 16, 1911-1918 (2006)
[5]	I.J. Bush, I.T. Todorov and W. Smith, Comp. Phys. Commun., 175, 323-329 (2006)
[6]	M.R.S. Pinches, D. Tildesley, W. Smith, 1991, Mol Simulation, 6, 51
[7]	Verlet, L. (1967). "Computer 'experiments' on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules". Phys. Rev. 159: 98-103.
[8]	Essmann U et al. "A smooth particle mesh Ewald method" J. Chem. Phys. 103 8577 (1995)
[9]	http://www.hector.ac.uk/cse/distributedcse/reports/DL_POLY03/DL_POLY03_domain/
[10]	C. Temperton, J. Comp. Phys. 52, 1, (1983)
11]	http://developer.amd.com/tools-and-sdks/cpu-development/amd-core-math-library-acm/
[12]	Ryckaert, J-P; Ciccotti G, Berendsen HJC (1977). "Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes". Journal of Computational Physics 23 (3): 327-341. Bibcode:1977JCoPh..23..327R. Doi:10.1016/0021-9991(77)90098-5.
[13]	Andersen, Hans C. (1983). "RATTLE: A "Velocity" Version of the SHAKE Algorithm for Molecular Dynamics Calculations". Journal of Computational Physics 52: 24-34. Bibcode:1983JCoPh..52...24A. doi:10.1016/0021-9991(83)90014-1.
[14]	"Modern Fortran Explained", Metcalf, Reid and Cohen, section 20.2.1
[15]	http://www.polyhedron.com/pb05-linux-language0html