蛋白质分子动力学模拟

来源：“分子动力学”公众号
链接：https://mp.weixin.qq.com/s/drbwbTbpjgE8Y4dCmM5Eig

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