前回は系のエネルギー最小化についての説明を行った。今回は次に行う昇温シミュレーションの解説とエネルギー最小化、昇温シミュレーションのスクリプトの説明を行う。

系のエネルギーを最小化した時点では、系は速度を持っていない。主にシミュレーションの速度(運動量)成分を決定するために昇温シミュレーションは行われる。 最初から目的の温度を与えれば良いでは無いかと思うかもしれないが、目的の温度までに急に加熱してしまうと非常に大きな運動量により時間積分が発散し、シミュレーションが停止してしまう場合がある。そのためある程度の時間ステップを用いて目的温度まで上昇させる。これを昇温(加熱)シミュレーションという。

以下にエネルギー最小化と昇温シミュレーションを行うNAMD configuration fileを示す。

#############################################################
## JOB DESCRIPTION                                         ##
#############################################################

# Minimization and Equilibration of 
# ralgds in a Water Box


#############################################################
## ADJUSTABLE PARAMETERS                                   ##
#############################################################

structure          /home/users/hogehoge/common/bpti_wb_ions.psf
coordinates        /home/users/hogehoge/common/bpti_wb_ions.pdb

set temperature    300
set outputname     bpti_heated

firsttimestep      0


#############################################################
## SIMULATION PARAMETERS                                   ##
#############################################################

# Input
paraTypeCharmm	    on
parameters          /home/users/hogehoge/topper/par_all27_prot_lipid.inp 
temperature         300


# Force-Field Parameters
exclude             scaled1-4
1-4scaling          1.0
cutoff              13.0
switching           on
switchdist          11.0
pairlistdist        20.0

#fixedAtoms on
#fixedAtomsFile /home/users/hogehoge/common/fix_bpti.pdb
#fixedAtomsCol B

#上記は拘束を行いたい場合
#sedやperl、vmd等のいずれかを用いて炭素原子のB値のみを1にしたファイルfix_bpti.pdbを用いることで
#炭素骨格を維持しながら立体障害を取り除くことが出来る。

#constraints         on
#consRef             /home/users/hogehoge/common/bpti_restrain_ca.pdb
#consKFile           /home/users/hogehoge/common/bpti_restrain_ca.pdb
#consKCol            B

#上記は調和振動子的に拘束を加えたい場合
#sedやperl、vmd等のいずれかを用いてアルファ炭素のみのB値を1にする。
#fixed atomsは指定した原子のエネルギーの計算を行わないことで位置の固定を行うが、constraintsは系に調和振動子拘束を加えることで固定をする。
#fixed atomsは本計算や昇温、平衡化シミュレーションで使用した例を見たことはない

# Integrator Parameters
timestep            2.0  ;# 2fs/step
rigidBonds          all  ;# needed for 2fs steps
nonbondedFreq       1
fullElectFrequency  2  
stepspercycle       20


# Constant Temperature Control
langevin            on    ;# do langevin dynamics
langevinDamping     5.0     ;# damping coefficient (gamma) of 1/ps
langevinTemp        $temperature
langevinHydrogen    off    ;# don't couple langevin bath to hydrogens


# Periodic Boundary Conditions
cellBasisVector1    60.0   0.0   0.0
cellBasisVector2     0.0  60.0   0.0
cellBasisVector3     0.0   0.0   60.0
cellOrigin           0.0   0.0   0.0

wrapAll             on

# PME (for full-system periodic electrostatics)
PME                 yes
PMEInterpOrder      6
PMEGridSizeX        64
PMEGridSizeY        64
PMEGridSizeZ        64

# Constant Pressure Control (variable volume)
useGroupPressure      yes ;# needed for rigidBonds
useFlexibleCell       no
useConstantArea       no

langevinPiston        on
langevinPistonTarget  1 ;#  in bar
langevinPistonPeriod  200.0
langevinPistonDecay   100.0
langevinPistonTemp    $temperature


# Output
outputName          $outputname

restartfreq         500     ;# 500steps = every 1ps
dcdfreq             50000
xstFreq             250
outputEnergies      100
outputPressure      100

#############################################################
## EXECUTION SCRIPT                                        ##
#############################################################

reassignFreq       100
reassignTemp       0
reassignIncr       1
reassignHold       300

minimize           2000

#fixedAtoms         off
#原子を固定してエネルギー最小化を実行していた場合は、上記のコメントアウトをはずさないと骨格が固定されたまま昇温シミュレーションが走る

run 100000 #you need more long time step,change me! In the current setting,simulation run 200ps.

#constraintScaling  0
#上記は拘束を外す命令。拘束を解いてから実行したい場合。

#run 15000 ;# 30ps

拘束を加えたい場合はfixedatomsやconstraintsまわりのコメントアウトを外せば良い。 また、この計算は長めの実行時間を要するため、席を外して珈琲を飲みながら待つと良いかもしれない。 長い本計算を行うのには必須であるのだが、上記の計算も待てないという方には、並列化という手が存在する。 NAMDをmpiバージョンでインストールしている前提で話を進めるが以下のコマンドにより実行することで、複数コアまたは複数ノードによる並列化を行うことが出来る。

mpiexec -np 10 namd2 bpti_heating.conf > bpti_heated.log
#もしくは
#charmrun +p10 namd2 bpti_heating.conf > bpti_heated.log 
#mpirun -np 10 namd2 bpti_heating.conf > bpti_heated.log
#-npの後ろや+pの後ろはコア数の指定である。

上記のような処理によって温度一定、圧力一定のシミュレーション用の準備が完了した。