蛋白质分子动力学模拟

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 andinpcrd文件）

> 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 afterncyc 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 tomdout 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 throughistep2；value1=0.1, value2=300.0  values of the change corresponding toistep1 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