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

＊ここに掲載するのは、STFC Daresbury LaboratoryのThomas W. Keal博士によるHECToRレポート「Microiterative QM/MM Optimisation for Materials Chemistry A dCSE Project, Thomas W. Keal, STFC Daresbury Laboratory, 31st May 2013」を要約したものです。

概要

ChemShell[1,2]は、量子力学計算$QM$と分子力学計算$MM$を組み合わせた計算に用いられる計算科学環境ツールです。QM/MM計算においては、正確な（だが高計コストの）QM法が系の反応サイトを扱うのに用いられ、一方、電子構造の記述が不要な低コストな古典的力場がそれを取り巻く環境を表現するのに用いられます。このやり方により、全てをQMで計算する場合に比べて、全体の正確性を損なうことなく全計算コストを削減できます

QM/MM法は不均一触媒の研究にとって特に有用です。ChemShellはembedded cluster計算をサポートしており、バルクな物質または表面を有限のMMモデルで表現して、バルクな静電性を模した点電荷を追加することも可能です。反応サイトは、環境による静電ポテンシャル効果を組み込んだQM計算でモデル化され、QMハミルトニアン内の点電荷としてクラスター環境が表現されています。

HECToRシステムにおいては通常、QM計算にはGAMESS-UK[4,5]が良く用いられ、MM計算にはGULP[6]が用いられています。これらは共に、並列ライブラリとしてChemShellに直接リンクされます。ChemShell内の座標最適化ライブラリDLFIND[7,8]は、反応物や、生成物、反応経路、遷移状態などを最適化することで表面上の反応を特徴づけるのに用いられます。

このプロジェクトの目的は、DLFINDへmicroiterative法を実装にすることにより、大規模座標最適化におけるChemShellの性能を改善することです。QM/MMによる最適化の計算時間は、通常は量子力学計算に支配され（この時数百原子を扱うとすると）、周囲のMM領域は通常はさらに大きな原子数（数千原子）を扱います。通常の座標最適化においてはは、その各ステップにおいて新たな座標におけるQMとMM計算が必要です。つまり、QMおよびMM領域の座標は同期して緩和される必要があります。Microiterative法[9]では、系はQM原子を含む内側領域と残り系を含む外側領域に分割されます。この時、内側領域の各最適化ステップ（macroiterativeサイクル）毎に、外側領域が最適化されます（microiterativeサイクル）。この最適化により、外側領域のMM計算回数は増加しますがQM計算の数を削減することが可能です。MM計算は通常は極めて低コストなため、こうすることで全計算時間を削減できます。

DLFINDに実装されているmicroiterative QM/MM最適化手法は、ChemShellのHDLCOptモジュールにおける古い実装を踏襲しています。HDLCOptは生体系に特化した設計がされており、物性化学には適しません。これは多くがDL-FINDに置き換えられましたが、microiterative法が不足しています。よってこのプロジェクトでは、HDLCOpt機能のDL-FINDへの移行を説明し、DL-FIND計算の高速化を行います。無論DL-FINDにはmicroiterative法も実装します。

2012年4月から2013年3月までのプロジェクト期間において9人月が開発に要しました。開発内容は以下の通りです：

• DLFINDにmicroiterative法実装、microiterative法、すなわちmicroiterativeおよびmacroiterativeサイクルの導入ためのフレームワーク実装。
• microiterative法の導入のためのChemShell/DLFINDインターフェイスの修正。
• 外側領域のエネルギー最小化へのmicroiterative法実装。
• 外側領域の最小化時の内側領域の鞍部点への最適化における、P-RFOおよびダイマー法に依る遷移状態のmicroiterative最適化実装。
• microiterative反応経路最適化の実装。ここで、各イメージの環境を緩和させる際の、nudged elastic band法（経路に沿った複数の座標を含む）を内側領域へ適用する。

これらは全てHECToRフェーズ3システムへ実装され、LiをドープしたMgO上の水素の解離についてテストされました。

実装

·エネルギー最小化

DL-FINDのmicroiterative最小化法の実装は、既存のHDLCOptに従いました。
既存のHDLCOptにおいては、macroiterativeステップエネルギー（即ち、内側領域ステップ前の正確なエネルギー）は、最初のmicroiterativeステップの検証のためのリファレンスエネルギーとして利用されていました。エネルギー最小化においてはmicroiterativeおよびmacroiterativeステップ共にエネルギーが減少していくため、このやり方は適しています。しかしながら、DL-FINDの実装はmacroiterativeステップで必ずしもエネルギーが減少しないケースに対して設計されるため、外側領域計算ステップの前にリファレンスとしてESP（静電ポテンシャル：QM領域ので電子密度を点電荷で近似したもの）による近似的なエネルギーを計算します。

Macroiterativeステップにおいては、ChemShellのQM/MMドライバを呼び出す通常の最適化が実行され、エネルギーE0とグラジェントg0が出力されます。この時QM計算では、QM領域周囲のファンデルワールス表面上の静電ポテンシャルが計算されます。これによりChemShellの最小二乗フィッティングコードを用いて、QM原子上に点電荷がフィッティングされてESPが生成されます。こうして、この点電荷表現を用いたQM/MMエネルギーとグラジェント計算のために再度ChemShellが呼び出され、E⁰_ESPおよびg⁰_ESPが算出されます。

文献[9]に従って、補正された現在のエネルギーE1とグラジェントg1を計算します。この補正は、microiterativeサイクルとmacroiterativeサイクル間を滑らかにつなぐものです

DL-FINDは複数の省メモリーBFGS$LBFGS$を起動可能で、そのうちの2つはmicroiterative最小化計算に用いられます。最初のLBFGSはmacroiterativeステップで用いられ、全座標が用いられますが、内側領域のステップのみが対象です。これが全座標を用いる理由は、内側領域のステップでは全原子のグラジェントが考慮されるからです。よって、microiterativeサイクルで外側領域は最適化を進めるのに完全に収束する必要がありません。第二のLBFGSは、エネルギーのmicroiterativeループの最初で初期化され、外側領域の座標のみを用います。補正されたQM/MMエネルギーとグラジェントは、内側領域を固定したまま環境を緩和させるのに用いられます。このLBFGSデータはmicroiterative法モジュール内に別の配列として保存されます。

Microiterative最小化は、古典力場で記述された水中に溶解したQMで扱われるグリシン分子のケースで正確性をテストされました。また、既存のHDLCOptと同等の性能改善を示すことを検査するために、文献[9]の水球ケースをテストしました。これは平衡した水バルクから切り出した半径25Áの水球です。中央の3つの水分子をQM領域のmicroiterative内部領域としました。QM計算は、AM1ハミルトニアンを用いた半経験的QMパッケージMNDOで実行しました。半径5,10,15,20Áの活性領域を定義し、その中の水分子をTIP3P力場で記述しました。その外側は最適化計算では固定しました。

テストにおいては、カルテシアン座標を用い（結合の拘束は無し）ました。よって自由度の数は文面[9]よりも多くなっているため、最適化ステップ数は直接比較ができませんが、似た傾向が見出されました。以前のHDLCOptテストのように、最適化は複雑なエネルギー表面上で多くの接近した最小値上で行われたため、実行によっては異なる経路をとり、若干異なる収束エネルギーに達する可能性があります。サイクル数は自由度の数に単純にQM計算回数は5から12倍減少しました。この小さな経験的なQM計算では全体の実行時間に影響を与えませんが、QM計算が支配するような実用計算では大きな削減が期待できます。

·遷移状態の最適化

P-RFO
DL-FINDのP-RFO最適化では全ヘシアン計算が必要なため、自由度の多い系には適しません。HDLCOptの実装は柔軟なやり方が可能で、環境をLBFGSで緩和しつつコア領域を鞍部点へ最適化することが可能です。このやり方はmicroiterative最小化に似た方法であり、各P-RFOステップ毎の後に環境のグラジェントをゼロに最適化しなくてはなりません[10]。
DL-FINDでもmicroiterative　P-RFO最適化は、P-RFOで最適化される内側領域とLBFGSで最適化される外側領域に系を分割します。外側領域はmicroiterative最小化と同じやり方で最適化されます。
P-RFOおよびヘシアン計算ルーチンは修正され、内側領域で定義されたサブセットを対象にします。この点が、macroiterative最適化が全系を対象とする最小化スキームと異なる点です。なぜなら、遷移状態の最適化では内側領域と外側領域が別々に最適化されるためです。これは環境が各ステップ後にフルに最適化されるP-RFOケースでは重要な点で、特に最適化の終了に近い時にはこうしないとmacroiterationステップは正確性を失うことになります[10]。
P-RFOは、系の併進と回転に相当する‘ソフトモード’を検出し、これを排除することで収束性に影響を与えないようにします。しかしながら、これはmicroiterative法の最適化では適切ではなく、外側領域に対する内側領域の並進および回転は正確な収束に必須です。

テストは最小化の場合と同じくグリシン分子に対して行われ、グリシンとその創生イオン型間の遷移状態を最適化しました。
QM領域はグリシン分子で構成され、AM1ハミルトニアンを用いたMNDOで計算しました。内側領域（macroiterative）はQM領域と同じものとして定義しました。外側領域はTIP3P力場の水分子62個で構成され、そのうち13個を活性としました。また、ESPフィッティングをmicroiterativeで用いました。
標準的な遷移状態最適化とmicroiterative法によるものは直接比較は難しいですが（何故なら、標準的な手法では外側領域も鞍部点最適化対象だからです）、内側領域を適切に選べば両者はよい一致を示します。これは遷移状態が一般によく局所化しているためです。溶解したグリシンケースでは、両者のエネルギーは極めて良い一致を示しました。

ダイマー法
ダイマー法[11]はDL-FINDにおける別の遷移状態最適化手法で、ヘシア計算が不要なため大規模系に適します。ダイマー最適化では、固定した距離分離れた場所にあるポテンシャルエネルギー表面上の接近した2点（ダイマー）の計算を行います。最適化サイクル毎に、ダイマーは最初に末端の2点のエネルギーの和が最小になるようにその中間点に関して回転させられ、最小ヘシアン固有値モードに沿って配置されます。次にダイマーを、このモードに沿ってエネルギーが最大になるよう、かつその方向に垂直な方向にエネルギーが最小になるように移動させます。
DL-FINDに実装された標準ダイマー法最適化では、ダイマーの向きと回転力の計算時に個々の座標に重みをつけることが出来ます。重みをゼロにすると、座標は回転ステップで無視されて並進ステップのみで最小化されます。この排除された座標は最適化時に緩和させられる環境のように振る舞います。
よって、重み付けのやり方をダイマー法最適化のmicroiterative形式へ拡張することは自然な考え方です。そこで、microiterativeP-RFOと同様に系を分割し、座標リストを最初にイメージ、次に領域により順序付けます。
外側領域座標については遷移状態最適化に考慮する必要がないため、重みをゼロにします。Macroiterativeサイクルは、内側領域のダイマー回転とその後の並進ステップから構成されます。Microiterativeサイクルは、microiterative P-RFOと同じく環境を緩和します。LBFGSを用いる最適化は3つ（ダイマー回転、ダイマー並進、環境エネルギー最小化）です。環境エネルギー最小化は、ダイマーの中央点の移動に相当するため、microiterativeサイクルではこの座標のみを用います。

Microiterativeダイマー法は前回と同様に、水に溶解したグリシンの遷移状態でテストしました。同様なQMとMMのセットアップを行い、ESPフィッティングをmicroiterativeサイクルで用いました。結果は、標準ダイマー法とmicroiterative法、およびmicroiterativeP-RFOの間で極めて良い一致を示しました。

·反応経路の最適化

DL-FINDに実装されているnudged elastic band $NEB$法[12,13]は、既知の反応物と生成物の座標間の反応経路を最適化します。NEB最適解においては、末端間の初期パスを推定し（線形に推定あるいはユーザ指定の座標のどちらか）、これに沿って複数の座標イメージを定義します。各イメージのエネルギーとグラジェントは独立に計算され、全NEBグラジェントが計算されます。これは、各イメージ間の局所的な接線上に働くバネ力とそれに垂直に働く真の力から構成されます。NEBグラジェントを最適化することで、イメージは最小エネルギー反応経路を表現します。

DL-FINDのNEB実装では、ダイマー法での陽に座標も重み付けが可能です。重みをゼロにすると、その座標にはバネ力は計算されず、単純に最小化されます。こうしてNEBの内側領域とは独立して、環境領域を定義できます。

Microiterative NEB最適化では、この重み付け手段への拡張が実装されました。ダイマー法のように系を、NEBグラジェントを計算する内側領域と重みをゼロにした外側領域に分割します。Macroiterativeステップにおいて、内側領域はNEBグラジェントに従って移動し、続くmicroiterativeサイクルでは外側領域が緩和されます列。外側領域のグラジェントは全イメージに対して定義されますが、実質的には独立した一連の最小化です（固定された内側領域を通して間接的に連動します）。しかしながらDL-FINDの実装では、一連の外側領域グラジェントは合わせられて、1つのLBFGSによる一つの最適化として扱うことが出来ます。microiterativeサイクルでは、RMS、最大ステップ、グラジェント収束条件は通常に計算されて、全イメージの平均エネルギーがエネルギー収束条件として用いられます。

Microiterative法も、frozenイメージとclimbingイメージを扱う必要があります。イメージは、それにかかる力が指定の閾値以下になれば固定されます。ここでこのイメージの外側領域の座標も固定されます。Climbingイメージは、他の外側領域イメージに関して最小化されたもう一つの外側領域座標のセットを導入します。

ESPフィッティングはNEB法ではさらに複雑です。Microiterativeサイクルを正確に実行するためには、macroiterativeステップの中で各NEBイメージに対してESPフィッティングを実行しなくてはなりません。これらは、イメージ番号を付与された別のChemShellフラグメントに保存され、microiterativeグラジェント計算の中で適切に呼び出されなければなりません。

テストとして、溶解したグリシンを用いて、最小化された標準および双性イオン性座標間の反応経路の最適化を行いました。他の手法と同じようにQMをセットアップします。イメージ固定とclimbingイメージの両方を用いました。しかしながらこれまでのように標準的最適化とmicroiterative法を直接比較することは困難です。何故なら、microiterative環境が最小化されているためで、最適化されたclimbingイメージのエネルギーと座標は、microiterative最適化が標準的な最適化と同じ点に収束していることを示します。
環境でNEBのバネ力を排除することにより注意すべき点は、各環境イメージ最小化が独立しているため、環境によって反応物から生成物に至る経路が最適化中に不連続になるというリスクを伴うことです。この問題を回避する方法がXie et al. [14]により提案されていますがこれは高コストです。今回のテストケースではこの問題が生じるケースは確認されなかったためこの方法を採用しませんでした。

·シェルモデルの最適化

通常の静電場埋め込み型のelectrostatic embedding QM/MM計算では、古典MM原子はQMハミルトニアンへ点電荷として影響を与えてQM領域を分極させますが、QM原子はMM領域を分極させません。シェルモデル力場はMM領域を分極させることが可能です。各MM原子はコアとシェルから構成され、これらはバネ力で結合します。MM電荷の一部はシェルへ割り当てられ、QM領域の静電ポテンシャルから影響されて動くことが可能です。逆にQM領域の分極もシェルの位置により影響を受けるため、QM/MMエネルギーとグラジェントの計算では、その状態がセルフコンシステントになるまでQMとMM計算を繰り返すことになります。

Microiterative法の観点からは、DL-FINDに特別な変更は不要です。シェルは最適化手法としては意識されず、最終的なQM/MMエネルギーとグラジェントのみが要求されます。ChemShellは、ESPフィッティングによるmicroiterativeシェルモデル最適化をサポートするように変更されました。完全なQM / MMエネルギーおよびグラジェント計算を行う場合は常に、QMおよびMM領域の分極はセルフコンシステントになるまで繰り返されます。これはmacroiterativeケースで、ESPフィッティングが使用されない場合は（QM領域が内側領域にある場合は固定されるが、MM領域によるQM領域の分極は固定されない）microiterativeケースでもあります。 しかしながら、QM領域がESPフィッティング電荷によって近似される場合、QM領域の分極は、microiteration中は固定されます。次にESP電荷はGULPに直接渡され、シェルは一回のGULP計算で緩和されます（反復する必要がある力そのものを供給するのではなく、シェルが動くときに力がどのように変化するかを考慮します）。

テストはHDLcOpt実装時[9]のMgO欠損例で行い、同じ最小値に収束することが確認されました。
性能測定のために、LiがドープされたMgO上の水素の解離反応を用いました。この系は6349原子のクラスターで、クラスター周囲に87個の点電荷を配置してバルク物質の静電性を模しました。QM領域はクラスター中心の33原子で、水素分子とLiドーパント、酸素原子5個、マグネシウム原子25個を含みます。このマグネシウムの内、4個をQM計算し、残りをQM周囲の境界に配置して擬ポテンシャルでモデル化して電子密度が周囲に漏れるのを防ぎます。834個の最適化対象の活性原子を定義しました。Microiterative計算では、内側領域は（境界の原子を含めて）QM領域と同じものを定義しました。
QM計算は、Def2 TZVP基底関数$274 球面調和基底関数$とB97-2汎関数を用いてGAMESS-UKで行いました。多重度2の非制限計算を行いました。MM計算は、シェルモデルポテンシャルを用いてGULPで行いました。
計算では、ESPフィッティングにおいて、多くのQM原子に非物理的な電荷が割り当てられたため、ESPフィッティングなしでテストを行いました。

テストは、反応末端の物理吸着および解離した水素のエネルギー最小化、P-RFOおよびダイマー法による遷移状態最適化、全反応経路のNEB最適化から構成されます。エネルギー最小化では、Macroiterativeサイクル数は標準の最適化に比べ劇的に減少しました（約10倍）。これは初期構造の良さと系が本来持つ性質によるもので、より複雑なエネルギー面ではこのレベルの改善は見られないと考えられます。
遷移状態最適化では、ダイマー法のmacroiterativeサイクル数は標準手法よりも減少しました。このケースでは、P-RFO法はダイマー法のmacroiterativeおよびmicroiterativeサイクル数を下回りました。
NEBにおいては、10個のイメージを用いて（freezing無し）比較的緩やかな収束条件を課しました。また10ワークグループによるタスクファーミング並列を行いました。ダイマー法と同様なサイクル数の減少が観測されましたが、NEBの場合はその数はより小さなもので、初期パスが良好だったことが示唆されます。

謝辞

このプロジェクトは、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]	Sokol AA, Bromley ST, French SA, Catlow CRA, Sherwood P $2004$ Int J Quantum Chem 99:695-712
[4]	GAMESS-UK, a package of ab initio programs. See http://www.cfs.dl.ac.uk/gamess-uk/
[5]	Guest MF, Bush IJ, van Dam HJJ, Sherwood P, Thomas JMH, van Lenthe JH, Havenith RWA, Kendrick J $2005$ Mol Phys 103:719-747
[6]	Gale JD, Rohl AL $2003$ Mol Simul 29:291-341
[7]	DL-FIND, a geometry optimiser for quantum chemical and QM/MM calculations. See ccpforge.cse.rl.ac.uk/projects/dl-find/
[8]	Kästner J, Carr JM, Keal TW, Thiel W, Wander A, Sherwood P $2009$, J Phys Chem A 113:11856-11865
[9]	Kästner J, Thiel S, Senn HM, Sherwood P, Thiel W $2007$ J Chem Theory Comput 3:1064-1072
[10]	Turner AJ, Moliner V, Williams IH $1999$ Phys Chem Chem Phys, 1:1323-1331
[11]	Kästner J, Sherwood P $2008$ J Chem Phys 128:014106
[12]	Henkelman G, J´onsson H $2000$ J Chem Phys 113:9978-9985
[13]	Henkelman G, Uberuaga BP, J´onsson H $2000$ J Chem Phys 113:9901-9904
[14]	Xie L, Liu H, Yang W $2004$ J Chem Phys 120:8039-8052