Gromacs / Plumed 活用・改造ガイド

現代において、MDソフトウェアを自作することは全くコストに見合わなくなりつつあります。シミュレーションソフトウェアに要求される内容（ポテンシャル関数等）が増えた一方で、ハードウェアが複雑化したため性能を出す為に高度なプログラミングが要求されるようになりました。しかしながら、MDシミュレーションの新しい手法を模索する場合、しばしばMDソフトウェアに手を入れる必要が出てきます。このページでは、GromacsとPlumedを改造する（あるいは改造なしでどうにかする）ことで、MDシミュレーション手法に手を入れたい人のためのガイドを説明します。読者としてGromacsでのMDシミュレーションをすでにやったことがある、学生さんやポスドクさんを想定しています。

GromacsかPlumedか

改造を行う際に考えておく必要があるのは、「その手法はポテンシャルの変更だけで実装できるか？」と「その手法を活用するにはトップスピードが必要か？」の2点です。

• ポテンシャルの変更だけで実装できる場合→Plumed

PlumedはGromacsなどのMDソフトウェアに対し、何らかの追加のポテンシャルを追加する時に汎用的に使えるソフトウェアです。実際の所、MDシミュレーションの手法のかなりの部分はこの方法で実装可能です。Plumedは様々な反応座標を用いることができ、さらに自動微分を用いて反応座標の関数などをポテンシャル計算に用いることができます。ポテンシャル関数への変更は、ほぼすべてのケースでこのPlumedをただ使うだけで実現でき、それが駄目な場合でもPlumedに少しの実装を加えることで実現可能です。

• 速度の問題があるときや、積分子など非ポテンシャル部分に手を入れたいとき→Gromacs

Plumedはポテンシャル以外の部分を変更することを苦手としています。ポテンシャル関数としての記述が難しく、MDソフトウェアの外側を大きく弄る必要がある場合、Gromacs本体を修正する必要があります。

また、Plumedを使う場合、どうしても速度が遅くなります。たとえば全近接原子間に相互作用を加える、といった操作をすると、Plumedでこれを実現するには全原子間に相互作用を入れる必要があり、毎ステップO(N^2)の時間が掛かり実用上使えません¹。速度を落とすこと無く、多数の原子に相互作用を追加したければ、Gromacsなどのソフトウェアに手を入れるのが適切です。

なお研究上は、速度が多少遅くなるだけであればまずPlumedでproof-of-conceptを実装し、MDの現実と理想との様々なギャップ(例: PBC)に依る問題が起こらないことを確認してからGromacs本体の改造に移ることをおすすめします。

Plumed

Plumed活用ガイド

Plumedができることは、PlumedのTutorialをざっと見て頂くのが早いと思います。Plumedではちょっとファイルに書き足すだけで、MDシミュレーションに対して追加のポテンシャル関数を容易に加えることが可能です。

• ポテンシャル関数は、原子座標から直接計算されるもの以外に、複数原子の重心、角度、二面角、といった値についての関数を使うことが出来ます。
• • Plumed内部ではこれは自動微分によって実現されます。計算する際には、まず原子座標から反応座標(角度、重心間距離、等)を計算し、しかる後にポテンシャル関数とその反応座標偏微分(反応座標方向の力)を計算します(自動微分のforward計算に相当)。次に、この反応座標での偏微分をもとに、chain ruleを用いて元の原子に掛かる力を計算します。
• ユーザーはf(g(x))のように、関数を自由にネストさせることが出来ます。たとえば、「タンパク質の3つのドメインのそれぞれの重心がなす角度」のようなものも指定できます。
• 使える関数は非常に多岐にわたっており、代表的なものは大方実装されています。さらに、簡単な関数であれば表引きやソースコードの修正無くユーザーが外から与える方法、PyTorchのへのインターフェースが存在します。

Plumed改造ガイド

ソースコードの修正が必要になりますが、ユーザーはPlumedに新たな反応座標を追加したり、新たな関数を追加したりすることが出来ます。それぞれのケースを簡単にまとめます。

• 反応座標の追加(執筆中)
• 関数の追加(執筆中)

Gromacs

Plumedでは難しい、となった場合、まずは既存のGromacsの表引き等のメカニズムで対処できるか、不可能である場合はさらにどういうメカニズムを足すかを考慮する必要があります。

Gromacs活用ガイド

Gromacsでもポテンシャル関数をある程度自由に弄ることが出来ます。Plumedに比べて自由度は下がりますが、自由エネルギー摂動計算などPlumedの対応に制約がある系を扱う場合、重要な選択肢となります。

Gromacsに入れられるユーザー定義ポテンシャル関数の話をするためには、まずGromacsのポテンシャル関数の分類を説明する必要があります。Gromacsのポテンシャル関数には以下の分類が存在します。

• listed-force

“listed-force” は[ bonds ], [ angles ]などで使われる、近距離相互作用の総称です。topology fileに最初に列挙されたlistだけを使って相互作用を計算し、勝手に増えたり減ったりはしません。制約として、空間分割を用いている場合、同じ空間分割のセルに2つの粒子が存在するか、あるいは隣接する2つの空間分割にまたがっている状況でなければなりません。それ以上離れた場合はエラーになります。

• nonbond

近距離でのみ働く相互作用の総称です。最も一般的には、Lennard-Jones相互作用と電荷相互作用(PMEの場合、Ewald和の実空間部分)がこれに当たります。相互作用が働くかどうかは動的に決定されます - 内部的にはnstlistステップごと²に粒子が相互作用しうる距離に存在するかを記述するneighbor listを更新し、 rlist以下の粒子ペア³をリストアップします。リストアップされた粒子間の相互作用ポテンシャルは各ステップ計算されて足されます。この機構から、rcoulombやrvdwの距離でポテンシャルが0になるような関数系にすることが望ましくなっています(そうでないとリストに入るかどうかで値が変わってくるため)。

Gromacsには、次のような方法でユーザー指定の関数をポテンシャルに入れる方法が存在します。

• [ bonds ]のfunctype = 8または9で、tabulated bondsを使うと、任意のポテンシャル関数がテーブル引きで実装できます。functype = 8と9の違いは、原子間のexclusion(nonbond相互作用の除去)を行うのがfunctype=8です。使い方はこんな感じです:

; (.topファイル)

[ bonds ]
;atoma atomb func tab
     1     2    8   0 ; 原子1と2の間にtabulated bond potentialを入れる。
; 原子1-2の間のnon-bond相互作用は切る
; table number は0→bondのbを合わせて(なんとか)-b0.xvgがポテンシャル関数となる

# table-b0.xvg
# x(nm), U(x), F(x) = -dU/dx(x) を書く。U(x)はkJ/mol, F(x)はkJ/mol/nm の次元になる。xは等間隔にすること
0.00    0.00   0.0
0.01    0.01   -2.0
0.02    0.04   -4.0
...
0.10    1.00   -20.0
...

この手法を用いる場合の注意点は、ドメイン分割による並列化(主にCPU runの場合に適用されます)を利用していると、隣のセルを超える長さの相互作用は入れられない点です。実質的に、数十Åぐらいの比較的近距離の長さまでの相互作用だけが可能です。

• [ angles ]のfunctype=8, [ dihedrals ]のfunctype=8で同様に角度、二面角の任意ポテンシャルが追加可能。

• nonbond相互作用についても、coulombtype=User, vdwtype=Userで同様にtable引きが可能です。

Gromacs改造ガイド

→Gromacs改造ガイド(執筆中)で詳述します。