Amber学习第七天：蛋白分子动力学的模拟

1.系统的准备：

假设有一个PDB文件含有500多个残基并且构成两条链。首先修饰PDB文件：除去remarks, connectivity data 和 HETATM lines的数据(不包括晶体水)，然后增加TER标签在蛋白链间和改变一些氨基酸的名字如：HIS to HIE, HID or HIP；和CYS。

tleap -s -f leaprc.ff99SB    (告诉tleap将载入ff99SB力场，接下来载入蛋白）

> mol = loadpdb protein.pdb     (这将增加氢到结构中)

> solvatebox mol TIP3PBOX 12.0  （加12埃的水盒子）

> charge mol                     （检查分子的电荷）

> addions mol Na+ 0           （增加抗衡离子来中性化系统）

> saveamberparm mol protein.prmtop protein.inpcrd    （创建prmtop and inpcrd文件）

> quit

protein.prmtop文件含有分子拓扑性，力场参数，原子和残基名字。rotein.inpcrd文件含有起始的坐标文件，用这两个文件产生PDB文件。

ambpdb -p protein.prmtop < protein.inpcrd > protein_solvated.pdb

2.能量最小化:

将运行两步最小化：第一步固定蛋白最小化水分子位置，第二步最小化整个系统。当运行sander时，应该用控制文件：min1.in 和 min2.in

min1.in：

energy minimization stage 1
 &cntrl
  imin=1, maxcyc=5000, ncyc=2500,
  cut=10.0, ntb=1,
  ntc=1, ntf=1,
  ntpr=10,
  ntr=1,
  restraintmask=':1-500',
  restraint_wt=2.0
 /

Information in the input file:
imin=1  perform minimization
maxcyc=5000  maximum number of minimization cycles
ncyc=2500  method of minimization will be switched from steepest descent to conjugate gradient after ncyc cycles
cut=10.0  specify the nonbonded cutoff, in Angstroms
ntb=1  periodic boundary, constant volume
ntc=1  SHAKE is not performed (for better energy convergence)
ntf=1  force evaluation, complete interaction is calculated
ntpr=10  every ntpr steps energy information will be printed in human-readable form to mdout file
ntr=1  flag for restraining specified atoms using harmonic potential
restraintmask=':1-500'  string that specifies the restrained residues
restraint_wt=2.0  the weight (in kcal/mol-A^2) for the positional restraints

运行：sander -O -i min1.in -p protein.prmtop -c protein.inpcrd -o protein_min1.out -r protein_min1.rst -ref protein.inpcrd

mdin  control data for the min/md run ；prmtop  molecular topology, force field, atom and residue names ；inpcrd  initial coordinates
mdout  user readable state info and diagnostics ；rstrt  final coordinates and velocities ；refc  reference coords for position restraints

min2.in:

energy minimization stage 2
 &cntrl
  imin=1, maxcyc=10000, ncyc=5000,
  cut=10.0, ntb=1,
  ntc=1, ntf=1,
  ntpr=10,
 /

运行：sander -O -i min2.in -p protein.prmtop -c protein_min1.rst -o protein_min2.out -r protein_min2.rst

3.加热

给系统加热从0K到300K并在体积不变下带有位置限制运行50ps动力学。

heat.in

heating
 &cntrl
  imin=0,irest=0,ntx=1,
  nstlim=25000,dt=0.002,
  ntc=2,ntf=2,
  cut=10.0, ntb=1,
  ntpr=500, ntwx=500,
  ntt=3, gamma_ln=2.0,
  tempi=0.0, temp0=300.0,
  ntr=1, restraintmask=':1-500',
  restraint_wt=1.0,
  nmropt=1
 /
 &wt TYPE='TEMP0', istep1=0, istep2=25000,
  value1=0.1, value2=300.0, /
 &wt TYPE='END' /

imin=0  do molecular dynamics ；irest=0  flag to restart the run, no effect ；ntx=1  no initial velocity information ；nstlim=25000  number of MD-steps to be performed
dt=0.002  time step (psec) ；ntc=2  flag for SHAKE, bonds involving hydrogen are constrained ；ntf=2  force evaluation, bond interactions involving H-atoms omitted
ntwx=500  every ntwx steps the coordinates will be written to mdcrd file ；ntt=3  use Langevin thermostat ；gamma_ln=2.0  collision frequency in temperature regulation
tempi=0.0  initial temperature ；temp0=300.0  reference temperature ；nmropt=1  varying conditions ；TYPE='TEMP0'  varies the target temperature
istep1=0, istep2=25000  change is applied over steps istep1 through istep2 ；value1=0.1, value2=300.0  values of the change corresponding to istep1 and istep2, respectively

运行：sander -O -i heat.in -p protein.prmtop -c protein_min2.rst -o protein_heat.out -r protein_heat.rst -x protein_heat.mdcrd -ref protein_min2.rst

4.平衡

为了平衡系统，将运行500ps的分子动力学在常压下没有位置的限制

equil.in：

equilibration
 &cntrl
  imin=0, irest=1, ntx=5,
  nstlim=250000, dt=0.002,
  ntc=2, ntf=2,
  cut=10.0, ntb=2, ntp=1, taup=2.0,
  ntpr=500, ntwx=500, ntwr=5000,
  ntt=3, gamma_ln=2.0,
  temp0=300.0,
 /

Information in the input file:
irest=1, ntx=5  restart calculation, requires velocities in coordinate input file ；ntb=2  periodic boundary, constant pressure ；ntp=1  flag for constant pressure dynamics, md with isotropic position scaling
taup=2.0  pressure relaxation time (in ps) ；ntwr=5000  every ntwr steps during dynamics, the restrt file will be written, ensuring that recovery from a crash will not be so painful

运行：sander -O -i equil.in -p protein.prmtop -c protein_heat.rst -o protein_equil.out -r protein_equil.rst -x protein_equil.mdcrd

5.Production dynamics

在常压下运行10ns动力学。获得的轨迹被用于分析

prod.in:

production dynamics
 &cntrl
  imin=0, irest=1, ntx=5,
  nstlim=5000000, dt=0.002,
  ntc=2, ntf=2,
  cut=10.0, ntb=2, ntp=1, taup=2.0,
  ntpr=1000, ntwx=1000, ntwr=50000,
  ntt=3, gamma_ln=2.0,
  temp0=300.0,
 /

prod.in文件于equil.in的差异主要在于nstlim, ntpr, ntwx, and ntwr值的不同

当运行five million MD-steps时，他是合理的平行的运行：

mpirun -np 32 sander.MPI -O -i prod.in -p protein.prmtop -c protein_equil.rst -o protein_prod.out -r protein_prod.rst -x protein_prod.mdcrd

http://enzyme.fbb.msu.ru/Tutorials/