AMBERを用いたNVT, NPT平衡化

AMBERでNVT, NPT平衡化をする場合の設定方法です。この記事はAMBERを用いた簡単なエネルギー最小化の続きです。

流れ

GROMACSを用いたNVT, NPT平衡化の場合と同じく、9段階でposition restraintsを徐々に外していくやり方です。初めての方はまずそちらの記事を読んできてください。md1ではNVT一定、md2~md9ではNPT一定条件のシミュレーションになっています

Position Restraintをかけるべき重原子の設定はrestraintmaskで定義することができます。これはAMBER MASK記法で書かれており、この文法の詳細はAMBERマニュアルを参照してください。ここでの:1-170 & !@H=は残基番号1-170の水素原子以外という意味です（:マークが残基番号の参照,&は論理記号and, !@H=はHで始まる原子以外、つまり水素原子以外ということになります。全体としては、「残基番号1〜170までの全原子のうち、水素原子以外」という意味になります）。

スクリプト例

md1.in

Heat system (constant volume)
 &cntrl
  imin=0,                         ! Molecular Dynamics
  irest=0,                        ! DO NOT restart MD simulation from a previous run.
  ntx=1,                          ! Coordinates and velocities will not be read
!  irest=1,                       ! Restart MD simulation from a previous run.
!  ntx=5,                         ! Coordinates and velocities will be read from a previous run.
  nstlim=100000,                  ! Number of MD steps ( 200 ps )
  dt=0.002,                       ! Timestep (ps)
  igb=0,                          ! No generalized Born term is used (Default)
  ntp=0,                          ! No pressure scaling (Default)
  ntb=1,                          ! Constant Volume. NVT simulation.
  ntc=2,                          ! SHAKE on for bonds involving hydrogen atoms
  ntf=2,                          ! No force evaluation for bonds with hydrogen
  cut=8.0,                        ! Nonbonded cutoff (Angstroms)
  iwrap=1,                        ! the coordinates written to the restart and trajectory files will be "wrapped" into a primary box.
  ntpr=5000,                      ! Print to mdout every ntpr steps
  ntwx=5000,                      ! Write to trajectory file every ntwc steps
  ntwr=5000,                      ! Every ntwr steps during dynamics, the "restrt" file will be written
  ntt=3,                          ! Langevin thermostat
  gamma_ln=2.0,                   ! Collision frequency for thermostat
  ig=-1,                          ! Random seed for Langevin thermostat
  tempi=100.0,                    ! Initial Temperature. For the initial dynamics run, (NTX < 3) the velocities are assigned from a Maxwellian distribution at TEMPI K.
  temp0=300.0,                    ! Reference temperature at which the system is to be kept.
  ntr=1,                          ! Harmonic position restraints ON. The restrained atoms are determined by the restraintmask string.
  restraintmask=':1-170 & !@H=',  ! String that specifies the restrained atoms when ntr=1.
  restraint_wt=10.0,              ! The weight (in kcal mol-1 Å-2) for the positional restraints.
  nmropt=1,                       ! NMR restraints and weight changes will be read.
  ioutfm=1,                       ! Binary NetCDF trajectory
 /
 &wt TYPE='TEMP0', istep1=0, istep2=100000,
  value1=100.0, value2=300.0, /
 &wt TYPE='END'
 /

restraintmask=':1-170 & !@H='の部分をご自身の系によって適切な原子の範囲に置き換えましたか？例えば1~408がタンパク質＋リガンドで残りが水やイオン原子ならば、ここはrestraintmask=':1-408 & !@H='になっていると思います。

md2.in ~ md9.in

md2からmd9までは、下記のrestraint_wt=xxxの値を10.0, 5.0, 2.0, 1.0, 0.5, 0.2, 0.1,
0という風に徐々に小さな値に設定します。

(constant volume)
 &cntrl
  imin=0,                         ! Molecular Dynamics
!  irest=0,                        ! DO NOT restart MD simulation from a previous run.
!  ntx=1,                          ! Coordinates and velocities will not be read
  irest=1,                        ! Restart MD simulation from a previous run.
  ntx=5,                          ! Coordinates and velocities will be read from a previous run.
  nstlim=50000,                   ! Number of MD steps ( 100 ps )
  dt=0.002,                       ! Timestep (ps)
  igb=0,                          ! No generalized Born term is used (Default)
  ntp=1,                          ! MD simulations with isotropic position scaling)
  ntb=2,                          ! Constant Pressure. NPT simulation.
  ntc=2,                          ! SHAKE on for bonds involving hydrogen atoms
  ntf=2,                          ! No force evaluation for bonds with hydrogen
  cut=8.0,                        ! Nonbonded cutoff (Angstroms)
  iwrap=1,                        ! the coordinates written to the restart and trajectory files will be "wrapped" into a primary box.
  ntpr=5000,                      ! Print to mdout every ntpr steps
  ntwx=5000,                      ! Write to trajectory file every ntwc steps
  ntwr=5000,                      ! Every ntwr steps during dynamics, the "restrt" file will be written
  ntt=3,                          ! Langevin thermostat
  gamma_ln=2.0,                   ! Collision frequency for thermostat
  ig=-1,                          ! Random seed for Langevin thermostat
  temp0=300.0,                    ! Reference temperature at which the system is to be kept.
  ntr=1,                          ! Harmonic position restraints ON. The restrained atoms are determined by the restraintmask string.
  restraintmask=':1-170 & !@H=',  ! String that specifies the restrained atoms when ntr=1.
  restraint_wt=xxx,              ! The weight (in kcal mol-1 Å-2) for the positional restraints.
  ioutfm=1,                       ! Binary NetCDF trajectory
 /

restraint_wt=xxxは正しく全部置き換えましたか？
これらmd1.in ~ md9.inを用意した上で、次のシェルスクリプトrun.shを実行します。

# !/bin/sh

# OpenMPIを使うようにPATHの設定をする。
MPIROOT=/Users/${USER}/apps/openmpi/3.1.4_gcc9.1.0
export PATH=$MPIROOT/bin:$PATH
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$MPIROOT/lib
export MANPATH=$MANPATH:$MPIROOT/share/man
# AMBER 18のPATH設定
source /Users/${USER}/apps/amber18/amber.sh

# main process
# AMBERの最小化を実行し終わったときに生成されるmin2.rst7のPATHを指定する
prev_rst7="../minimize/min2.rst7"

for (( i=1;i<10;i++ ));do
    echo "Start the ${i} cycle."
    #まだmd${i}ステップが終わっていなければ(md${i}.rst7が存在していないなら)実行
    if [ ! -f md${i}.rst7 ];then
        pmemd.cuda.MPI -O \
            -i md${i}.in \
            -o md${i}.out \
            -p ../../top/leap.parm7 \
            -c ${prev_rst7} \
            -ref ${prev_rst7} \
            -r md${i}.rst7 \
            -x md${i}.nc \
            -inf md${i}.info || exit 1
    #すでにmd${i}ステップが終わっていれば次のサイクルへ移動
    else
        echo "md${i}.rst7 exists. Starts next cycle."
    fi
    prev_rst7="md${i}.rst7"
done

# 不要であればmd1~md9.ncを削除
rm md[1-9].nc

GROMACSと大体同じ流れです。GPU＋MPI版のPMEMDをインストールしており、利用可能な場合、pmemd.cuda.MPIを利用した方がとても速くシミュレーションが行われます。GPU Accelerationされていない場合はここのコマンドをmpirun -np ${NCPUS} pmemd.MPIに設定します。なぜかわからないですけれどGPUのときはmpirunを使おうとすると逆に遅くなってしまいました。