以下是分子动力学模拟的一般性步骤, 具体的步骤和过程依赖于你研究的体系和所用的软件, 但这并不影响我们把它当作一个入门指南.
1. 评估体系 首先需要对我们要进行模拟的体系做一个简单的评估, 有三个问题是我们必须要明确的:
2. 选择工具 选择合适的模拟工具, 大前提是它能够实现你所感兴趣的目标. 这需要你非常广泛谨慎地查阅文献, 看看别人用这些工具都做了些什么, 有没有和你的研究体系类似的, 相关的研究. 千万不要做到一半才发现原来你用的工具根本就不能实现你所感兴趣的idea, 切记!
在选择合适的模拟工具时, 主要考虑下面两点:
- 软件的选择. 通常与软件主流使用的力场有关, 软件本身也具有一定的偏向性.
- 蛋白体系: GROMACS, AMBER, NAMD均可
- DNA, RNA体系: 首选AMBER
- 界面体系: DL_POLY比较强大
- 材料体系:LAMMPS是不错的选择
- 力场的选择. 力场用来描述体系中最小单元间的相互作用, 是对实验性质或量子化学计算结果拟合后生成的经验式. 有人会嫌它粗糙, 但它确确实实为我们模拟大系统提供了可能, 只能说关注的切入点不同罢了. 常见的有三类力场: 全原子力场, 联合力场, 粗粒化力场. 此外还有所谓第一代, 第二代, 第三代力场的说法, 这里就不一一列举了.
再次提醒注意: 必须选择适合于我们所关注体系和所感兴趣的性质及现象的力场.
3. 初始结构 通过实验数据或者某些工具得到体系内的每一个分子的初始结构坐标文件. 之后, 我们需要按我们的想法把这些分子按照一定的规则或是随机的放在一起, 从而得到整个体系的初始结构, 这也是我们模拟的输入文件.
4. 输入参数 得到了结构输入文件, 我们还需要力场参数输入文件, 也就是针对我们体系的力场文件. 这通常由所选用的力场决定, 包括电荷, 键合参数和非键参数等势能函数的输入参数.
5. 确定盒子 体系的大小通常由你所选用的盒子大小决定. 我们必须对可行性与合理性做出评估, 从而确定体系的大小, 这依赖于具体的体系.
6. 能量最小化 由于初始构象可能会存在两个原子靠得太近的情况(称之为bad contact), 所以需要在正式模拟开始的第一步对体系进行能量最小化. 比较常用的能量最小化方法有两种, 最速下降法和共轭梯度法. 最速下降法是快速移除体系内应力的好方法, 但是接近能量极小点时收敛比较慢, 而共轭梯度法在能量极小点附近收敛效率高一些. 所以一般做能量最小化时都是先利用最速下降法进行优化, 完成之后再对得到的构象利用共轭梯度法优化一次, 这样做能有效地保证后续模拟的进行.
7. 平衡模拟 你需要设置适当的模拟参数, 并且保证这些参数的设置与力场的构造过程相一致. 举个简单的例子, GROMOS力场是用范德华势双截断来定义范德华参数的, 如果你用GROMOS力场的话也应该用双截断来处理范德华相互作用. 常见的模拟思路是, 先在NVT下限制住你的溶质(剂)做限制性模拟, 这是一个升温的过程, 当温度达到你设定的温度后, 接着做NPT模拟, 此过程将调整体系的压强进而使体系密度收敛.
如何判断体系达到平衡是比较技术性的问题. 简单地讲可以通过以下几种方式:
- 看能量(势能, 动能和总能)是否收敛
- 看体系的压强, 密度等等是否收敛
- 看体系的RMSD是否达到你能接受的范围
- 其他经验
8. 成品模拟 经过一段时间的平衡模拟, 在确定体系已经完全弛豫之后, 就可以开始采集数据了. 运行足够长时间的模拟以确定我们所感兴趣的现象或是性质能够被观测到, 并且务必确保此现象出现的可重复性.
9. 数据分析 数据拿到手后, 很容易通过一些可视化软件得到轨迹动画, 但这并不能拿来发文章. 真正的工作才刚刚开始——分析数据. 你所感兴趣的现象或性质只是表面, 隐含在它们之中的机理才是文章的主题.
小分子过程：
reduce test.pdb > test_h.pdb
antechamber -i test_h.pdb -fi pdb -o test.mol2 -fo mol2 -c bcc -s 2
parmchk -i test.mol2 -f mol2 -o test.frcmod
#use tleap
tleap -f oldff/leaprc.ff99SB
source leaprc.gaff
AAA = loadmol2 test.mol2
list
loadamberparams test.frcmod
check AAA
saveoff AAA test.lib
saveamberparm AAA aaa.prmtop aaa.inpcrd
一个简单模拟过程：
Ssh n1~n21 #进入节点
>tleap #xleap有图形页面，这个没有，服务器上只能运行tleap
>source leaprc.ff14SB ##加载力场（官网上的那个protein应该去掉）
>foo=sequence { ACE ALA NME } ##存储在foo单元中一个丙氨酸二肽
>desc foo ##可以查看该二肽
＞source leaprc.tip3p ##加载水盒
(setbox com "vdw" #按照原子位置设置盒子边界
set com box {30 30 30} #自定义盒子边界)
＞solvatebox foo TIP3PBOX 10.0 ##使用TIP3PBOX水模型，水盒表面距离蛋白10埃
>saveamberparm foo prmtop inpcrd ##保存参数拓扑文件和坐标文件
>quit ##退出tleap
PS：如果体系电荷不为0，需要添加抗衡离子，顺序最好是先溶剂化后中和电荷。
中和电荷的时候有addions和addions2两个算法。addions是在溶质周围添加离子，如果与水分子重合，则删除该位置的水分子；addions2是把溶剂也看做溶质去处理，其它与addions一样。
个人理解：生物分子模拟一般情况下用addions就行。
solvatebox com TIP3PBOX 0
charge com
addions com Na+ 0
addions com Cl- 0
一、能量最小化的脚本 (01_Min.in)
Minimize
&cntrl #起始符
imin=1, #进行最小化
ntx=1, #从inpcrd文件中读取坐标（而不是速度）
irest=0, #不重新启动模拟（不适用于最小化）
maxcyc=2000, #能量最小化的最大周期
ncyc=1000, # 0-ncyc 循环用最速下降算法，然后在 ncyc-maxcyc循环#转换为共轭梯度算法
ntpr=100, #每ntpc个循环打印到amber输出文件
ntwx=0, #没有写入amber轨迹文件（其实也可以写入，=1）
cut=10.0, #非键截止距离，应该小于整个盒子边长的一半；如果大于盒子边长的一半，计算能量的时候，待考察的原子会和周围的镜像相互作用。amber中一般设置为10，namd中设置为12。因为范德华相互作用能在8埃的时候已经基本为0
/
二、升温的脚本 (02_Heat.in)
Heat
&cntrl
imin=0, ##不进行最小化
ntx=1, #无初始速度信息
irest=0, #不重启run，run as a new MD
nstlim=10000, #跑MD的步数
dt=0.002, #时间步长（ps）
ntf=2, #不计算震动约束键的力
ntc=2, #允许震动来约束所有和氢键相连的键
tempi=0.0, #初始温度
temp0=300.0, #最终温度
ntpr=100, #每100步打印一次输出
ntwx=100, #每ntwx步写一次amber轨迹文件
ntwr=100, #每100步写入一次rst文
cut=8.0, ntb=1, #定容的周期边界 ntp=0, #没有压力控制 ntt=3, #Langevin恒温器温度控制
gamma_ln=2.0, #Langevin恒温器碰撞频率
nmropt=1, #核磁共振限制和读取重量改变
ig=-1, #为伪随机数生成程序随机化种子 ！！！！！如果要跑多副本动力学，这个参数必须设置！不然每次都一样！！！！
ntxo=1, #rst文件中的坐标，速度，盒子尺寸信息的格式，=1（ASC），=2（NetCDF）
ioutfm=1, #dcd文件中的坐标和速度的格式，=0（ASC），=1（NetCDF）
ntwv=0, #控制dcd文件中的速度信息是否写入，=0（不写入），=1（每隔ntwx步写入一次速度）
iwrap=1, #处理镜像
/
&wt type='TEMP0', istep1=0, istep2=9000, value1=0.0, value2=300.0 /
&wt type='TEMP0', istep1=9001, istep2=10000, value1=300.0, value2=300.0 /
&wt type='END' /
##最后三行允许恒温器在整个模拟过程中改变目标温度。对于前 #9000步，温度将从0K增加到300K。对于步骤9001到10000，温度将保持在300K。
##因为这个模拟非常简短，所以NTPR和NTWX的设置非常低。使用这些设置进行长时间的MD模拟将产生非常大的输出和轨迹文件，比普通的MD要慢。对于真正的跑MD，你需要增加NTPR和NTWX，不然时间太长。
三、跑动力学的脚本 （03_Prod.in）
Production
&cntrl
imin=0,
ntx=5, #从未格式化的inpcrd坐标文件中读取坐标和速度
irest=1, #重启之前的MD（意味着inpcrd文件提供初始原子速度）
nstlim=30000,
dt=0.002,
ntf=2,
ntc=2,
temp0=300.0, #恒温器的温度在300k
ntpr=100,
ntwx=100,
cut=8.0,
ntb=2, #使用常压的周期性边界条件
ntp=1, #对恒压模拟使用Berendsen气压计
ntt=3,
gamma_ln=2.0,
ig=-1,
/
四、跑amber动力学sander（sander是amber的通用MD引擎）
现在我们已经有了所有的元素:参数和拓扑文件prmtop，坐标文件inpcrd，以及输入文件01min。02 _heat。03 _prod。在此，我们准备运行实际的最小化、加热和生产MD。
1、运行最小化
$ $AMBERHOME/bin/sander -O -i 01_Min.in -o 01_Min.out -p prmtop -c inpcrd -r 01_Min.rst -inf 01_Min.mdinfo
-O #如果输出文件已经存在，就覆盖他
-I 01_Min.in #选择输入文件
-o 01_Min.out #写输出文件
-p prmtop #选择参数和拓扑文件
-c inpcrd #选择坐标文件inpcrd
-r 01_Min.rst #用坐标和速度写输出重新启动文件
-inf 01_Min.mdinfo #用模拟状态写MD信息文件
#####在sander完成之后，应该会有一个输出文件01_Min.out，一个重启文件01_Min.rst（将使用它继续为系统加热），和一个MD信息文件01_Min.mdinfo
Ps：提示prmtop不存在，我又做了一遍tleap
2、跑升温动力学 （使用初始最小化的重新启动文件，我们将对系统进行加热）
$ $AMBERHOME/bin/sander -O -i 02_Heat.in -o 02_Heat.out -p prmtop -c 01_Min.rst -r 02_Heat.rst -x 02_Heat.mdcrd -inf 02_Heat.mdinfo
-c 01_Min.rst #现在对于输入坐标我们从最小化中选择重新启动文件
-x 02_Heat.mdcrd #输出动力学模拟的轨迹文件
#####在heating完成之后，我们将得到三个文件：02_Heat.out 02_Heat.rst 02_Heat.mdinfo ps：文件的内容同最小化
我们打开02_Heat.out，我们能够看到系统信息，包括时间步长，能量和温度。举个例子（1000时间步长）：
NSTEP = 1000 TIME(PS) = 2.000 TEMP(K) = 30.20 PRESS = 0.0
Etot = -6945.8055 EKtot = 115.0391 EPtot = -7060.8446
BOND = 0.8441 ANGLE = 1.0099 DIHED = 9.0591
1-4 NB = 2.7105 1-4 EEL = 46.1439 VDWAALS = 1440.3323
EELEC = -8560.9443 EHBOND = 0.0000 RESTRAINT = 0.0000
Ewald error estimate: 0.6018E-03
NSTEP #所在的时间步长
TIME #模拟的总时间（包括重启）
TEMP #系统温度
PRESS #系统压力
Etot #系统总能量
Ektot #系统的总动力学能量
EPtot #系统的总势能
3、跑动力学
$ $AMBERHOME/bin/sander -O -i 03_Prod.in -o 03_Prod.out -p prmtop -c 02_Heat.rst -r 03_Prod.rst -x 03_Prod.mdcrd -inf 03_Prod.info &
### !!!注意：最后的&号表示让程序在后台运行（因为模拟时间过长），期间可以通过$ tail -f 03_prod.out 查看跑的进程，也可以通过$ cat 03_Prod.info 查看精确的系统数据和估算的完成时间。 Ctrl+c 结束tail
4,使用vmd观看轨迹
从amber跑出来的轨迹，官方推荐使用以.nc结尾的二进制格式，但是这种格式Windows版本的vmd识别不了，只有linux版本的才可以（NetCDF）。
我们需要先使用
cpptraj -p name.prmtop -y name.nc -x namd.dcd
将轨迹转换成DCD格式
再用vmd打开（CHARMM，NAMD，XPLOR DCDTrajectory）