5. 分子動力学計算¶
分子動力学計算は以下の手順で行います。
Amberがインストールされているディレクトリはロードしたバージョンによって異なりますが、環境変数としてはすべて、$AMBERHOMEに設定されます。
このディレクトリ下の主なディレクトリを以下に示します(例ではamber/16up10)。
$AMBERHOME/内の主なディレクトリ
| ディレクトリ | 説明 |
|---|---|
| dat/ | Amberのデータベース |
| examples/ | デモンストレーションデータ |
| doc/ | マニュアルのファイル |
| bin/ | Amber実行モジュール |
| test/ | 各モジュールを実行するためのスクリプト例 |
Minimization、equilibration、Molecular Dynamicsは、binディレクトリ内にある sander というモジュールで行います。sander の使用方法は下のようになっています。
usage: sander [-O|A] -i mdin -o mdout -p prmtop -c inpcrd -r restrt
[-ref refc -x mdcrd -v mdvel -e mden -frc mdfrc -idip inpdip -rdip rstdip -mdip mddip
-inf mdinfo -radii radii -y inptraj -amd amd.log -scaledMD scaledMD.log] -cph-data <file>
それぞれのファイルは以下の通りです。
sanderによる計算に必要なファイル
| ファイル | 説明 |
|---|---|
| mdin | min/mdにおけるコントロールデータ(input)。 |
| mdout | アウトプットファイル(output)。 |
| prmtop | 分子のトポロジー、パラメータの情報のデータ(input)。 |
| inpcrd | 初期構造と(オプションで)初期速度のデータ(input)。 |
| restrt | 最終の座標、速度、ボックスサイズ(for Constant Pressure)のデータ(output)。 |
| refc | ポジションを束縛するときのデータ(input)。 |
| mdcrd | 座標値のトラジェクトリファイル(output)。 |
| mdvel | 速度のトラジェクトリファイル(output)。 |
| mden | エネルギーのトラジェクトリファイル(output)。 |
| mdinfo | mdoutの一番最新の部分(output)。 \ |
インプットファイルのオプション指定は&cntrlから/の間に記載し、オプションが複数ある場合は「,」で区切り、&cntrl から/までの行は必ず行頭に空白を入れて下さい。
インプットファイルで使用できる全オプションの詳細説明はマニュアルを参照して下さい。
インプットファイルのオプション
| オプション | 説明 |
|---|---|
| imin | 0の場合はMD計算、1の場合はエネルギー極小化計算、5の場合は振動解析を行います。 |
| nmropt | 0の場合はNMRを考慮しない計算、1,2の場合はNMR計算を行います。 |
| ntx | 初期座標、速度および基本セルの読み込み方法の指定フラグです。1がデフォルトとなっており、初期座標のみが呼び出されます。MD計算のリスタートでは5を利用して、初期座標、速度を呼び出します。基本セルはntb>0のときに行われます。なお、irestが1のときのみ、速度が呼び出されます。 |
| irest | リスタートフラグです。0の場合はリスタートを行わず、1の場合は座標や速度をリスタートファイルから読み込みます。 |
| ntpr | mdoutおよびmdinfoを更新するステップ数、デフォルトは50 |
| ntb | 周期的境界条件の設定 |
| cut | カットオフ値、単位はオングストローム |
| maxcyc | 極小化計算の最大サイクル数 |
| ntr | 0より大きい場合に位置束縛条件の設定を行います。デフォルトは0 |
| ntc | SHAKE法による分子固定設定、ntfと合わせて利用します。 |
| ntf | 力の評価設定、ntcと合わせて利用します。 |
| tempi | 初期温度設定 |
| temp0 | 定常温度の設定 |
| gamma_ln | 衝突頻度因子 |
| ntwx | mdcrd |
| ntwr | リスタートファイルを更新するステップ数、デフォルトはnstlim |
| nstlim | MD計算のステップ数 |
| pres0 | 初期圧力設定、デフォルトは1.0 |
| taup | 圧力緩和時間、デフォルトは1.0 |
| ntt | 定温シミュレーションの条件 |
| dt | トラジェクトリー計算における時間ステップ |
5.1. Minimization(構造最適化)¶
5.1.1. Minimizationのインプットファイル¶
#mdin
#minimization (minimizationの実行)
&cntrl
imin = 1, maxcyc = 200,ntpr = 10,
/
以下の実行項目を実施する場合はmdinという名称で上記サンプルを保存してください。
5.1.2. Minimizationの実行(インタラクティブ実行)¶
qrshコマンドで計算ノードにログイン後に以下のコマンドを実行してください。
*サンプルでは最小化サイクルの上限までまわるため、最小化計算は完了しません。
参考エネルギー極小化するためには14223回のステップが必要となります。
$ cp $AMBERSAMPLE/dna_pol/* ./
$ ls
inpcrd inpdip mdin prmtop
$ sander -i mdin -o mdout -p prmtop -c inpcrd
$ ls
inpcrd inpdip mdin mdinfo mdout prmtop restrt
並列計算を行う場合は下記コマンドを実行ください。
なお、サンプルは28並列計算時に5分以内で終了します。
$ mpirun -np [並列数] sander.MPI -O -i mdin -o mdoutp -p prmtop -c inpcrd
計算の詳細を確認する場合はmdoutファイル、概要を確認する場合はmdinfoを確認ください。
5.2. 溶媒を考慮したMM計算(構造最適化)¶
以下のサンプルはAmberのチュートリアルを元にTSUBAME3環境にあわせて各計算は10分以内に終了するように作成しております。
溶媒を考慮する場合は様々な手法がありますが、本書では、溶質を固定した計算(Step1)、全体計算(Step2)を行います。
計算時間は各手法ともに並列化なしで5分程度となります。その他手法については各自論文をご確認ください。
qrshコマンドで計算ノードにログイン後に以下のコマンドを実行してください。(以下のサンプルではAmber10Update10を利用します)
*ログインノードでの計算実行は行わないでください。
$ module load amber/16up10
サンプルファイルをダウンロードし、スクリプトを実行します。
$ cp $AMBERSAMPLE/dna_wat/step1/* ./
$ sh step1.sh
$ cp $AMBERSAMPLE/dna_wat/step2/* ./
$ sh step2.sh
step2.shを実行すると2つのPDBファイルが生成されます。
生成されたPDBを端末にダウンロードすることでBIOVIA DiscoveryStudio等の可視化ツールにて確認ができます。
下図はBIOVIA DiscoveryStudioで可視化した溶質固定後の計算結果、固定化を外した全体計算の結果、両者の重ね合わせを示します。
**溶質を固定した計算結果**
**固定化を外した全体計算の結果**
**溶質固定後の計算結果(赤)、固定化を外した全体計算(青)の結果の重ね合わせ(DNAのみ)**
5.3. 溶媒を考慮したMD計算(CPU/GPU)¶
5.3.1. ヒーティングおよび平衡化¶
MD計算時にも溶質の拘束計算後(Step3)に拘束を解いて全体計算(Step4)を行います。
また、Step3の20psでヒーティングを行い、Step4の100psで平衡化を行います。
Step3が並列化なしで6分、Step4は2ノード、CPU56並列でsander.MPIでは8分以上かかりますが、並列化効率の高いpmemd.MPIでは3分、cudaによる2ノード8GPU並列(pmemd.cuda.MPI)では2分以内に計算が終了します。
Step4ではqsubコマンドをご利用ください。
なお、下記のサンプル計算を行う場合はstep1,2の計算結果が必要となります。
pmemdは真空条件での計算はできませんので、ご注意ください。
非並列オプション実行時
$ cp $AMBERSAMPLE/dna_wat/step3/* ./
$ sh step3.sh
並列オプション実行時
$ cp $AMBERSAMPLE/dna_wat/step3/* ./
$ sh step3_parallel.sh
STEP4の並列計算(2ノード/CPU56並列利用)
$ cp $AMBERSAMPLE/dna_wat/step4/* ./
$ qsub step4_bacth.sh
STEP4の並列計算(CPU2ノード/8GPU並列利用)
$ cp $AMBERSAMPLE/dna_wat/step4/* ./
$ qsub step4_GPUbacth.sh
5.3.2. 結果ファイルの確認¶
Molecular Dynamics計算の出力ファイルmdoutファイルの内容を見て、計算結果を確認して下さい。 mdoutの主な内容は、出力された系の温度や圧力、エネルギー、体積などの情報です。 ここでは例としてCPUで計算した場合のstep4の計算結果の一部を表示します。
NSTEP = 40000 TIME(PS) = 80.000 TEMP(K) = 296.73 PRESS = 96.5
Etot = -28367.6754 EKtot = 5896.0101 EPtot = -34263.6855
BOND = 148.7758 ANGLE = 359.7502 DIHED = 458.8581
1-4 NB = 167.4830 1-4 EEL = -380.6307 VDWAALS = 3886.5010
EELEC = -38904.4229 EHBOND = 0.0000 RESTRAINT = 0.0000
EKCMT = 2718.6134 VIRIAL = 2514.8148 VOLUME = 97839.7312
Density = 1.0402
Ewald error estimate: 0.2648E-03
------------------------------------------------------------------------------
A V E R A G E S O V E R 40000 S T E P S
NSTEP = 40000 TIME(PS) = 80.000 TEMP(K) = 297.10 PRESS = -104.6
Etot = -28396.2451 EKtot = 5903.2761 EPtot = -34299.5212
BOND = 154.1835 ANGLE = 341.9533 DIHED = 458.8789
1-4 NB = 168.4548 1-4 EEL = -360.6335 VDWAALS = 3954.1168
EELEC = -39016.4751 EHBOND = 0.0000 RESTRAINT = 0.0000
EKCMT = 2722.3002 VIRIAL = 2965.6486 VOLUME = 98706.1131
Density = 1.0319
Ewald error estimate: 0.1054E-03
------------------------------------------------------------------------------
R M S F L U C T U A T I O N S
NSTEP = 40000 TIME(PS) = 80.000 TEMP(K) = 19.45 PRESS = 402.2
Etot = 724.1996 EKtot = 386.4596 EPtot = 352.5658
BOND = 11.1978 ANGLE = 15.7098 DIHED = 10.5140
1-4 NB = 4.9522 1-4 EEL = 13.0426 VDWAALS = 201.2018
EELEC = 526.0166 EHBOND = 0.0000 RESTRAINT = 0.0000
EKCMT = 176.9888 VIRIAL = 816.2426 VOLUME = 2935.8983
Density = 0.0276
Ewald error estimate: 0.8083E-04
主なパラメータの内容を以下に示します
mdoutのパラメータ
| パラメータ | 説明 |
|---|---|
| NSTEPS | ステップ数 |
| TIME(PS) | MDシミュレーション時間 |
| TEMP(K) | 温度 |
| PRESS | 圧力 |
| Etot | トータルエネルギー(単位はkcal/mol)( =EKtot + Eptot ) |
| EKto | 運動エネルギー |
| EPtot | ポテンシャルエネルギー( = BOND + ANGLE + DIHED + 1-4NB + 1-4EEL + VDWAALS + ELEC + EHBOND) |
| BOND | 結合の弾性エネルギー |
| ANGLE | 角度の弾性エネルギー |
| DIHED | 二面角によるエネルギー |
| 1-4 NB | 1-4非結合相互作用エネルギー |
| 1-4 EEL | 1-4静電エネルギー |
| VDWAALS | 非結合相互作用エネルギー |
| EELEC | 静電エネルギー |
| EHBOND | 水素結合のエネルギー |
| VOLUME | 体積(単位はA3) |
| Density | 密度(単位はg/cm3) |
5.3.3. mdout_analyzer.py¶
5章で計算を行ったサンプルデータを利用してmdout_analyzer.pyの利用方法をご説明します。
サンプルデータを読み込み、起動する場合は下記コマンドを実行ください。
$ mdout_analyzer.py dna_water_min.out
自分が計算したデータを解析したい場合は下記コマンドを実行ください。
$ mdout_analyzer.py <解析したいmdoutファイル>
コマンドを実行すると下記の画面が立ち上がります。
**mdout_analyzer起動画面**
グラフを作成するにはMdoutAnalyzerのSelect Data to Analyzeよりグラフ化したいボタンを選択し、GraphThem!をクリックしてください。 下図は5.1で計算した構造最適化におけるエネルギーの極小化グラフです。
**Step毎のエネルギーの推移**