分子動力学ソフトウェアパッケージAMBERを用いて系のエネルギー最小化を行うスクリプト例の紹介です。

GROMACSを用いた簡単なエネルギー最小化というのは以前に紹介していたのですが、本来AMBERが想定しているような、AMBERだけでMDシミュレーションを行うときのスクリプト例を置いておきたいと思います。

環境

macOS, Linux OS. Windows 10の場合はWindows Subsystem for Linuxをインストールしておいてください。
AMBER22とAmberTools22がインストールされている
AMBERのtleapを使ってMDシミュレーションを始めるのに必要なトポロジーファイル(leap.parm7)と初期座標ファイル(leap.rst7)をすでに作成してある

ではここからエネルギー最小化をやっていきます。

エネルギー最小化

エネルギー最小化その1 （min1.in）

1段階目は、水・イオンだけを動かすエネルギー最小化を200ステップくらいやります。

min1.in

moleculr dynamics minimization run 1
 &cntrl
    imin=1,            ! Do Minimization
    nmropt=0,          ! No restraints
    ntx=1,             ! Coordinates, but no velocities, will be read (default)
    irest=0,           ! Do NOT restart the simulation
    ntxo=1,            ! ASCII Output Format of the final coordinates
    ntpr=20,           ! Every ntpr steps, energy information will be printed.
    ntwr=2000,         ! Every ntwr steps during dynamics, the “restrt” file will be written
    ntwx=0,            ! Every ntwx steps, the coordinates will be written to the mdcrd file. if 0, no output.
    ioutfm=0,          ! ASCII format of coordinate and velocity trajectory files (mdcrd, mdvel and inptraj).
    ibelly=1,          ! If ibelly=1, the coordinates except the bellymasked atoms will be frozen.
    bellymask=":WAT,Na+,Cl-",   ! mask for ibelly = 1
    ntr=0,             ! Flag for restraining specified atoms in Cartesian space using a harmonic potential. No restraints.
    ntmin=1,           ! Method of minimization. For NCYC cycles the steepest descent method is used then conjugate gradient is switched on
    maxcyc=200,        ! The maximum number of cycles of minimization. Default = 1.
    ncyc=100,          ! If "ntmin" is 1, the method of minimization will be switched from steepest descent to conjugate gradient after NCYC cycles. Default 10.
    nstlim=2000,       ! number of MD-steps to be performed.
    nscm=0,            ! Flag for the removal of translational and rotational center-of-mass
    jfastw=0,watnam='TIP3',OWTNM='OH2',hwtnm1='H1',hwtnm2='H2',
 /

最後に1行以上の空行をはさまないとあとのsanderを使った実行時にエラーになります。ibelly=1, bellymask=":WAT,Na+,Cl-"あたりが重要な設定。これにより、水・イオン原子以外を完全に固定した状態で他の原子を最小化することができます。ちなみにGROMACSではこれが実装されていないので、10 kcal/molの位置拘束をかけて1段階目の最小化をやっています。GROMACSと違い、固定した対象の原子は本当に座標が変動しません。
この条件で最大200ステップのminimizationを行います（maxcycで指定しています）。

エネルギー最小化その2（min2.in）

min2.in

moleculr dynamics minimization run 2
 &cntrl
    imin=1,            ! Do Minimization
    nmropt=0,          ! No restraints
    ntx=1,             ! Coordinates, but no velocities, will be read (default)
    irest=0,           ! Do NOT restart the simulation
    ntxo=1,            ! ASCII Output Format of the final coordinates
    ntpr=20,           ! Every ntpr steps, energy information will be printed.
    ntwr=2000,         ! Every ntwr steps during dynamics, the “restrt” file will be written
    ntwx=0,            ! Every ntwx steps, the coordinates will be written to the mdcrd file. if 0, no output.
    ioutfm=0,          ! ASCII format of coordinate and velocity trajectory files (mdcrd, mdvel and inptraj).
    ntr=0,             ! Flag for restraining specified atoms in Cartesian space using a harmonic potential. No restraints.
    ntmin=1,           ! Method of minimization. For NCYC cycles the steepest descent method is used then conjugate gradient is switched on
    maxcyc=200,        ! The maximum number of cycles of minimization. Default = 1.
    ncyc=100,          ! If "ntmin" is 1, the method of minimization will be switched from steepest descent to conjugate gradient after NCYC cycles. Default 10.
    nstlim=2000,       ! number of MD-steps to be performed.
    nscm=0,            ! Flag for the removal of translational and rotational center-of-mass
    jfastw=0,watnam='TIP3',OWTNM='OH2',hwtnm1='H1',hwtnm2='H2',
 /

このファイルも最後に1行以上の空行をはさまないとあとのsanderを使った実行時にエラーになります。今度は水・イオンも入れて200ステップのエネルギー最小化を行います。

最小化の実行

無償でダウンロード可能なAmberToolsしか持っていない場合はsanderまたはsander.MPIしか利用できませんが、有償のAMBERをがある場合はpmemdやpmemd.MPIが利用できます。pmemd, pmemd.MPIの方が非常に速いです。

最小化についてはGPU Accelerationプログラムの対象外なので、pmemd.cudaが利用できないことに注意してください。

sanderの実行スクリプト1 （1コアのみ）

sanderしかない場合はこちら

run.sh

#!/bin/sh

sander -O \
  -i min1.in \
  -o min1.out \
  -p /path/to/leap.parm7 \
  -c /path/to/leap.rst7 \
  -r min1.restrt \
  -inf min1.info
sander -O \
    -i min2.in \
    -o min2.out \
    -p /path/to/leap.parm7 \
    -c min1.rst7 \
    -r min2.rst7 \
    -inf min2.info

/path/to/leap.parm7と/path/to/leap.rst7の部分にはそれぞれtleapで作成したAMBERのトポロジーファイルと座標ファイルを入れます。

sander.MPI/pmemd.MPI の実行スクリプト1 （マルチコア）

DO_PARALLEL="mpirun -np ${NCPUS}"でコア数並列を実現することが可能。

run.sh

#!/bin/sh

#コア数の指定
NCPUS=16

DO_PARALLEL="mpirun -np ${NCPUS}"

${DO_PARALLEL} pmemd.MPI -O \
  -i min1.in \
  -o min1.out \
  -p /path/to/leap.parm7 \
  -c /path/to/leap.rst7 \
  -r min1.rst7 \
  -inf min1.info
${DO_PARALLEL} pmemd.MPI -O \
    -i min2.in \
    -o min2.out \
    -p /path/to/leap.parm7 \
    -c min1.rst7 \
    -r min2.rst7 \
    -inf min2.info

小ネタ

imin=1（Read in a trajectory for analysis.）とmaxcyc=1を同時に指定すると、1ステップめのアウトプットファイルにはインプットの座標（トラジェクトリ）データからのエネルギー値を返してくれる。例えば量子計算のエネルギー値の結果と比較して、新しい化合物の力場のパラメータを自身で作成するときの指標として使える。