GROMACS.mdpオプション¶
ここではGROMACSでMDシミュレーションを実行する時に必要となるmdpファイルの設定について、その内容を日本語で解説する。 より詳しく知りたい場合は公式ドキュメントを参照してほしい。
一般情報・前処理¶
define¶
プリプロセッサに渡す定義を指定する。基本的には次の2つが指定できる。
define | 説明 |
---|---|
-DFLEXIBLE | 水モデルに剛体モデルではなくフレキシブルモデルを使用する |
-DPOSRES | posre.itpなどに設定された位置拘束を適用する |
define = -DFLEXIBLE
;define = -DPOSRES コメントアウトして必要な時のみ有効にする
実行制御¶
integrator¶
演算アルゴリズムを指定する。MDシミュレーションを実行する時は基本的にmd
を使用する。 一方、エネルギー最小化にはsteep
を使用することが多い。
integrator = md
integrator = steep # 最急降下法アルゴリズム
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]
エネルギー最小化¶
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
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