AMBERを用いたNVT, NPT平衡化

AMBERでNVT, NPT平衡化をする場合の設定方法です。この記事は[AMBERを用いた簡単なエネルギー最小化](https://qiita.com/Ag_smith/items/842a7946c35b4bd276e8)の続きです。

## 流れ

[GROMACSを用いたNVT, NPT平衡化](https://qiita.com/Ag_smith/items/24fbb9467787f69fbf1c)の場合と同じく、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`を実行します。

``` bash
#!/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`を使おうとすると逆に遅くなってしまいました。