3. Equilibration

MDシミュレーションでは、システムが定常状態に達するまでに巨大分子とその周囲の溶媒は通常10から数100psをかけて緩和される。初期の非定常状態セグメントの simulated trajectory は平衡化時の性質を計算する際には普通 は廃棄される。MDシミュレーションではこのステージは平衡時ステージと言われる。

コールは現在でも個人の好みの問題である。いくつかのプロトコルは非常に精巧な進め方を採用しステップ状に温度を上げるが、別の方法では直線的に温度と上昇させ想定する温度にシステムがなるまで熱する。 

本説明では、 AMBER マニュアルで示唆されている方法に従い、２段階平衡を使用する。最初に100Kから300Kまで10psかけて熱を上げる。その間は容積一定で平衡化をおこなう。その入力ファイルは以下の通りである。:

Heating up the system equilibration stage 1
 &cntrl
  nstlim=5000, dt=0.002, ntx=1, irest=0, ntpr=500, ntwr=5000, ntwx=5000,

  tempi =100.0, temp0=300.0, ntt=1, tautp=2.0, ig=209858,
 
  ntb=1, ntp=0,

  ntc=2, ntf=2,
 
  nrespa=2,
&end

一行ずつ説明すると

Line 1:
5000ステップを2fs刻みで実行する(nstlim=5000 dt=0.002)。新しく実行するので参照する velocity information (ntx=1, irest=0)を使用せず、 5000 steps (ntwr=5000)ごとにエネルギー出力を書き 5000 step(ntwx=5000) ごとにcoordinateを保存する。

Line 2:
初期温度を 100K (tempi=100.0) 、参照温度 300K (temp0=300.0) で Berendsen coupling アルゴリズムを使って定常温度にしている(ntt=1)。温度カップリングのtime constantは 2 ps (tautp=2.0).です。初期速度 assignment に対して random number seed を (ig=209858) としている。同じseedを用いると全く同じtrajectoryが得られます。もし初期速度の異なるmultiple trajectoryで実行したいなら、別のrandom seed number を使う必要がある。

Line 3:
一定の体積(ntb=1)と圧力制御をしない (ntp=0) 周期境界を用いる。

Line 4:
水素を含む結合のみを拘束するために SHAKE algorithm を 使い (ntc=2)、その結合の評価を無視する (ntf=2) 。 coordinateのリセットにたいして強いトレランスを設定する。

Line 5:
計算を少しだけスピードアップするため、 force field の slow-varying terms を評価により大きなステップ時間を使う (nrespa=2) 。この設定でステップ時間はnrespa dtと等しくなる ( 2 * 2 fs = 4 fs.)。

note: nrespa オプションは AMBER 7 以上で可能で、それより古い AMBER を用いる場合はinput ファイルからnrespaを除く。

この変数は force fieldの遅い slowly-varying termsを評価することを可能にする。 Particle Mesh Ewald (PME)で, "slowly-varying" (now) は逆数和を意味（mean）する。generalized Born (GB)では "slowly-varying" forcesは 効果的半径や "inner" cutoffより大きい距離にあるペアの相互作用に関係したforceである。もし NRESPA >1 であれば、slowly-varying force　は各ステップごとに計算される。forece は適正に調整されそのステップのインパルスを導く。 nrespa x dt <= 4 fs ならエネルギー保存は厳密には調整されない。しかし、 nrespa x dt>4ならばシミュレーションは不安定になる。エネルギーと関連した値はnrespaステップごとに利用できるが他の時間の値は無意味であることを心しておく。------とamberのマニュアルには書いてあるが、現在のところ私には理解不能である。
このファイルを "eq_v.in" として保存する。準備が出来たので実行を開始する。（vは容積一定で平衡化したと言う意味でつけてある。）

$ sander -O -i eq_v.in -p wt1mg.parm7 -c wt1mg_min_all.rst -r wt1mg_eq_v.rst -x wt1mg_eq_v.crd -o wt1mg_eq_v.out

この計算にはデュアルCPU Pentium IIIで２時間以上かかる。それで結果のファイルをダウンロードして使っても良い。 (付記：現在ではpmemd.MPIを用いてcpu 8個を使うと3分、sander (iMac i7) で10分、cuda では12秒で終わる。） download page
----------------------------------------------------------------------------------------------------------------
Conflex社　AMBER チュートリアルでは以下のように記されている。 温度制御は AMBER 8 から新 しく導入されているランジュバン（ Langevin） 法(NTT = 3)を使用し、温度を 300 K に保ちます。この温度制御法では、GAMMA_LN で与えられる衝突頻度によりランジュバン 動力学計算を行います。この温度制御法は、古いバージョンの AMBER で推奨されていた Berendsen 温度カップリング(NTT = 1)よりも、 系の温度を平衡化するのに著しく効果的です。古い Berendsen 法での最 大の問題は、運動エネルギーが単純に目的の温度に適合しているとしてい るアルゴリズムでした。このアルゴリズムでは、分子の全体で温度が超え ていても、何もしません。これにより、溶媒の温度が高く溶質が冷たくなっ ているという現象の原因になったりします。これを避けるため、シミュ レーションの全体を通して分子をゆっくりと暖めるという手の込んだ温度 スケーリング法が推奨されていました。新しい Langevin 法は温度の平衡 化に関してかなり効果的で、AMBER 8 以降では温度平衡化の推奨される 手法になっています。しかし、Langevin 温度制御法は注意深く使用して ください。効率よく系の温度が平衡に達しますので、系の動力学が速くな ります。特に実際に溶媒を置いた動力学計算では、一旦 ntt=3 を使用して 系を平衡化させ、ntt=1 に切り替えるかあるいは純粋なニュートン動力学 計算(ntt = 0)を行った方がいいでしょう。しかしチュートリアルの例では、 シミュレーション全体を ntt=3 にしたままにします。この Langevin 温度 制御法は効果的で、今回のようにかなり良い構造の場合、系をゆっくりと(0 K から室温までおよそ 20 ps かけて)上げる必要はなく、300 K から実際 に始めることができます。ですからこのシミュレーションでは、NTT = 3 と GAMMA_LN = 1 をセットします。そして、初期温度と最終温度を 300K にセットします(TEMPI = 300.0, TEMP0 = 300.0)。これは、系の温度が 300 K 近傍にとどまる事を意味します。またここの2例ではそれぞれ、時 間間隔 1 fs(DT = 0.001)で 100,000 ステップ(NSTLIM = 100000)走らせ、 100 ps(100,000 x 1 fs)の長さのシミュレーションにします。 このことはntt=1（Berendsen coupling アルゴリズム）が入力ファイルに記されていることに関係している。

[Previous][Index][Next]