＊ここに掲載するのは、STFC Rutherford Appleton LaboratoryのDominik B. Jochym博士によるHECToRレポート「Efficient massively-parallel tools for the study of catalytic chemistry A dCSE Project, Thomas W. Keal, STFC Daresbury Laboratory, 29th June 2010」を要約したものです。
概要
計算化学環境ツールChemShell [1,2]は、量子力学（QM）と分子力学（MM）ソフトウェアパッケージを統合してQM/MMを組み合わせた計算を行う手段を提供します。外部パッケージはエネルギーとグラジェントの計算のみに利用され、ChemShellのルーチンはQM/MMのグラジエントを形作り、座標最適化や分子動力学などの高度なタスクを処理するために使用されます。ChemShellの計算は、Tcl言語[3]スクリプトで制御され、ユーザから外部プログラムの詳細を可能な限り隠蔽して、一般的なコマンドが用いられます。
QM/MMアプローチは、不均一触媒の研究において特に有用です。ChemShellは、バルク材料または表面が有限MMクラスタモデルで表され、バルク静電気の影響を模倣するためにオプションで追加のポイントチャージがある、埋め込みクラスター計算[2,4]をサポートしています。アクティブサイトは静電埋め込みによる高レベルQM計算によってモデル化され、クラスタ環境はQMハミルトニアンの点電荷によって表現されます。この設定は、アクティブサイトの精度と計算コストの間で最適なバランスを取るようにします。
ChemShellは、現在のリリース（v3.3）には独自の並列アルゴリズムは持ちませんが、シリアルまたは並列プログラムとして構築することができます。並列バージョンでは、エネルギーとグラジェント計算に用いるパッケージに含まれる並列性を利用します。通常は、計算の評価ステップが計算コストが最も高いので、このステップは良くスケールします。しかしながら、何千ものコアを含む計算が日常的なHECToRなどの大型スーパーコンピュータでは、外部プログラムの並列アルゴリズムが効率的にスケールしないことがしばしばあります。この領域では、ChemShellレベルへ第2階層の並列処理を追加することで大きな利得が期待されます。
プロジェクト概要
このプロジェクトの目的は、材料のクラスタシミュレーションで一般的に採用されるタスクを並列化するために、ChemShellにタスクファーミング並列フレームワークを追加することです。2009年1月から2010年6月までの18ヶ月の期間にわたって、9カ月の開発が行われました。このプロジェクトのコード開発目標は以下の通りです：
- ChemShell（および関連するTclライブラリ）は、コミュニケータ分割により並列処理を行う。Tclインタプリタ内の条件テストを使用して特定のワークグループでコマンドが実行可能にする。
- ChemShell内の有限差分へシアン行列コードを修正してタスクファーミングを可能にする。既存の実装（1,000プロセッサを超える大規模なHECToRパーティションで実行される単一の並列タスク）とタスクファーミングバージョンを比較したときのテスト計算において2倍高速化。
- ChemShell内で別のタスクを行うDL-FINDの並列実装とのインタフェース実装を行う。表面化学問題における遷移状態決定のための並列NEBの実証において2倍高速化。
- DL-FIND並列最適化機能を、タスクファーミングバージョンにインターフェイスする。確率的探索アルゴリズムと遺伝的アルゴリズムを用いた大域的最小探索の実証計算を行う。ベースラインコード（ChemShellの非タスクファーミングバージョン）コードにおいて、確率的検索を3倍高速化。
以下に詳述するように、すべての目標が達成され、ベンチマーク目標を上回りました。
このレポートにおけるベンチマーク計算は、フェーズ2a構成のHECToRのXT4コンポーネントで実行されました。これは、1コアあたり2GBのメモリの2.3GHzのクアッドコアAMD Opteronプロセッサを搭載したCrayマシンです。ここで用いた並列化手法は単純な特性のため、最近導入されたHECToRのPhase 2b構成および他の並列プラットフォームでも効率的に実行されると期待されます。ここで用いたChemShell、GAMESS-UK、GULPは、Crayコンパイララッパーftnを用いたPGI CコンパイラおよびFortranコンパイラ（v9.0-4）によりHECToRでコンパイルされたものです。
タスクファーミング並列
タスクファーミング並列処理は独立した計算セットを実行する場合に効果があります。「タスクファーム」は利用可能な全てのプロセッサで構成され、ワークグループと呼ばれるプロセッサのサブセットに分割されます。タスクは、並列動作するワークグループ間で分割されます。これら計算は独立しているため、計算中にワークグループ間で情報を交換する必要はなく、全てのタスクが完了するまで結果の共有を待つことが可能です。ChemShellでタスク・ファーミングを実装するには、外部プログラムと部分的なコードの修正が必要でした。
· Chemshell
ChemShellはMPIを用いて並列実装されています。既存コードでは、1つのノードが「マスター」として機能し、Tclインタープリタが実行され、入力スクリプトが実行されます。他のノードは「スレーブ」で、並列コマンドが実行されるまでアイドルします。これが通常の外部プログラムを使用するグラジェント計算のやり方です。最大限の効率を達成するために外部プログラムはChemShell実行形式ファイルに直接リンクされており、追加プロセスの生成を必要とせず直接呼び出すことができます。したがって、外部プログラムはChemShellと同じMPI環境を共有し、すべてのノードを自身の並列計算に利用することができます。コマンドが完了すると、制御がマスターノードに戻ります。
タスクファーミング並列処理では、プロセッサにもう一つのレイヤーを追加します。これはMPIコミュニケータを用いた独立したワークグループにグループ分けされるます。既存の並列フレームワークでは、ノードはデフォルトのMPI_COMM_WORLDコミュニケータにグループ化されていました。独立したワークグループを生成するために、MPI_COMM_WORLDは、MPI_COMM_WORKGROUPと呼ばれるワークグループコミュニケータを持つ小さなプロセッサセットに分割されます。ChemShellの実行時には、コマンドライン引数nworkgroupsで生成するワークグループ数を指定します。
各ワークグループでは、1つのノードがマスタとして動作し、残りのノードはスレーブとして動作します。つまりマスターノードは全体で複数存在します。
· DL-FIND
DL-FINDは座標最適化ライブラリで[5,6]、ChemShellディストリビューションに含まれる標準最適化ドライバです。ChemShellとDL-FINDは、2つのコード間でオプション、座標、エネルギー、グラジエントを渡すために使用され、明確に定義されたインターフェースを介して通信します。
DL-FINDは独自のタスクファーミング並列実装を持ちます[6]。この並列処理は、プロセッサ間でデータを共有するためのMPIライブラリ呼出しのラッパー関数群と、呼出し元プログラムとの間でタスクファームの詳細を渡すためのパラレルインタフェースで構成されます。この並列化は2種の実装を持ちます：最初のものでは、呼出しプログラムはタスクファームをセットアップし、セットアップ情報とMPIコミュニケータをDL-FINDに渡し、2番目のやり方では、DL-FINDがタスクファームをセットアップし、セットアップ情報を返します。
ChemShellにリンクする場合には、ChemShellが初期化された際にタスクファーム環境を設定し、DLFINDがこの中で動作しなければならないため、最初の戦略が最も適切です。DL-FINDは、MPI呼び出しを使用してワークグループ間で直接通信することに注意してください。これは通常、各最適化サイクルの終わりにグラジェントデータを共有するために発生します。ChemShellによって提供されるMPI_COMM_COUNTERPARTSというコミュニケータはこの目的で使用され、マスターノードをグループ化するものです。
ChemShell環境をDL-FINDにアクセス可能にするために、インターフェイスサブルーチンをChemShellに追加し、タスクファームのセットアップに関する必要な情報を提供するようにしました。
· 外部プログラム
ChemShellインタプリタレベル以下の全ての並列作業は、ワークグループ内で実行する必要があります。したがって、外部プログラムがエネルギーおよびグラジェント計算で呼び出された場合、MPI_COMM_WORLDではなくMPI_COMM_WORKGROUPを使用する必要があります。よってChemShellにインターフェイス機能を追加し、ワークグループコミュニケータを外部プログラムに渡しました。
ChemShellからMPI_COMM_WORKGROUPコミュニケータを受け入れるように、GAMESS-UK [7,8] QMパッケージとGULP [9] MMパッケージを変更しました。これらはライブラリとしてChemShell実行形式ファイル内へコンパイルすることができ、単純な関数呼び出しを使用して情報を渡すことができます。
GAMESS-UKにおいては、行列固有値と固有ベクトル計算に多大な計算時間が掛かるため、大規模並列プラットフォームの利点を生かすにはPeIGS[10]といった並列ソルバの利用が重要となります。PeIGSライブラリもMPIを用いるため、タスクファーミング環境で実行可能にするために修正が必要でした。これにはMPI_COMM_WORLDのインスタンスをChemshellから渡されたワークグループコミュニケータへ置き換える作業を含みます。
原理的には如何なる並列プログラムも、MPI_COMM_WORLDの代わりにMPI_COMM_WORKGROUPを用いるようにすれば、タスクファーミング並列が可能です。また、ワークグループに1プロセッサのみが存在するように設定すれば、シリアルプログラムでも利用可能です。
並列性能
Chemshellの計算の中で3つのタイプを、タスクファーミング実装しました。これらは複数の独立したエネルギーあるいはグラジェント計算を含むため、ハイレベル並列化の利得を得られることが期待されます。
· 有限差分へシアン
へシアン行列の各要素は、2つのグラジェントの差から得られます。前進差分を用いたN原子系では、3N+1個の独立なグラジェント計算が必要です。中心差分では6Nです。
既存のChemshell実装では、グラジェントとへシアン計算は一つのTclコマンドで実行されます。これはタスクファーミング版では、3つのステージに分かれます。最初のステージでは、グラジェントが計算されてCehmshellオブジェクトとしてディスクに書き込まれます。この作業はワークグループに渡って分割されて並列に実行されます。第二ステージでは、Chemshellグラジェントオブジェクトはコマンドtaskfarm_globalise_forcegradientsにより全ワークグループで利用可能にされます。最後に、へシアン行列はこれら（既存の）グラジェントを用いて計算されます。必要であればこの計算を、1つのワークグループでの実行に制限することもできます。
57原子のケイ酸VO3クラスターを用いて、タスクファーミング性能を計測しました。エネルギーとグラジェントは、B3LYP汎関数[11]を用いたGAMESS-UKで計算しました。ここでは2種の基底関数を用いました：LANL2有効コアポテンシャル基底[12](413基底関数)および全電子はTZVP [13, 14](1032基底関数)。より大きな基底関数のテストはPeIGS版GAMESS-UKで可能です。PeIGSは全プロセッサ数より対角化する行列が大きい場合にのみ用います。
最初に、シングルポイントでのへシアン計算の並列性能を異なるプロセッサーに対して測定すると、近似的にスケールするのは128プロセッサまででした。
次に様々なワークグループ数に対して、前進差分へシアンを1024プロセッサで評価しました。負荷バランスは単純に、ワークグループに均等に同数のグラジェントを割当てるやり方です。
LANL2 ECP基底では、64ワークグループでの7倍が最も高速となり、GAMESS-UKの並列化のみの場合よりも効率的です。これはTZVP基底でも同様な結果が見られました。
· Nudged elastic band (NEB)法の性能
これは2つの構造間の最小エネルギーパスを見い出すためのアルゴリズムです。通常はエネルギー障壁を含む反応経路を分類するのに用いられます。DL-FINDではimproved-tangent variant of NEB[15]が利用可能です。
各NEB最適化サイクルは、2つのエンドポイント間のパスに沿って配置された一連の構造（イメージ）のエネルギーとグラジェントの計算から構成されます。最終的なNEB勾配は、イメージを接続するばね力を用いて構築されます。しかしながらイメージのグラジェント計算は独立しているため、並列に計算することが可能です。
DL-FINDのNEB法は静的な負荷バランス法で並列化されています。イメージは、ワークグループ数を法とするイメージ番号がワークグループIDと等しい様にアサインされます。この方法では常に特定のイメージは同じワークグループに割当てられることから、外部プログラムがリスタートだいるを用いる際に便利です。サイクルの最後には、DL-FINDのMPI関数により全ワークグループ内で、エネルギーとグラジェントが共有されます。
NEBサイクルの最初だけは、ワークグループは全パスに沿って順番にグラジェントを計算していきます。これは、QMプログラムの収束性のためで、前回のイメージの波動関数を次のイメージの初期値として用いるためです。引き継づくサイクルでは、各イメージにおいて前回の計算結果を初期値として用います。
アルミニウムがドープされた酸化亜鉛表面上に吸着した二酸化炭素と2つの水素分子から成る、3207原子のクラスターのQM/MM計算をベンチマークしました。2つの異なるサイトに在るH-CO2のHの交換に関するエネルギー障壁の計算に、NEB法を用います。QM領域はPVDZ基底関数(Zn以外。ここでStuttgart ECPを用いました[16])を用いた32原子から成ります。これは全194基底関数です。文献[4]に従って、これをQM計算対象の内部の19原子と、これを覆う擬ポテンシャルを持つ13個の亜鉛原子によるQM/MMインターフェイス領域に分割します。このインターフェイス領域はlocalised embedding potentialを与え、内部領域からMM領域の正電荷へ電子が移動するのを防ぐ役割を果たします。
QM計算には、B97-1汎関数を用いたGAMESS-UKを用います。GULPは、シェルモデル原子間ポテンシャルを用いたMMエネルギーとグラジェント計算に用いられます[17]。パスは10個
エネルギーとグラジェントの1点計算は、QMとMMの繰返しサイクル計算です。これは、QM領域がMM原子により分極し、逆にMM系のシェルはQM領域により分極するためです。よってQM/MMグラジェントは計算されるたびにこれが収束するまで繰り返されます。のイメージで表現され、両端を固定してサイクル当り8回のグラジェント計算を行います。
プロセッサを増やした際の性能は極めて貧弱です。1024プロセッサでの実行時間は256プロセッサの場合よりも大きくなります。これは、追加のプロセッサによる利得を上回る通信オーバーヘッドのためです。1点計算の全体に掛かる時間は極めて大きいですが、個々のQMやMM計算はそれほどでもなく、大規模プロセッサの利点を生かせません。QM領域を大きくすればその利得を得ることが出来ますが、この場合は実行コストが膨大になります。この結果からワークグループ数を可能な限り大きくすること(つまり、8に設定する)が示唆されます。
NEBベンチマーク計算を50サイクル実施しました。
1ワークグループ計算と比較して、1024プロセッサの場合は8倍以上の高速化が示されました。これは128プロセッサでの1点計算が1024プロセッサの場合よりも高速なためで、そうした場合のみ生じます。
· 大域的最適化法
DL-FINDでは2つの並列最適化法を選択できます[6]。即ち、遺伝的アルゴリズムと確率的探索法です。これらは一般に、ポテンシャルエネルギー曲面上の大域的最小値を求める際に用いられます。また共に各サイクル中の構造の分布上の計算を含み、これらは並列に実行可能です。 これらはNEB法の並列処理のベースとして用いられます。
DL-FINDは原子座標のみを用い、シェルを無視します。これはChemShell特有のものです。以前はシェル座標も絶対座標として保存されて、座標の更新に用いられていました。この場合に座標が大きく変化すると収束が遅く、あるいは困難になり得ます。これは大域的最適化に起こりがちな現象です。ChemShellにおいてシェルは相対座標として保存し、座標が大きく変化しても親原子の近傍に置かれるように改良されました。この改良はDL-FINDの他の最適化手法にも利点がありますが、特に大域的最適化に重要です。
ベンチマークは文献[19]に従い、ZnOナノクラスターを用いました。そこではベンチマーク用に固定したイオンのMM計算が使用され、より計算負荷の高いQM計算はGAMESS-UKが用いられました。クラスターは(ZnO)28です。B97-1汎関数を用いて、PVDZ基底関数(560基底)とZn原子にStuttgartECPを用いました。ここで構造は32個を用いました。これはこうした手法での典型的なサイズであり、かつタスクファーミング並列に十分な量です。
遺伝的アルゴリズムと確率的手法は似たスケール性を持つことから、確率的探索をベンチマークに使用しました。一つのエネルギーとグラジェントの計算は極めて速く終わってしまうので、十分な計測時間を得るためにフルサイクル計算を行いました。結果は256プロセッサが最高速で、プロセッサ数が増えると遅くなります。
これらの手法は数百数千サイクル実行することで、大域的最小を見つける機会が増します。しかしながら計算コストを鑑みれば、一つのワークグループでこうした計算を行うのは非現実的です。実際にはスケール性を見るために、20サイクルのベンチマークを実行しました。
結果として、16ワークグループを用いるとほぼ最高速で5倍以上に高速化されます。
· 全体性能
全ての手法にタスクファーミング並列が実装され、最も高性能なワークグループ数（負荷分散が発生する手前の限界値）が求まりました。これは当然のことながら、エネルギーやグラジェントの並列処理とは異なり、ChemShellの並列処理ではオーバーヘッドが極めて小さいためです。単一配置での計算は最終的な性能の良い指標たりえます。これを用いてその先のワークグループの最適な数を決定することが推奨されます。
謝辞
このプロジェクトは、NAG Ltd.が運営するHECToRの分散計算科学および工学(CSE)サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびNAG Ltd.のCSEサポートサービスにより管理運営されています。
文献
| [1] | ChemShell, a computational chemistry shell. See www.chemshell.org |
| [2] | Sherwood P, de?Vries AH, Guest MF, Schreckenbach G, Catlow CRA, French SA, Sokol AA, Bromley ST, Thiel W, Turner AJ, Billeter S, Terstegen F, Thiel S, Kendrick J, Rogers SC, Casci J, Watson M, King F, Karlsen E, Sjovoll M, Fahmi A, Schafer A, Lennartz C (2003) J Mol Struct (THEOCHEM) 632:1-28 |
| [3] | Tcl, the tool command language. See www.tcl.tk |
| [4] | Sokol AA, Bromley ST, French SA, Catlow CRA, Sherwood P (2004) Int J Quantum Chem 99:695-712 |
| [5] | DL-FIND, a geometry optimiser for quantum chemical and QM/MM calculations. See ccpforge.cse.rl.ac.uk/projects/dl-find/ |
| [6] | Kästner J, Carr JM, Keal TW, Thiel W, Wander A, Sherwood P (2009) J Phys Chem A 113:11856-11865 |
| [7] | GAMESS-UK, a package of ab initio programs. See http://www.cfs.dl.ac.uk/gamess-uk/ |
| [8] | Guest MF, Bush IJ, van Dam HJJ, Sherwood P, Thomas JMH, van Lenthe JH, Havenith RWA, Kendrick J (2005) Mol Phys 103:719-747 |
| [9] | Gale JD, Rohl AL (2003) Mol Simul 29:291-341 |
| [10] | The PeIGS parallel eigensolver library. See http://www.emsl.pnl.gov/docs/global/peigs.html |
| [11] | Stephens PJ, Devlin FJ, Chabalowski CF, Frisch MJ (1994) J Phys Chem 98:11623-11627 |
| [12] | Hay PJ, Wadt WR (1985) J Chem Phys 82:270-283 |
| [13] | Dunning TH (1971) J Chem Phys 55:716-723 |
| [14] | McLean AD, Chandler GS (1980) J Chem Phys 72:5639-5648 |
| [15] | Henkelman G, Jonsson H (2000) J Chem Phys 113:9978-9985 |
| [16] | Igel-Mann G ECP28SDF. See http://www.theochem.uni-stuttgart.de/pseudopotentials/clickpse.en.html |
| [17] | Whitmore L, Sokol AA, Catlow CRA (2002) Surf Sci 498:135-146 |
| [18] | Henkelman G, Uberuaga BP, Jonsson H (2000) J Chem Phys 113:9901-9904 |
| [19] | Al-Sunaidi AA, Sokol AA, Catlow CRA, Woodley SM (2008) J Phys Chem C 112:18860-18875 |