GROMACS.mdpオプション¶

ここではGROMACSでMDシミュレーションを実行する時に必要となるmdpファイルの設定について、その内容を日本語で解説する。より詳しく知りたい場合は公式ドキュメントを参照してほしい。

一般情報・前処理¶

title¶

シミュレーションの名前を定義する。任意の名前を指定できる。

title = OPLS NVT SIMULATION

define¶

プリプロセッサに渡す定義を指定する。基本的には次の2つが指定できる。

define	説明
-DFLEXIBLE	水モデルに剛体モデルではなくフレキシブルモデルを使用する
-DPOSRES	posre.itpなどに設定された位置拘束を適用する

define = -DFLEXIBLE
;define = -DPOSRES    コメントアウトして必要な時のみ有効にする

実行制御¶

integrator¶

演算アルゴリズムを指定する。MDシミュレーションを実行する時は基本的にmdを使用する。一方、エネルギー最小化にはsteepを使用することが多い。

integrator = md
integrator = steep    # 最急降下法アルゴリズム

nsteps¶

シミュレーションのステップ数を指定する。

nsteps = 50000

dt¶

シミュレーションのタイムステップを指定する。**単位は[ps]**であることに注意。一般的な全原子のシミュレーションでは、タイムステップは1[fs]や2[fs]とすることが多い。

dt = 0.002    # 2[fs]

シミュレーション時間（どのくらいの時間のMD計算を行うか）はnstepsとdtの組み合わせによって決定する。次の例は50[ns]のシミュレーションを実行する時の設定である。

# 50[ns]のシミュレーション設定
nsteps = 25000000
dt     = 0.002         # 2[fs] x 25000000[steps]

エネルギー最小化¶

emtol¶

エネルギー最小化において、系の最大力がこの値よりも小さくなった場合に収束したと判定する。デフォルトは10.0[kJ mol-1 nm-1]である。

emtol = 1000.0

emstep¶

エネルギーの初期ステップサイズを指定する。デフォルトは0.01[nm]である。

emstep = 0.01

エネルギー最小化は以下のような設定とともに、非結合相互作用を定義して行うことが一般的である。具体的なmdpファイルの内容はGROMACSのチュートリアルを確認してほしい。なお、初期構造に原子のオーバーラップがあった場合、この段階で系のエネルギーが跳ね上がるので、原因となる原子を取り除くか、初期構造を操作する必要が出てくる。

# 一般的なエネルギー最小化のパラメータ設定
md     = steep
emtol  = 1000.0
emstep = 0.01
nsteps = 50000

出力制御¶

nstxout¶

座標をトラジェクトリファイルに書き込むステップ間隔を指定する。出力ファイルはtrr形式となる。デフォルトでは0（出力しない）設定なので注意すること。しかし、最終構造の座標は常に出力される。

nstxout = 500

nstvout¶

速度をトラジェクトリファイルに書き込むステップ間隔を指定する。 nstxoutと同様にtrr形式であり、デフォルトでは出力しない。最終速度は常に出力される。

nstvout = 500

nstfout¶

力をトラジェクトリファイルに書き込むステップ間隔を指定する。デフォルトでは出力しない。

nstfout = 0

nstlog¶

エネルギーをログファイルに書き込むステップ間隔を指定する。デフォルトは1000 [steps]であり、最終エネルギーは常に出力される。

nstlog = 1000

nstcalcenergy¶

エネルギーを計算するステップ間隔を指定する。次のnstenergyでファイルに出力することができる。デフォルトは100 [steps]であり、0 [steps]を指定することはできない。

nstcalcenergy = 100

nstenergy¶

エネルギーをバイナリのエネルギーファイルに書き込むステップ間隔を指定する。デフォルトは1000 [steps]であり、nstcalcenergyの整数倍の値のみが指定できる。最終エネルギーは常に出力される。

nstenergy = 1000

nstxout-compressed¶

非可逆圧縮をした座標をトラジェクトリファイルに書き込むステップ間隔を指定する。出力ファイルはxtc形式となる。デフォルトでは0（出力しない）となっているが、trr形式よりもファイルサイズが小さいため、座標のみの解析を行う場合は同時に出力しておくことが多い。

nstxout-compressed = 500

結合情報（拘束条件）¶

constraints¶

結合または角度の拘束を指定する。構造制御を行う場合や、結合のダイナミクスが目的のタイムスケールで必要ない場合などに使用する。指定できるものの一部を次に示す。

constraints	説明
none	拘束を行わない結合や角度は全てポテンシャル関数で記述される
h-bonds	水素原子との結合に拘束条件を適用する
all-bonds	全ての結合に拘束条件を適用する
h-angles	全ての結合に加えて、水素原子を含む角度を結合拘束として適用する
all-angles	全ての結合と角度に拘束条件を適用する

基本的なシミュレーションの場合はconstraints = h-bondsを指定することが多い印象である。

constraint_algorithm¶

結合拘束を行う時に使用するアルゴリズムを指定する。以下の2つがある。

constraint_algorithm	説明
LINCS	LINCSアルゴリズム
SHAKE	SHAKEアルゴリズム

LINCSアルゴリズムがデフォルトとなっており、こちらの方が高速であるが、角度の拘束条件を扱う場合はSHAKEアルゴリズムを使用する必要がある。

continuation¶

拘束条件を新しく適用するか、適用せずに以前のものを引き継ぐかを指定できる。新しいシミュレーションの場合にはcontinuation = no、以前からの継続または再計算（rerun）の場合にはcontinuation = yesを指定すればよい。

以上3つが必須である。次の例は、NVT平衡化が終わった後にNPT平衡化をしたい場合の設定例である。

# 結合拘束条件
constraints           = h-bonds
constraint_algorithm  = LINCS
continuation          = yes

近接原子情報¶

cutoff-scheme¶

近接原子を探索し、相互作用のペアリストを作成するための方法を指定する。 Verletとgroupがあるが、基本的に前者を使用する印象。

nstlist¶

ペアリストの更新頻度を指定する。デフォルトは10[steps]となっている。 0[steps]にするとリストは一度だけ作成され、その後は更新されない。エネルギー最小化ではデフォルトでよいが、本計算では20[steps]や40[steps]にするとパフォーマンスがよいらしい。

ns_type¶

近接原子の探索タイプを指定する。 gridとsimpleが選択できる。

ns_type	説明
grid	ボックス内にグリッドを作成し、近接グリッドの原子をチェックする。大きな系では`simple`よりも速い。
simple	単純にペアリストが生成した時、ボックス内の全ての原子をチェックする。

rlist¶

近接原子ペアリストのカットオフ距離を指定する。デフォルトは1[nm]である。 cutoff-scheme = Verletの場合は無視されるらしい。

ペアリストに関してはnstlistの数値だけ調整して、次の3つを設定すればよいだろう。

# 近接原子情報
cutoff-scheme = Verlet
nstlist       = 20
ns_type       = grid

非結合情報（LJ相互作用、クーロン相互作用）¶

coulombtype¶

クーロン相互作用の計算タイプを指定する。特にこだわりがなければPMEを指定しておこう。

coulombtype	説明
Cut-off	`rlist`と`rcoulomb`を使用した通常のカットオフ。 "rlist >= rcoulomb"である必要がある。
Ewald	古典的なエワルドの方法。実空間のカットオフ`rcoulomb`は`rlist`と同じにする必要がある。
PME	高速でスムーズなパーティクルメッシュエワルド法。実空間はエワルドに似ているが、逆空間はフーリエ変換を使用する。

rcoulomb¶

クーロン相互作用のカットオフ距離を指定する。デフォルトは1[nm]。

vdwtype¶

LJ相互作用の計算タイプを指定する。基本的にCut-offを指定する。

vdwtype	説明
Cut-off	`rlist`と`rvdw`を使用したツインレンジカットオフ。"rvdw >= rlist"である必要がある。
PME	LJ相互作用のための高速スムーズパーティクルメッシュエワルド法。

vdw-modifier¶

LJポテンシャルを修正する。カットオフで力をスムーズに0とするためのForce-switchやポテンシャルを0とするためのPotential-switchが利用できる。これらのオプションはrvdw-switchからrvdwまでの範囲に適用される。

rvdw-switch¶

LJポテンシャルスイッチングの開始点（開始距離）。デフォルトは0[nm]で、スイッチングを使用する場合のみ指定する。

rvdw¶

LJ相互作用のカットオフ距離を指定する。デフォルトは1[nm]。

DispCorr¶

長距離分散補正の指定。エネルギーと圧力に使用するEnerPresとエネルギーのみに使用するEnerがある。

クーロン相互作用とLJ相互作用の例を示す。LJはスイッチングを使用することが多い。

# クーロン相互作用
coulombtype  = PME
rcoulomb     = 1.2
# LJ相互作用
vdwtype      = Cut-off
vdw-modifier = Force-switch
rvdw-switch  = 1.0
rvdw         = 1.2

温度制御・圧力制御¶

tcoupl¶

温度制御の方法を指定する。参照温度はref_t、カップリング時定数はtau_tで指定する。 tc_grpsを指定すると、選択したグループごとに別々に温度制御を適用することができる。下記はtcouplに指定できる温度制御法の一部である。

tcoupl	説明
no	温度制御を行わない
Berendsen	Berendsenサーモスタットによる温度制御
Nose-Hoover	Nose-Hooverサーモスタットによる温度制御
V-rescale	速度スケーリングを用いた温度制御

# Nose-Hooverサーモスタットによる温度制御
tcoupl   = Nose-Hoover
tc_grps  = Other    Water
tau_t    = 1.0      1.0
ref_t    = 310      310

pcoupl¶

圧力制御の方法を指定する。参照圧力はref_p、カップリング時定数はtau_pで指定する。また、compressibilityも同時に指定する。一般には、300[K]における水の値である4.5x10-5[bar-1]が用いられる。 pcouplで指定できる圧力制御法の一部を以下に示す。

pcoupl	説明
no	圧力制御を行わない体積一定）
Berendsen	Berendsenバロスタットによる圧力制御
Parrinello-Rahman	Parrinello-Rahmanバロスタットによる圧力制御

ここが重要であるが、圧力制御を行うにはpcoupltypeで制御タイプを指定する必要がある。主に使用するのは次の2つである。

pcoupltype	説明
isotropic	等方性の制御タイプ
semiisotropic	半等方性の制御タイプ xy方向とz方向で独立している

isotropicを使用する場合はref_pとcompressibilityのそれぞれに1つずつの値を設定する。一方でsemiisotropicを使用する場合は、それぞれに2つずつの値を設定する必要がある。

# 等方性圧力制御
pcoupl          = Parrinello-Rahman
pcoupltype      = isotropic
tau_p           = 5.0
ref_p           = 1.0
compressibility = 4.5e-5

# 半等方性圧力制御
pcoupl          = Parrinello-Rahman
pcoupltype      = semiisotropic
tau_p           = 5.0
ref_p           = 1.0       1.0
compressibility = 4.5e-5    4.5e-5