Sensenbobo Blogging

关于两个非平衡分子动力学模拟的命名

分析蛋白质，核酸等生物分子结构机械特性时，常常可以在模拟中对这些分子添加外来拉力，以模拟原子力学显微镜（AFM）的做法。

现在这中模拟方法的使用越来越多，叫法也比较混乱。原则性上，可以分为常速拉伸分子动力学模拟（constant velocity pulling MD simulation) 以后常力拉伸分子动力学模拟（ constant force pulling MD simulation)。其实就是模拟AFM实验中力学探针拉伸生物分子的动作。

在英文名称中，很多人把所有这些受力模拟称为“Steered MD simulation"，这个称谓主要来自美国某为比较知名教授，本人觉得莫名其妙，同时十分难翻译。另外的一种称为与AFM联系起来，非别称常速拉伸与常力拉伸为”Force-probe MD simulation（FPMD）“和”Force-clamp MD simulation（FCMD）“。这种跟模拟原理相结合的称为比较合理，可以翻译为“力学探针分子动力学模拟”和力学“探钳分子动力学模拟”。

在此强烈建议各位帅哥靓妹使用这种与原理相结合的命名方法，不要在杂志上造成混乱。嘿嘿。

进模拟的基本步骤行

http://hi.baidu.com/picb

以下是进行分子动力学模拟的最基本步骤，来自Gromacs维基。具体的步骤可能会因为模拟目标的不同而存在差异，这里只做基本的概述。

1.

十分清楚要什么样的模拟系统性质和现象，明白模拟的目标。
2.

选择适当的工具以实现模拟目标。这需要熟知其他研究人员对相似模拟系统的模拟方法等资料，要好好读别人的论文，主要包括：

-模拟使用的软件，当然还要考虑使用的力场，有时因为力场关系，需要更改使用的软件。

-力场描述了系统中原子的性质和相互作用，选择一个适合自己研究目标的力场至关重要，还是需要多看论文。根据sensenbobo血泪史和别人的经验，模拟核酸的话，可以采用amber力场，水分子可以使用TIP3P；模拟蛋白质的话，可以采用OPLS-AA，水分子使用TIP4P（重申一下，我本人没有跑过测试）。
3.

获取模拟分子的坐标文件（PDB），或者自己生成一个。自己生成很好啊，强烈鼓励，但是要十分明白有一定的危险性，要明白每一步模拟系统的现象是为什么。
4.

根据模拟软件，生成模拟系统初步坐标文件。
5.

获取或者生成系统的拓扑文件，比如Gromacs的pdb2gmx命令，或者使用网络服务器PRODRG（google一下这几个字母），amber的leap命令和antechamber，NAMD的psfgen或者VMD。请参考相关手册。
6.

给模拟系统添加一个盒子（盒子的翻译方法我十分不喜欢，箱子显得太大，边界显得太难理解，但是盒子听起来想到蛋糕，我们是严肃的科学家。娃哈哈。）在盒子中添加水分子，添加离子如NA和CL之类（添油加醋，生活就是这样进行的）。添加任何东西到你的系统中，就要相应的修改一下拓扑文件，以保持与系统对应，不然会出错哦。
7.

进行能量最优话。这样做是消去系统中原子之间距离太近造成的高能量。必要的时候，分子可以现在真空中能量最优话，然后再做溶剂环境中能量最优话。
8.

选择恰当的模拟参数模拟平衡模拟系统。需要慎重考虑力场的作用，必要的时候保持系统体积温度不变（NVT系综）进行位置限制性平衡系统，然后在保持系统温度压强不变（NPT系综）修正系统的密度。然后再选择正确的系综（NPT，NVT，NET，BT？？）进行数据收集。（真烦！）
9.

跑足够长时间的模拟，使系统达到理想的平衡。使用模拟软件的工具或者其他软件观察系统的性质，如温度，能量瞪瞪。
10.

使用适当的参数进行数据采集模拟，注意模拟衔接过程中各个系统性质保持不变，如果原子速度之类。这里还要考虑力场的作用，同时要考虑怎么测量描述模拟目标。
11.

模拟足够长的时间使模拟目标明显出现（可能也出现跟你期望相反的结果，出现率应该百分之五十，加上个人感情因素，可能会达到百分之八十。保重！）。
12.

分析模拟结果，解释得到结果。要好好利用视图工具，同时要十分注意漂亮的三维图说明不了什么问题。图表，比较，统计，一个都不能少。

就这样吧。

LEAP是一个好东西！！

前天突然想找一个D-ALA的残基，想接一段含有D-ALA的肽段，太PDB库里有很多，但是最后还是leap帮我搞定啦。

leap是amber的分子建模软件，有tleap和xleap两个版本。 tleap是命令行界面；xleap是图形界面的。leap本身的命令很多，可以看看amber手册。dock软件现在也包含了leap了，不知道他们两家人怎么不打架。唉唉，看看国内的QQ告珊瑚虫，真的没法比。

好了，xleap很好用，要做一个D氨基酸很简单，把侧链移动一下就搞定。当然，用的是amber力场。只要你不要保存你修改过的残基，下次启动leap的时候，它会使用L残基。

嗯，好玩！

Gromacs的mdp文件参数浅析（一）

在使用Gromacs进行分子动力学模拟时，控制分子动力学模拟的参数文件称为mdp (Molecular Dynamic Parameters)文件。该文件中定义分子动力学模拟的各种行为，参数颇多。以下为各个参数的基本解析，详细请参考gromacs相关文件。

转载请注明出处，欢迎访问南开之星主页：www.nkstars.org

mdp文件分为几个部分，每次使用grompp读入mdp文件以生成tpr模拟文件时，总会生成一个mdout.mdp文件，在该文件中有你定义的参数，也有其他没有定义的默认参数。mdp文件各个部分如下：

预处理

title：即文件头，你想定义你模拟文件的名字，可以随意定义，是一个刷个性的地方。

cpp：预处理器，于C/C++的预处理器一样，默认为（/lib/cpp）

include：引用文件，即你拓扑文件中引用其他文件的路径，引用方式与C/C++引用一样：
格式：-I/home/../include -L/home/../lib

define；预定义，没有默认预定义。可以使用如何模拟的预定义方法控制模拟进程，gromacs目前手册可供选择的有以下两种（其实还有其他，可以查看邮件列表）：
" -DFLEX_SPC "：即让gromacs为SPC水模型引入软性性质。这样可以在进行能量最优话是使共轭梯度法产生作用，同时也可以更进一步进行最速下降法进行能量最优话。
" -DPOSRE "：让gromacs把posre.itp文件包含到拓扑文件中，这样做是为了进行坐标限制性动力学模拟。

模拟控制

integrator：动力学模拟方法，即整合牛顿力学定理的方法，根据模拟目的不同选择不同。
" md: "使用跳蛙算法(leap-frog)整合牛顿定律。
" sd: "另外一种跳蛙法统计整合。使用这个选项时，某个或者某些原子组(tc_grps)的温度设定为某特定温度(ref_t[K])，这些组运动反方向的摩擦常数可以设定为某一个值(tau_t[ps])。tcoupl参数在这个选项中被忽略。这个参数的随机算子由ld_seed设定。(NOTE: 这个方法中温度偏差的回原要比使用Berendesen热浴方法快一倍，即使使用相同的tau_t值。)
" bd: "使用Euler整合方法处理Brownian或者坐标Langevin动力学模拟，模拟中的粒子的速度为所受力除以摩擦因子(bd_fric[amu ps^-1)，加上一个随机的热力学噪音(bd_temp[K])。当bd_fric=0时，模拟粒子的摩擦因子为其质量除以tau_t，这与sd方法一致。随机算子由ld_seed指定。

以下几种模拟方法不是整合牛顿力学定律，但是也在此处指定，主要用于能量最优话模拟等。
" steep: " 使用最速下降法进行能量最优话，能量最优话最大位置移动为emstep[nm]设定，能量最大容忍度为emtol[kJ mol^-1nm^-1决定。
" cg: "使用共轭梯度法进行能量最优化，能量最大容忍度为emtol[kJ mol^-1nm^-1决定。在进行最速下降法能量最优化之后再进行一次共轭梯度法能量最优化是十分有效的能量最优化综合方法，可以使用nstcgsteep设定。在要对能量最优化进行常态分析时，最好使用双精度的GROMACS，以保证较高的精确度。
" nm: "对tpr文件中的系统结构文件进行常态分析。GROMACS必须为双精度。

tinit: 模拟计算开始时间，默认为0，即开始新模拟计算，单位为ps（该参数只对md，sd和bd有效）。

dt: 模拟步长，默认为0.001ps（只对md，sd和bd有效）。

nsteps: 整合算法的最大模拟步数。

comm_mode: 对系统或者系统中各个组质心的操作，有三种选项：
" Linear: "移除质心的平动;
" Angular "移除质心的转动和平动;
" No: "不对质心进行任何操作。

nstcomm：对质心进行操作的频率，默认为一个模拟步长，

comm_grps：对质心进行操作的组，可以是索引文件中的一个，或者多个组。

Gromacs的genion命令

在给蛋白质添加水环境之后，一般要在水环境中添加金属离子，使模拟系统更加接近真实系统。如果系统中蛋白质本身已经带了静电量，那么就更要给系统加几个带相反电量的金属离子，使系统处于电中性。

gromacs中添加金属离子的命令是genion，使用" genion -h "可以得到其使用的参数，其中有几个比较常用：
---------------------------
" -s: "指定系统tpr文件。
" -p: "指定系统拓扑文件，在往系统中添加金属离子时，genion会往拓扑文件最后的分子类型中写入添加的离子数，并修改拓扑文件中系统原子数。
" -o: "指定输出文件，genion的输出是pdb文件或者gro等结构文件。也就是说你产生这个文件之后，还要再用这个文件产生tpr文件。
" -np/-nn: "带正/负电金属粒子的数目。这个数目有一点讲究，一般需要看个人的应用。假如想要得到" 0.1mol/L " 的离子浓度到底要加多少，可以自己算一下（很简单，方法很多，比如看课本）。也可以直接使用" -conc " 参数直接指定离子浓度，在使用" -conc "参数时，建议使用" -neutral "参数配合，即使系统的最后处于电中性。嗯，gromacs开发组想得不要太周到哦（南京话）。
" -pn/-nn "指定正负金属离子的名字，比如" NA+ "或者" CL- "。可以看看gromacs安装途径" share/gromacs/top/ " 下面你用的力场文件中离子到底用什么名字，也可以使用新的离子，但是要在力场中定义，或者把新离子的itp文件使用" include "添加到系统拓扑文件中。
" -random "随机位置添加离子。这个比较有说法，如果不用该参数，那么离子就会添加在势能最低处，即靠近蛋白质相反电量的部位。有的人说这样可以节省时间，但是这样一跑MD，离子就会抱跟着蛋白，死死不放，不是很好。
" -seed "有随机，就是随机种子，如果发现使用" -random "添加离子自由，离子离蛋白太近（比如说小于0.1nm），那么可以指定新的seed。
---------------------------
嗯，给一个例子：
------
" genion -s topol.tpr -o system_ion.pdb -p system.top -np 100 -pname Na -nn 100 -nname Cl -random "
------
说明一下：加了100个Na和100个Cl，随机加，输出文system_ion.pdb文件。

再使用grompp生成一个新的tpr，就可以计算了。

Gromacs轨迹分析工具g_traj

几乎gromacs的所有分析数据都可以输出为xmgrace的数据文件，g_traj可以产生gromacs轨迹的各个组的坐标，速度，受力和边界等等。使用” -com “参数可以求出轨迹中各个组的质心的坐标，速度，受力等；使用” -mol “则可以求取系统中各个分子的信息；” -ot “则可以求出系统中各个组的温度。

还有几个其他参数，比如” -cv “可以求平均速度，” -cf “可以求平均受力等。其他参数同其他命令无异，参加说明文件。

Gromacs的g_energy另外一种用法

Gromacs的各个工具都很有个性，如果互相结合，可以做很多事情。

g_energy求系统轨迹各个能量的，一般跑完MD之后，使用g_energy处理ener.edr只能得到系统的各个能量项。但是如果想求系统中两个不同部分在模拟过程中的相互作用能量，那就要使用一些小窍门。

以下是实现的一个方法：
第一，根据原来的tpr文件建立一个新的tpr，在这个新的tpr中，明确定义感兴趣的组。这要用索引文件，见上文。
第二，用mdrun的" -rerun "参数指定原来的轨迹文件再跑一次模拟，这个过程很快。如果还想更快，可以使用trjconv把水分子去掉。这一个重复的模拟也产生轨迹文件，重要的是，还产生一个新的ener.edr文件，这个文件中包含了tpr文件中定义的各个组能量及相互作用能量（库伦相互作用能，范德华相互作用能等）。
第三，再使用g_energy把各个能量项提出来，想要什么提什么。

嗯，结果非常好。不信你试试。

使用grompp提取上一次模拟最后速度和能量

在上面提到的使用gromacs程序包中tpbconv命令制作新的tpr文件中，最后提到新制作的.tpr文件只能使用跟原来.tpr文件一样多的CPU数目。还抱怨说这是tpbconv一个不足的地方。使用grompp可以制作一个新的.tpr文件，从上一步模拟的轨迹文件中提取速度，并从上一步能量文件中提取能量，也可以无缝的链接重启模拟计算。

要做到从上一步的最后的一个系统状态开始新的模拟计算。首先要在.mdp文件中把“ gen_vel ”参数定义为" no "，这样做是为了告诉grompp不要重新为系统中的原子指定随机速度。指定新模拟开始的时间，即修改" tinit "参数。然后可以使用一下命令制作一个从上一步模拟文件中提取速度和能量的.tpr文件：
----------------------------------
grompp -f [.mdp文件] -c [上一步模拟最后的系统坐标文件] -p [拓扑文件] -t [上一步的trr轨迹文件] -e [上一步能量文件] -time [坐标文件对应的模拟时间] -o [输出tpr文件] -np [CPU数目]
----------------------------------
提取上一步模拟系统的速度时使用trr文件，是因为xtc为单精度，没有trr文件精确。" -time "参数告诉grompp在上一步模拟文件中提取该时间的能量和速度，所以该时间要和系统的坐标文件相一致。

看起来好像要比tpbconv命令复杂一点，但是可以改变CPU数目，还算十分灵活。Gromacs是灵活的人的MD工具。^_^

Gromacs索引文件概述

Gromacs的索引文件，即index文件，由make_ndx程序生成，文件后缀为.ndx。

索引文件是gromacs最重要的概念文件，使用它可以在模拟过程中为所欲为。举一个简单的例子，比如想详细了解HIV整合酶切割DNA的反应机理，使用量子力学模拟反应位点的反应过程，而分子其他部位使用一般分子动力学模拟。于是我们就面临一个对模拟系统进行分割定义的问题，在gromacs中，就要用到索引文件。

基本的思路是这样的，在索引文件中，定义一个独立的组，这个组包括反应位点处所有原子。在模拟的.mdp文件中，对这个组定义量子力学模拟，事情就是这么简单。对蛋白进行量子力学模拟时，一般使用洋葱模型。所谓洋葱模型，就是对反应位点使用量子机制，在反应位点一定的半径内，使用半量子力学机制，然后分子部分使用分子机制。那么索引文件就定义一个使用量子力学的组，把需要引进量子机制的原子都放到这个组中；再定义一个半量子机制的组，同时放进需要半量子力学机制模拟的原子，再在.mdp文件中独自定义即可。

再举一个例子，比如说在进行SMD（Steered Molecular Dynamics，这个我一直没有想到或者找到切恰的中文翻译方法，或许可以叫做牵引分子动力学？？别扭！！）中，要对蛋白莫一个原子或者残基作用力，那么可以建立一个索引文件，在该文件中定义一个组，把要施力的残基或者原子放到该组中。然后在.ppa文件中使用该组就行了。

如果我还没有说明白，那么看看gromacs的参考文件吧。如果还是不明白，可以来找我，我免费培训。^_^

索引文件使用make_ndx命令产生，"make_ndx -h"可以看到全部的参数。运行make_ndx后，可以使用" r "命令选择残基，" a "命令选择原子，" name 命令多组进行改名。可以使用" | "表示或运算，" & "表示与运算。下面是几个简单的例子：
-----------------------------
1.选择56号残基
r 56
2.选择3至45号残基
r 3-45
3.选择3至15，23至67号残基
r 3-15 | r 23-67
4.选择3至15号残基的主干链原子
r 3-15 & 4 #在索引文件中，4号组为默认的主干链。
-----------------------------
组合是灵活的，使用的时候好好发挥聪明的大脑啦。

个人觉得gromacs把索引文件概念做得非常好，并独立成一类文件是一个不小的创举啊。在这个概念很值钱的年代，基本上使gromacs多一个大大的卖点。（再说明一下，gromacs是免费的GNU产品，质优无价。）

使用genbox命令为Gromacs模拟分子添加水环境

使用editconf把分子放进一个盒子之后，接下来要做的就是往盒子里面添加水分子。Gromacs中添加水环境的命令是genbox，这是一个较其他命令简单一点的命令，因为参数不多。

通常用到的genbox参数有一下几个：
-------------------------------------------
-cp ：带盒子参数的分子坐标文件，也就是editconf的输出文件；
-cs ：添加的水分子模型，如spc216、spce、tip3p、tip4p等，关于各个模型的区别，请参考scholar google；
-o ：输出坐标文件，就是添加水分子之后的分子坐标文件，默认是.gro文件，但是也可以输出其他文件格式，如pdb；
-p ：系统拓扑文件，genbox会往里面写入添加水分子的个数，这个不要忘记，不然在进行下一步计算时，会出现坐标文件和拓扑文件原子数不一致的错误；
-ci ：索引文件，genbox可以为分子特定部位添加水环境，这样可以减少原子数，只有研究的分子部位添加水环境，节省计算时间。
-seed ：随机种子数，添加水分子时，各个水分子的位置是随机的，可以改变这个随机数子使水分子重新分布。
-------------------------------------------

添加水分子后可以用VMD等软件看看结果，一般完成这个之后事情开始变得复杂。比如某一个水分子出现在蛋白结构中，而这个位置本来就是不希望水分子一开始就存在，那么可以找出这个水分的残基标号，进行删除，同时删除拓扑文件中水分子的数目等。

使用tpbconv重启gromacs模拟

在使用gromacs的mdrun进行模拟计算过程中，很多因素可以是模拟计算终止。比如突然断电，断网或者磁盘空间满，或者windows死机（^_^）等等。重启gromacs模拟计算是一件十分方便的事情，因为gromacs众多的程序里面就有一个专门（或者吧）用来修改tpr文件的，就是tpbconv。

gromacs把模拟需要的所以文件都打包成一个tpr二进制文件，里面包含了分子坐标，各个原子在给定温度下速度和能量的分布。当模拟突然终止时，只要将终止时候系统的状态，即各个原子的位置、速度、坐标等装入tpr文件即可。tpbconv的参数也不少，可以使用"tpbconv -h "查看，但是制作一个重启tpr文件的参数和格式一般如下：

tpbconv -s topol.tpr -f traj.trr -e ener.edr -o newtopol.tpr

其中topol.tpr为原来的tpr文件，traj.trr为双精度坐标文件（不要用xtc文件，因为精度不够），ener.edr为系统能量输出文件，newtopol.tpr是重启模拟文件。以上的命令得到的是在计算突然终止前一个系统构象的信息。也可以在命令中加上一个"-time "参数来指定从那一个时间重新开始，如一下指定从一纳秒处重新开始模拟：

tpbconv -s topol.tpr -f traj.trr -e ener.edr -time 1000 -o newtopol.tpr

同时，如果模拟正常结束，而模拟时间让人觉得不够长时，可以使用tpbconv写一个延长模拟的tpr文件，一般格式如下：

tpbconv -s topol.tpr -f traj.trr -e ener.edr -extend 1000 -o newtopol.tpr

其中"-extend 1000"表示延长1000ps的模拟时间。呵呵，非常好用。

这样断了又开始，就会产生很多轨迹文件，分析的时候非常不方便，gromacs有其他常用的命令把坐标文件，能量文件连接成一个文件，其中比较常用的如trjcat和eneconv，格式分别如下：

trjcat -f traj1.trr traj2.trr.... -o traj_all.trr
eneconv -f ener1.edr ener2.edr... -o ener_all.edr

即使用"-f "读入所有轨迹或者能量文件，使用"-o "输出完整的轨迹和能量文件。

最后说说一个tpbconv的弱点。tpbconv不能更改你原来tpr文件中并行计算的节点数，比如你原来的tpr文件是8个节点的，那么使用tpbconv得到的重启tpr文件也是8个节点的。如果想更改使用节点数，那只能用grompp重新做一个了。但是使用grompp做重启模拟文件时，就算你指定了原来的轨迹文件和能量文件，它还是会根据麦克斯韦分布重新给各个原子指定速度，真气人。

嗯，如果你觉得这是一个大问题，那就伸长脖子等gromcas新版本出来吧。

Gromacs的editconf命令

在使用pdb2gmx创建模拟分子系统之后，可以使用editconf为你的分子画一个盒子。也可以认为使用editconf把分子放进一个盒子中，这样，你就可以往盒子里面添加水分子，离子，或者其他溶剂等等了。

使用“ eidtconf -h ”可以看到editconf的参数，其中比较常用的有以下几个：
----------------------------------------------------------
-f    指定你的坐标文件。
-n    分子系统的索引文件。索引文件是gromacs一个十分突出的功能，刚接触有点复杂，但是功能相当强大。
-o    输出文件，即放进盒子里面的分子系统。
-bt    盒子类型，有正方型，长方形，八面型等等，看个人需要跟癖好啦。
-box    自定义盒子大小，需要三个长度，即X、Y、Z三个方向的长度。
-d    分子离盒子表面的最短距离。这个跟-bt一起使用，基本就足够了;如果蛋白在模拟过程尺寸变化很大，那就用-box吧。
-center    确定哪一部分分子放在盒子中心，如果和索引文件一起使用，可以非常详细的定义分子的位置。
-translate    平移分子，跟X、Y、Z三个方向，随便移动。可以参考我前面一篇关于加大盒子的e文文章，练习e文，不好意思。
-rotate    转动分子，还是三个方向。我觉得editconf开发者太有才，什么都想到了，我想要什么，他就有什么，以后遇到他，我会让他给我签名。
-princ    这个参数可以用来对齐分子，比如使分子沿X轴对齐。举一个例子吧，比如你想将分子中两个残基沿Y轴对齐，那么就在索引文件中将这俩个残基标记以下，然后使用-princ，根据提示走就能对齐分子啦。
----------------------------------------------------------

由于gromacs的很多命令都可以接受不同的文件类型，editconf也有其他功能，如和trjconv一样进行gro和pdb的转化：
editconf -f sen.pdb -o sen.gro
或者
editconf -f sen.gro -o sen.pdb

我觉得我要好好看语文，都并不知道怎么结束这样一篇实用性这么强的作文。悲哀~~

Gromacs的pdb2gmx命令

使用gromacs做分子动力学模拟时，第一个要用到的命令一般都是pdb2gmx。这个命令吧pdb分子文件转化成gromacs独特的gro分子结构文件类型，同时产生分子拓扑文件。

Gromacs是典型的GPL软件，每一个命令都有很多命令参数。这对熟悉windows环境的人来说有一点烦，但是如果熟悉了Linux环境，也就慢慢喜欢啦。（建议多使用命令，就像VMD， Pymol， rasmol和Chimera等等分子可视化软件，如果接合命令使用，功能都非常强大。另外一个比较bt的软件叫做WHATIF的，完全建立在bt的命令菜单上，心理承受能力不强者多半吐血而终。

使用“pdb2gmx -h”可以得到pdb2gmx的所有参数及简单说明（gromacs的任何命令都可以使用-h参数得到类似帮助）。pdb2gmx的参数很多，但是常用的只有以下几个：
-----------------------------------------------------------------------
-f    指定你的坐标文件，可以是pdb、gro、tpr等等包含有分子坐标的文件;
-o    输出文件，也就是处理过的分子坐标文件，同样可以是pdb、gro、g96等文件类型;
-p    输出拓扑文件。pdb2gmx读入力场文件，根据坐标文件建立分子系统的拓扑;
-water    指定使用的水模型，使用pdb2gmx的时候最好加这个参数，不然后面会吃苦头。它会提前在拓扑文件中添加水分子模型文件;
-ff    指定力场文件（下文讨论），也可以不用这个参数，再自行选择;
-ignh    舍弃分子文件中的H原子，因为H原子命名规则多，有的力场不认;
-his    独个指定HIS残基的质子化位置。
------------------------------------------------------------------------
其他的参数还不少，可以好好看一个pdb2gmx的帮助文件，一般的pdb2gmx的执行格式如下（假设你的分子坐标文件为sen.pdb）：

pdb2gmx -f sen.pdb -o sen.gro -p sen.top -water tip4p -ignh -his

该命令读入分子文件，使用tipp水模型，等等，然后pdb2gmx会让你选择力场文件。然后它就很聪明的帮你建立初始模拟系统啦。

gromacs自带的力场有很多：
------------------------------------------------------------------------
0: GROMOS96 43a1 force field
1: GROMOS96 43b1 vacuum force field
2: GROMOS96 43a2 force field (improved alkane dihedrals)
3: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
4: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
5: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
6: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
7: [DEPRECATED] Gromacs force field (see manual)
8: [DEPRECATED] Gromacs force field with hydrogens for NMR
9: Encad all-atom force field, using scaled-down vacuum charges
10: Encad all-atom force field, using full solvent charges
------------------------------------------------------------------------
本人建议使用OPLS力场和tip4p水模型。tip4p水有四个粒子，分别是两个氢原子，一个氧原子和一个没有质量的电粒子。这个电粒子在其他三个原子中间靠近氧原子。tip4p水模型多了一个粒子，模拟代价高一点，但是结果要好一点。但是最近有报道说使用spce水模型和GROMOS力场的计算结果最好，嘿嘿，好一个百家争鸣的学术氛围啊，我真的号感动哦。Gromacs也可以使用其他力场，如AMBER力场等，使用方法请参考google。

pdb2gmx的输出基本可以做真空中模拟了，现在对MD模拟的要求高，一般都要有点水。为分子系统添加水环境和离子环境需要其他命令，要慢慢来。

化学键摩尔斯势能

在经验力场中，化学键势能表达方法是使用弹簧势能表达的。相对来说，摩尔斯化学键势能表示方法要准确一些。摩尔斯化学健势能表达为：

U(x) = D_e[1 - e^-(x-r_e)]²

其中，D_e为两个成健原子互相分离时的化学能量; x为两原子之间的相对距离; r_e为两原子的平衡位置。
摩尔斯势能图可以简单表示如下：

Adding Ions to MD systems

Before preforming you MD simulation, proteins should be solvated into a water box. Try not to run MD simulation in vacuum these days, update you computer or just apply for some computational resources to Nankai Stars. ^_^

If you use gromacs, you can use genion to add ions to you system. Actually, genion is quite stupid. It just place the ions around the protein, because it thinks that is the low potential energy solution. you can use the -random flag to place the added ions all around the system. if you find out some ions are too close to the protein, change the -seed flag until you think it is reasonable.

Amber suite have two commands to add ions into MD system: addions and addions2. Both of them have some drawbacks. addions just add ions around the protein, while addions2 add ions to the system at the edge of the bound box. I think addions2 is more reasonable, cause it does not cause some artificial effects. But if you add some positive ions and some negative ions at the same time using addions2, both positive and negative ions will be at the edge of the bound box, then after the energy minimization, they will bound to each other. that is too bad.

Too solve the addions/addions2 problems in amber suite, you can use addions to add some ions at first, and use addions2 at the second time to neutralize the system. Or you can solvate your protein a little bit, then use addions2 to add some ions, then fully solvate your system and use addions2 to neutralize the whole system.

And how about Namd? you can tell me. :D

How to Enlarge Your System Using Gromacs

When doing SMD, sometimes you do not want to do the energy minimizaton(em) and position restrain(pr) with a large system. So, you just solvate your protein into a water box a little bigger than your protein, like with 1.5 nanometer buffer distance.

So here is the question: what do you do after the free MD run before starting to pull? You can take you protein out of the water box and solvate if into a bigger water box of course, and do em and pr again. But still, you can do it in another way.

First of all, use the -translate flag of editconf to move your system, like 5 nm from it's original position, and elongate the side of the water box along which you want to pull your protein( You have your protein oriented right??). After that, use genbox to solvate all your system. ( Or you do not have to use the -translate flag, just use editconf to enlarge your water box.)

Second, use genion to add ions to your system. This is the most complicate party of the job. You can do it step by step like the following:

1. Change all the newly added water molecular in the .gro file to another group name, like WAT.

2.Edit the .top files. At the end of the .top file, you can see genbox have automaticly name your new added water molecular in a new SOL group just like the other water molecular. Change this SOL to WAT.

3. Make a topology file for residue WAT. Take tip4p for example, copy the tip4p.itp to some file name like tipppp.itp and change all SOL in it to WAT.

4.The last step, edit you system topology file to include tipppp.itp.

OK, You now can use grompp to make a tpr file and use genion to add ions to your system. Remember to choose the WAT group to ions.

And now, you can do much shorter time of em and pr , and just pull you protein. ^_^

New Blog, Old man......

I register a new blog site on Hi.baidu.com: Sensenbobo Hi Baidu.

You can find something help about basic MD and other thing there, and here will still be my home. ^_^

GROMACS程序DEMO例程

####################
### 概述 ###
####################
-----------------------------------------------------------------
-----------------------------------------------------------------
该例程来自Gromacs程序/share/tutor/目录下。整个例程大概只需要十分钟就可以完成，非常适合初学者学习。

该例程是一个完整的分子动力学模拟过程，涵盖了Gromacs程序基本的使用方法。模拟内容是一个水环境下的小肽链。模拟唯一要求是该小肽链的PDB文件。

文档翻译如果有错，请你给我发信：sen@nkbbs.org。相关内容请参阅Gromacs文档，或者给Gromacs开发组询问。

Gromacs官方网站：http://www.gromacs.org
Email: gromacs@gromacs.org
-----------------------------------------------------------------
-----------------------------------------------------------------

#########################
### 环境变量设置 ###
#########################
-----------------------------------------------------------------
-----------------------------------------------------------------
在以下的例程中，所有命令都直接运行，没有添加绝对路径。所以，必须将Gromacs安装路径的bin文件夹加入到系统PATH变量中。如果不加入PATH变量，那么运行时要加入命令的绝对路径。
-----------------------------------------------------------------
-----------------------------------------------------------------

################
### PDB2GMX ###
################
-----------------------------------------------------------------
-----------------------------------------------------------------
在进行如何分子动力学模拟之前，必须建立分子的拓扑文件。Gromacs的分子拓扑文件是用pdb2gmx命令生成，文件后缀名为 .top。攑db2gmx的输入唯一文件是分子的PDB文件，可以从www.pdb.org寻找，文件后缀名 .pdb。”

不是所有的PDB文件都含有氢原子，pdb2gmx将添加所有缺失的氢原子。在pdb2gmx输出的gromacs结构文件中，包含了蛋白质结构的每一个原子，并定义了分子结构的尺寸大小。文件后缀名为 .gro。”
-----------------------------------------------------------------
-----------------------------------------------------------------
假设分子PDB文件名为MOL.pdb，命令运行格式：

pdb2gmx -f MOL.pdb -o MOL.gro -p MOL.top > output.pdb2gmx

命令得到三个文件：MOL.gro是gromacs结构文件；MOL.top分子拓扑文件；output.pdb2gmx是pdb2gmx命令输出文件。

###############
### GENBOX ###
###############
-----------------------------------------------------------------
-----------------------------------------------------------------
真空条件下进行分子模拟输出结果误差较大，所以模拟之前，必须给分子添加水环境。Gromacs中使用genbox命令完成。

Genbox命令读入gromacs的结构文件，并读入水盒子尺寸的大小，输出文件包含分子文件和水盒子。同时genbox更改原来的拓扑文件，使其包含水分子。在使用genbox之前，要使用editconf命令定义水盒子大小。

-----------------------------------------------------------------
-----------------------------------------------------------------
命令格式：

editconf -f MOL.gro -o MOL.gro -d 0.5 > output.genbox

genbox -cp MOL.gro -cs -o MOL_b4em.gro -p MOL.top >> output.genbox

##############
### GROMPP ###
##############
-----------------------------------------------------------------
-----------------------------------------------------------------
原则上讲，到此为止已经可以开始进行模拟计算了。但是，添加在准备分子系统过程中，系统有很多原子距离太近，局部能量太高。这些相互距离太近的原子多是由于genbox程序产生的，溶剂分子与蛋白分子之间存在不稳定的高能量。如果现在开始模拟计算，而不进行能量最优化，系统将可能很不稳定。

去除这些局部高能量的办法是对系统进行能量最优化。能量最优化过程是改变系统中局部高能量的原子的位置，降低这些点的能量。

在进行能量最优化之前，我们先用GROMACS的预处理程序grompp处理所有输入文件。grompp预处理拓扑文件（.top），坐标文件（.gro）和一个参数文件（.mdp），然后输出一个二进制拓扑文件（.tpr）。这个二进制文件包含所有模拟需要的信息，利用这个文件即可以进行能量最优化和动力学模拟。
-----------------------------------------------------------------
em.mdp文件范例，详细请参考Gromacs手册
-----------------------------------------------------------------
title = MOL
cpp = /usr/bin/cpp
define = -DFLEX_SPC
constraints = none
integrator = steep
dt = 0.002 ; ps !
nsteps = 100
nstlist = 10
ns_type = grid
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
;
; Energy minimizing stuff
;
emtol = 1000.0
emstep = 0.01
---------------------------------------------------------------------------------

Grompp命令格式：
grompp -f em -c MOL_b4em -p MOL -o MOL_em >&! output.grompp_em

################
### MDRUN EM ###
################
-----------------------------------------------------------------
-----------------------------------------------------------------
有了二进制拓扑文件，现在可以开始进行能量最优化了。进行能量最优化的程序是mdrun，在
Gromacs程序中，所有模拟都是用mdrun程序进行。

在进行能量最优化过程中，注意查看mdrun程序的输出文件。在输出文件在中，从左到右第一个数字是模拟步数，第二个数字是计算步长，第三个数字是系统能量。如果模拟顺利利，可以看到系统能量从一个很高的值迅速降低，最后稳定在一个大负值。
-----------------------------------------------------------------
-----------------------------------------------------------------
命令格式：
mdrun -nice 4 -s MOL_em -o MOL_em -c MOL_b4pr -v >& ! output.mdrun_em

#################
### GROMPP PR ###
#################
-----------------------------------------------------------------
能量最优化之后，首先保持蛋白质不动，对蛋白质周围水环境进行动力学模拟，该过程称为位置限制性分子动力学（position restrained MD）

位置限制分子动力学保持蛋白质位置不变，对溶剂分子进行平衡计算，可以使溶剂分子填补空间空洞。这个可能存在的空洞是genbox程序产生的。

首先对输入文件进行预处理得到二进制拓扑文件，这些输入文件就是能量最优化得到的输出文件。然后再写一个参数配置文件和索引文件。

默认情况下，程序对模拟系统是分部分的。我们利用两个部分进行位置限制性模拟：Protein和
SOL部分，分别表示蛋白质和溶剂。这这个过程中，即保持protein位置不变。

参数文件（.mdp）包含了位置限制性模拟的所有参数，包括步长，步数，温度等等。同时参数文件告诉GROMACS模拟的类型，如能量最优化、位置限制性或者是分子动力学模拟。
-----------------------------------------------------------------
pr参数文件
-----------------------------------------------------------------
title = MOL position restraining
cpp = /usr/bin/cpp
define = -DPOSRES
constraints = all-bonds
integrator = md
dt = 0.002 ; ps !
nsteps = 500 ; total 1.0 ps.
nstcomm = 1
nstxout = 10
nstvout = 1000
nstfout = 0
nstlog = 10
nstenergy = 10
nstlist = 10
ns_type = grid
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
; Berendsen temperature coupling is on in two groups
Tcoupl = berendsen
tau_t = 0.1 0.1
tc-grps = protein sol
ref_t = 300 300
; Pressure coupling is not on
Pcoupl = no
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocites is on at 300 K.
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529

----------------------------------------------------------------------------------------
命令格式

grompp -f pr -c MOL_b4pr -r MOL_b4pr -p MOL -o MOL_pr >& ! output.grompp_pr

################
### MDRUN PR ###
################
-----------------------------------------------------------------
-----------------------------------------------------------------
现在进行1ps的位置限制性分子动力学模拟（the Position restrained Molecular Dynamics simulation）在这里进行1ps只是为了节省例程演示实践时间，在实际研究中，1ps模拟将太短，不符合实际情况。

-----------------------------------------------------------------
-----------------------------------------------------------------
命令格式
mdrun -nice 4 -s MOL_pr -o MOL_pr -c MOL_b4md -v >& ! output.mdrun_pr

#################
### GROMPP MD ###
#################
-----------------------------------------------------------------
-----------------------------------------------------------------
到此为止，分子系统已经可以进行正式的分子动力学模拟了。再一次利用grompp程序将输入文件转化成二进制拓扑文件（.tpb/.tpr文件后缀）。
分子动力学模拟配置文件（参数详细意义请参考Gromacs文档）
-----------------------------------------------------------------
title = MOL MD
cpp = /usr/bin/cpp
constraints = all-bonds
integrator = md
dt = 0.002 ; ps !
nsteps = 5000 ; total 5 ps.
nstcomm = 1
nstxout = 50
nstvout = 0
nstfout = 0
nstlist = 10
ns_type = grid
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
; Berendsen temperature coupling is on in two groups
Tcoupl = berendsen
tau_t = 0.1 0.1
tc-grps = protein sol
ref_t = 300 300
; Pressure coupling is not on
Pcoupl = no
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocites is on at 300 K.
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529
----------------------------------------------------------------

命令格式：

grompp -f md -c MOL_b4md -p MOL -o MOL_md >& ! output.grompp_md

################
### MDRUN MD ###
################
-----------------------------------------------------------------
-----------------------------------------------------------------
现在进行分子动力学模拟计算。详细观察模拟步数的增加。（整体模拟步数为2500，时间长度为5ps。再一次强调，模拟时间太短，仅适合教学演示。）
-----------------------------------------------------------------
-----------------------------------------------------------------
mdrun -nice 4 -s MOL_md -o MOL_md -c MOL_after_md -v >& ! output.mdrun_md

############
### NGMX ###
############
-----------------------------------------------------------------
-----------------------------------------------------------------
最后，可以使用Gromacs的ngmx程序观看动力学模拟轨迹（文件后缀.trj）。轨迹文件包含了分子动力学过程的坐标，原子速度和受力情况等等。可供分析使用。Ngmx的详细使用方法，请参考Gromacs文档。由于分子动力学模拟分析方法各异，这里不再讨论。

本文来自gromacs程序教学例程，仅供参考。^_^

分子动力学模拟概论

五、经验力场

对生物分子的理论研究主要是从原子水平探究其结构、功能与动力学之间的关系。这样的问题涉及大量原子，目前还无法直接使用量子力学解决。但是如果使用经验力场，问题则可以很好的解决。所谓经验力场，简单的说就是使用一系列势能方程描述原子之间的相互作用。使用经验力场解决生物分子中各个原子之间的相互作用，可以节约大量计算时间，但是是牺牲一定准确性为代价的。

经验力场发展到今天，已经很好的在计算精确度和计算效率之间找到一个平衡。这些经验力场一般经过一定使用时间后还不停与实验结果进行校准与修改。经验力场对物理实验的再现能力主要包括：X光衍射和NMR得到的生物结构信息、光谱学方法和无弹性中子散射得到的动力学信息、统计热力学信息等等。经验力场方程的各个参数的校订是一件费时费力的工作，需要不停的大量优化。在过去二十年中，经验力场研究已经形成一个独立的领域，不同的研究小组还不停的加入到这个领域中去。现在比较著名的经验力场有AMBER力场、CHARMM力场、GROMOS力场和OPLS/AMBER力场等。这些力场的持续发展重视基础研究同时更重视实际应用。

上文提到，使用经验力场有一定的限制性。经验力场最明显也最重要的一个限制性就是模拟的生物分子不能发生重大构象变化，如化学键的断裂和重建。为了突破这样的限制性，新型的经验力场已经开始引进量子力学机制。这一类力场还需要一定的完善时间，如有兴趣可以参阅其他资料。

经验力场之间一般大同小异，并一般公开使用，如NAMD软件使用CHARMM力场，也可以使用AMBER力场；DOCK实用部分AMBER力场最为评分函数。在此简单介绍CHARMM22力场，该力场可以用于模拟蛋白(MacKerell et al 1998)、核酸(MacKerell et al. 1995)、流体(Schlenkrich et al. 1996)和碳水化合物(Ha et al. 1988)等。

References
A. D. MacKerell Jr., J. Wiórkiewicz-Kuczera, and M. Karplus (1995). An All-Atom Empirical Energy Function for the Simulation of Nucleic Acids, J. Am. Chem. Soc. 117, 11946-11975.

A. D. MacKerell, Jr., D. Bashford, M. Bellott, R. L. Dunbrack Jr., J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, III, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin, M. Karplus (1998). All-atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. 102, 3586-3616.

M. Schlenkrich, J. Brickmann, A. D. MacKerell Jr., and M. Karplus (1996). An Empirical Potential Energy Function for Phospholipids: Criteria for Parameter Optimization and Applications, in Biological Membranes: A Molecular Perspective from Computation and Experiment, K. M. Merz, Jr. and B. Roux, editors (Birkhauser), 31-81.

S.N. Ha, A. Giammona, M. Field, J.W. Brady, A revised potential-energy surface for molecular mechanics studies of carbohydrates, Carbohydr. Res. (1988), 180(2), 207-221.

---

CHARMM经验力场

一个系统的能量E，可以表达为一个关于系统中所有原子笛卡尔坐标R的方程。系统能量值主要由两个部分组成：分子内部化学键、键角、二面角相关的能量值E_bonded；外部非成键相互作用E_non-bonded，如静电和范德华作用等。

其中E_bonded部分主要有三个能量相组成：



分别代表键长畸变、键角畸变和二面角畸变。

---

上述等式第一项键长畸变能是以弹簧振子弹性势能方程近似表达的，其中b₀表示成键原子之间的平衡位置，K_b为力学常数。这两个参数都是根据原子性质不同而不同的。






力学常数常常需要根据红外线频率延展实验或者量子力学计算进行调整，而键长平衡值则多依据高分辨率晶体结构或者微波光谱分析数据进行调整。

---

第二项键角畸变是指化学键角偏离平衡角度θ₀所具有的弹性势能。平衡键角θ₀和力学常熟K_θ根据不用化学键角不同而不同。






以上两项能量都是分子结构偏离平衡位置所产生的势能，在一个平衡分子体系中，该两项理想状态下总合应该接近零。

---

第三项是涉及二面角畸变所具有的势能，包含三个连续化学键连接的四个原子。二面角畸变能是化学键旋转造成的，势能变化具有周期性，一般使用一个余弦方程表示。

除了以上三项，CHARMM力场还包含其它两项与化学键相关的能量项：
i Urey-Bradley作用能：相隔两个化学键原子之间的相互作用；
ii 非正常二面角作用能：用于手性分子以及保持原子构成的平面。
而力学常数：K_b、K_q、K_f主要是对大量小分子复合物研究得来，并同时与量子力学计算结果相比较调整。

---

ChARMM力场中对非键对相互作用对势能的影响的考虑主要分两个部分：范德华相互作用能和静电相互作用能。在其他经验力场中，有的还独立计算氢键对势能的影响，而CHARMM力场则将氢键包含在范德华相互作用和静电相互作用中。



两个原子之间的范德华相互作用力主要是在两个原子距离小于平衡距离时的排斥力和距离大于平衡距离时表现出来的吸引力。排斥力与吸引力都与原子的电子云分布有关，详细内容可以参考相关资料。当排斥力与吸引力大小相同，互相抵消的时候，两原子的距离为平衡距离，总体势能最小。一个系统中原子之间都处于平衡位置时，系统最为稳定。系统最低能量值和原子之间平衡距离大小取决于不同原子类型。



范德华相互作用是保持生物分子系统稳定最为重要的因素，其计算模型一般为Lennard-Jones 6-12势能模型：



常量A和C主要取决于原子类型，可以通过如气相光散射的实验测量。
原子之间的静电相互作用可以使用库仑电势能表示：



其中D为有效电介质常数，r为两个原子电量中心之间q_i和q_k之间的距离。静电相互作用根据系统中各个原子的坐标位置可以各个求值，同时可以确定各个原子的受力情况以应用于分子动力学模拟计算。

---

经验力场是一种能量近似，存在一定不精确的局限性。

经验力场中各个原子的属性是同一的，原子属性用于定义与其它原子成键或者相互影响时表现的特性。但是，比如一个在sp3电子成键状态下成键的C_&alpha原子与一个HIS残基环结构上的C原子的性质是不同，这不是经验力场所能描述的。除了对各个原子的原子属性进行统一处理，经验力场还对不同原子进行分类，以尽可能减少原子类型，简化计算，这回导致原子类型识别错误。某一些原子，如脂肪质中的碳原子或者生物分子中的氢原子对周围环境比较不敏感，这类原子的统一属性定义一般不会出现较明显的误差；但是如氧原子和氮原子与周围其它原子影响十分强烈，这一类原子属性需要根据成键环境不同而不同。

为了降低计算资源的要求，经验力场对生物系统中势能计算使用原子势能成对叠加近似，即系统中一个原子与周围的相互作用是该原子与其他各个原子相互作用的叠加，两个原子自间的相互作用忽略其他原子的影响。三个及三个原子以上之间的同时间的相互作用在经验力场中是不予考虑的，所以极化作用并不能精确的表述。这些都有可能使计算结果与实验结果存在差异。

另外重要一点是经验力场势能方程并没有考虑熵效应，所以势能的最小值并不一定与系统稳定状态或者最有可能处在的状态互相对应，这些状态在实验条件下是与系统势能最低态相对应的。

Namd 2.5与Namd 2.6安装说明

安装简单指南：
Installation Guide

分子动力学模拟概论

分子动力学模拟概述

四、经典力学

分子动力学模拟方法是建立在牛顿第二定律之上的，即牛顿运动方程：

F=ma

其中F为作用在粒子上的力，m表示粒子的质量，a表示粒子的运动加速度。只要知道作用在每一个粒子上的力，即可知道系统中每一个粒子的加速度。利用这个方程，模拟过程记录系统中各个粒子的位置、速度和加速度随时间变化的系统轨迹文件；而从这个轨迹，即可以求出系统各个性质的平均值。这是一个确定性的方法，系统中的每一个粒子的任一时刻位置跟速度一旦知道，系统过去或者未来的性质都可以推算得到。分子动力学模拟方法虽然要求计算时间长，计算资源多；但是计算机已经变得越来越快，价格也越来越便宜。现在的分子动力学一般在纳秒级别，毫秒级别还鲜有报道。

--------------------------------
牛顿运动方程：

其中，F_i表示作用在质点i上的力，m_i表示质点i的质量，而a_i表示质点i的加速度。作用力可以表达成势能的梯度函数，即

将上面两个式子进行整合，可得

其中V表示系统的势能。牛顿运动方程即表达成系统势能与随时间变化的质子位置之间的表达式。

--------------------------------
牛顿第二定律的一个简单应用：

如果加速度恒定

有：

因为速度是位移随时间的变化

有：

将上列式子进行整个，可以得到某一时刻t与位置x的关系:

其中x₀和v₀分别表示初始位置和初始速度。加速度可以表示为势能对位置的导数：

所以，为了求取一个运动轨迹，必须知道系统中粒子的初始位置，初始速度分布和加速度，这些都可以从势能导数方程推导得到。比如，开始时间的粒子的位置和时间可以准确得到某一时刻t的位置和速度。初始位置可以系统实验结构数据得到，如Ｘ光衍射或者ＮＭＲ实验。系统初始速度分布主要来自计算机随机分布，同时根据系统温度进行修正以确保系统的整体动量为零：

其中，速度v_i是根据温度在一个波尔兹曼分布或者高斯分布中随机抽取，所以原子i可以得到在温度Ｔ在ｘ方向上的速度v_x。

反过来，系统温度可以从速度求得：

其中Ｎ表示系统中粒子数。

-------------------------------

算法综合
由于系统地势能是系统内所有原子的位置的函数，方程十分复杂而没有解析解，只能用数值方法求解。几种常用的数值算法有：

1、Verlet算法
2、Leap-frog算法
3、Velocity Verlet算法
4、Beeman's算法

在选择不用算法是，要注意算法是否能保持系统内部动量的恒定，并可以高效地解决问题，同时允许进行比较长的步长模拟。

-------------------------------

所有以上的算法都假定系统中原子的位置、速度和加速度可以利用泰勒张开近似求得：

其中，r表示位置，v表示速度（随时间变化的第一个因变量），a为加速度（随时间变化的第二个因变量）。

Verlet算法可以表达成以下形式：

将两式合并可得到：

Verlet算法利用在t和t - &delta t时刻的位置和加速度以计算在t + &delta t时刻的新位置。这个算法并没有计算原子的精确速度。这个算法主要优点有：i )可以直接求解；ii )计算中加饥结果少，可以节省内存。但是这个的缺点就是精确度比较小。

-------------------------------
Leap-frog算法

在这个算法中，先计算在t+1/2&delta t时刻的速度，然后再用其计算在t+&delta t时刻的位置。这样，速度根据速度进行计算，而位置又根据得到的速度进行计算，一步步跳跃向前。这个算法得优势是速度可以精确计算，但是缺点是不能在同一时刻求得速度和位置。t时刻的速度可以用下列关系大约求得：

-------------------------------
Velocity Verlet算法

这个算法在同一时刻求取位置、速度和加速度，但是精确度不高，主要表达为：

-------------------------------
Beeman算法

这个算法从Verlet算法演变而来，表达为：

这个算法比较精确的表达了速度，同时也确保能量守恒。但是由于算式复杂，计算量较大。
-------------------------------

分子动力学模拟概论

分子动力学概论

三、统计力学理论

分子动力学模拟可以得到原子显微水平的信息，包括原子位置和运动速度。要将这些微观信息转化为宏观水平信息，如压强、能量、热容量等等，需要用到统计力学理论。统计力学是分子动力学模拟生物系统研究的基本理论。以下讨论内容，主要讨论统计力学的最基本的概念，详细讨论请参考详细论述著作。

1、统计力学概论
在分子动力学模拟计算中，通常希望如计算一个分子系统的待定药物小分子的结合自由能或者构象转化的能量变化等得到这个分子系统的宏观性质。这种微观模拟以得到宏观性质的方法是以统计力学为基础的，因为统计力学用严格的数学方式表达了这种宏观状态与系统内部原子和分子运动之间的关系；分子动力学模拟方法就是根据这些严格的数学公式描述粒子的运动并求取这些数学公式的结果。利用分子动力学模拟，可以研究系统的热力学性质以及时间依赖性的动力学现象。
----------------------------------
统计力学参考书目
D. McQuarrie, Statistical Mechanics (Harper & Row, New York, 1976)
D. Chandler, Introduction to Modern Statistical Mechanics (Oxford University Press, New York, 1987)
R. E. Wilde and S. Singh, Statistical Mechanics, Fundamentals and Modern Applications (John Wiley & Sons, Inc, New York, 1998)
----------------------------------

统计力学是物理学的一个从分子角度研究一个宏观系统的分支。统计力学的研究目的主要是从构成系统的单独分子的性质出发，研究并预示系统的状态。这样的系统可以使液体分子系统，也可以是蛋白质与DNA复合物溶剂系统等等。为了把系统宏观性质与微观性质联系起来，必须引用系统随时间变化的统计平均量，这些主要的平均量的主要有：

热力学状态：一个系统的热力学状态是由一系列的参数组成，一般指系统的温度T，系统地压强P，和系统包含的粒子数量N三个参数。其他的系统参数可以由以上三个参量根据热力学系统状态方程和其他热力学基础方程推导得到。

微观（机械）状态：一个系统微观状态主要由系统内部各个原子的位置q，动量p决定，这些量可以组成一个多维空间即相空间的坐标。如果一个系统具有N个粒子，那么这个系统将有6N个维数。在相空间内的一个点，标记为&Gamma，描述系统的状态。一个系综就是相空间内这些满足可以描述系统特殊热力学状态的点的集合。分子动力学模拟可以利用时间为变量的方程产生相空间内一系列的点；这些点属于同一个系综，他们动量的变化的集合即表现为系统的变化。系综有多种，主要有以下几种。
----------------------------------
系综：一系列具有不同微观状态却表现出相同宏观热力学状态的系统的集合。主要系综类型有：
&oplus：微正则系综（NVE），该系综的系统的热力学状态表现为内部粒子数N恒定不变，体积不变，同时系统能量也不变。该系综主要见于孤立系统。
&oplus：正则系综（NVT），这类系统内部热力学状态表现为粒子数N恒定，同时系统的体积和温度都保持不变。
&oplus：等压等温系综（NPT），这类系统粒子数恒定，同时压强和温度保持不变。
&oplus：巨正则系综（&mu VT），这一类系统的化学势&mu保持不变，同时体积和温度保持不变。
----------------------------------

2、分子动力学模拟计算平均量

化学实验对象通常是一个可以观察其宏观状态的包含成千上万个原子或者分子的样品。在统计力学研究中，实验对象作表现的状态和现象对应于系综的平均值，这已经被实验所证明。一个系综的平均值是同时对系综所有系统同时取样得到。

在统计力学中，平均值定义为系综的平均值。系综平均值由

得到。其中有意思的是可以表达为一个以动量p和位置r为变量的方程。

系综的可能密度可以由

得到，其中H为哈密尔顿函数，T为温度，k_B为波尔兹曼常量，而Q为配分方程

以上的完整积分求解十分困难，因为必须求解系统的所有状态。在分子动力学模拟过程中，系统中的点随时间变化连续求解以计算系统状态的平均值，即必须求解系统的所有状态以得到某一个特定的热力学约束。

分子动力学模拟的另外一个方法是求解A随时间变化的平均值，方程表达为

其中&tua为模拟时间，M表示模拟的步数，A(pN,rN)表示A的瞬间值。

分子动力学模拟得到的是系综随时间变化的平均值，而实验观察到的是总体平均值，两者根据统计力学基础统计假设是相等的的，有

即：时间平均值=总体平均值。

这里可以这样理解这个基本概念：如果允许一个系统不确定地自由变化，那么系统将遍历所有可能的状态。所有分子动力学模拟的一个目的就是产生足够多的系统状态以确保等式成立。这样，实验要得到的结构、动力学和热力学信息就可以一举一定的计算能力得到，因为为了模拟一定的模拟时间，必须求解足够的相空间。

一些平均量的例子：

平均势能

其中M表示在模拟轨迹中构象数，Vi表示每一个构想的势能。

平均动能

其中M为模拟过程的构象数，N为系统中原子数，mi和vi分别表示原子i的质量和速度。

如上所述，一个分子动力学模拟必须要有足够长的模拟时间，才能得到足够多的模拟构象，使模拟更加接近现实。

Docking基本步骤

分子对接Docking基本步骤，使用UCSF分子对接程序Dock5.x或者6.x。

首先获得靶蛋白分子PDB文件，可以从PDB database网站上搜索。同时取得小分子的PDB或者mol2文件。Dock程序要求分子文件必须为mol2结构，因为mol2文件可以记录分子中各个原子的代电量等，包含信息比较多。准备分子对接文件时候，最后也使用UCSF提供的Chimera(http://www.cgl.ucsf.edu/chimera),好用免费，同时功能齐全。也可以使用VMD，但是要和tleap程序一起用。

首先对靶蛋白分子进行处理，去除靶蛋白PDB文件中包含的其他小分子和水分子，无关紧咬的金属离子等。以Chimera为例，可以使用Dock Prep功能一步搞定。然后利用Chimera给靶蛋白分子加上H原子，同时计算代电量。如果靶蛋白的分子量过大，再后续步骤中可能会出错，所以无关必要的要清除干净。保存靶蛋白为mol2格式。然后对要进行docking的小分子进行同样操作，加氢和加电，然后保存为Mol2格式。

准备完两个基本对接分子后，要准备靶蛋白分子的分子表面文件，需要使用dms程序，同时需要靶蛋白分子的PDB文件。用Chimera去除靶蛋白分子的氢原素，保存为pdb文件格式。然后利用dms生成靶蛋白分子表面文件。命令格式可以为：
-------------------------------------------------------------------
dms input_file [-a -d density -g file -i file -n -w radius -v] -o file
-a use all atoms, not just amino acids
-d change density of points
-g send messages to file
-i calculate only surface for specified atoms
-n calculate normals for surface points
-w change probe radius
-v verbose
-o specify output file name (required)
-------------------------------------------------------------------
例，如果靶蛋白无氢PDB文件为rec_noH.pdb，可以使用一下命令生成表面文件：dms rec_noH.pdb -n -w 1.4 -v -o rec.ms。其中参数都为默认值。rec。ms为靶蛋白分子表面文件，可以用Chimera监查。dms可以从一下地址下载：http://www.cgl.ucsf.edu/Overview/software.html#dms

下一步使用Dock程序包自带的sphgen程序生成靶蛋白的球面文件，即计算靶蛋白分子表面的性状各个不分的曲率半径。sphgen要求一个INSPH文件在目录下，格式如下：
-------------------------------------------------------------------
rec.ms #molecular surface file
R #sphere outside of surface (R) or inside surface (L)
X #specifies subset of surface points to be used (X=all points)
0.0 #prevents generation of large spheres with close surface contacts (default=0.0)
4.0 #maximum sphere radius in Angstroms (default=4.0)
1.4 #minimum sphere radius in Angstroms (default=radius of probe)
rec.sph #clustered spheres file
-------------------------------------------------------------------
其中各项为默认值，细节氢参考Dock文档。rec。ms为上一步生成的分子表面文件，rec。sph为球面文件。

下一步需要使用Dock自带的showspere程序显示有可能接合的位置，showspere可以是问答式，也可以使用一个配置文件以：showspere < showsphere.in 命令方式运行。showsphere。in 文件格式如下：
-------------------------------------------------------------------
rec.sph #sphere cluster file
1 #cluster number to process (<0 = all)
N #generate surface as well as pdb file
selected_cluster.pdb #name for output file
-------------------------------------------------------------------
运行后得到selected_cluster.pdb文件，即为有可能接合位置文件，可以使用Chimera查看。再用DDock自带的sphere_selector程序参考小分子mol2文件选择有可能接合的位置，命令如下。
sphere_selector rec.sph lig_charged.mol2 10.0
其参数为默认值。lig_charged。mol2为小分子带电荷文件。该程序会生成sphere_selected.sph文件，即为小分子有可能接合位置文件。

接下来使用Dock自带的showbox程序划出有可能接合的几何区域，dock将在这个Box里面进行，所以很重要。命令方式如下：
showbox < box.in
box.in文件格式为
-------------------------------------------------------------------
Y
5
sphere_selected.sph
1
rec_box.pdb
-------------------------------------------------------------------
命令得到rec_box。pdb文件，即为小分子与靶蛋白进行对接的区域。到此Docking的工作完成一半。剩下的就是Grid和docking。如果没有出错，后面做docking一般不会出现错误。

Dock自带的grid和dock程序的输入格式基本一样： grid -i grid.in -o grid.out/dock -i dock.in -o dock.out.
grid.in和dock.in为输入文件，参数较多，请参考dock手册。

AMBER分子动力学简例（三）

分子动力学（2）

水环境中的分子动力学模拟

溶剂环境中的分子动力学模拟分为以下四步进行：

1、溶剂环境能量最优化。这一步保持溶质（蛋白）不变，去除溶剂中能量不正常的范德华相   互作用。

2、整系统能量最优化。去除整个系统中能量不正常的相互作用。

3、有限制的分子动力学。保持蛋白质不动，溶解溶剂的不同层，同时逐渐将系统温度从0K提升到300K。

4、整系统分子动力学模拟。在一个大气压，300K的环境下整个系统分子动力学模拟。可以得到成果的分子动力学模拟。
#############################

1、溶剂环境能量最优化。

该步骤的配置文件min1.in如下：

---

oxytocin: initial minimisation solvent + ions
&cntrl
imin = 1,
maxcyc = 1000,
ncyc = 500,
ntb = 1,
ntr = 1,
cut = 10
/
Hold the protein fixed
500.0
RES 1 9
END
END

---

该过程保持肽链不动，其中500.0单位是kcal/mol，表示作用在肽链上使其不动的力。“RES 1 9”表示肽链残基数目，因为我们学习使用的oxytocin有9个残基。

模拟命令如下：

sander –O –i min1.in –o min1.out –p oxy.top –c oxy.crd –r oxy_min1.rst –ref oxy.crd &

2、整系统能量最优化。

配置文件min2.in如下：

---

oxytocin: initial minimisation whole system
&cntrl
imin = 1,
maxcyc = 2500,
ncyc = 1000,
ntb = 1,
ntr = 0,
cut = 10
/

---

命令如下：

sander –O –i min2.in –o min2.out –p oxy.top –c oxy_min1.rst –r oxy_min2.rst &

3、有限制的分子动力学。

第一步分子动力学保持蛋白分子位置不变，但是不是完全固定每个原子，同时缓解蛋白分子周围的水分子，是溶剂环境能量优化。在这个步骤中，我们将主要目的是对特定的原子使用作用力使其能量优化。我们要优化溶剂环境，至少需要10ps，我们将使用20ps用来优化我们上两步制作的分子系统的周期性边界的溶剂环境。

命令配置文件md1.in如下：

---

oxytocin: 20ps MD with res on protein
&cntrl
imin = 0,
irest = 0,
ntx = 1,
ntb = 1,
cut = 10,
ntr = 1,
ntc = 2,
ntf = 2,
tempi = 0.0,
temp0 = 300.0,
ntt = 3,
gamma_ln = 1.0,
nstlim = 10000, dt = 0.002,
ntpr = 100, ntwx = 500, ntwr = 1000
/
Keep protein fixed with weak restraints
10.0
RES 1 9
END
END

---

上述参数解释如下：
ntb = 1：表示分子动力学过程保持体积固定。
imin = 0：表示模拟过程为分子动力学，不是能量最优化。
 nstlim = #：＃表示计算的步数。
dt = 0.002：表示步长，单位为ps，0.002表示2fs。
temp0 = 300：表示最后系统到达并保持的温度，单位为K。
tempi = 100：系统开始时的温度。
gamma_ln = 1：表示当ntt＝3时的碰撞频率，单位为ps^-1（请参考AMBER手册）
ntt = 3：温度转变控制，3表示使用兰格氏动力学。
tautp = 0.1：热浴时间常数，缺省为1.0。小的时间常数可以得到较好的耦联。
vlimit = 20.0：保持分子动力学稳定性速度极限。20.0为缺省值，当动力学模拟中原子速度大于极限值时，程序将其速度降低到极限值以下。
comp = 44.6： 溶剂可压缩单位。
ntc = 2：Shake算法使用标志位。1表示不实用使用，2表示氢键将被计算，3表示所有键都将被计算在内。
tol = #.#####：坐标位置重新设置的几何位置相对容忍度。

我们将使用一个较小的作用力，10kcal/mol。在分子动力学中，当ntr＝1时，作用力只需要5-10kcal/mol（我们需要引用一个坐标文件做分子动力学过程的比较，我们需要使用"-ref"参数）。太大的作用力同时使用Shake算法和2fs步长将使整个系统变得不稳定，因为大的作用力使系统中的原子产生大频率的振动，模拟过程并步需要。

运行命令如下：

sander –O –i md1.in –o md1.out –p oxy.top –c oxy_min2.rst –r oxy_md1.rst –x oxy_md1.mdcrd –ref oxy_min2.rst –inf md1.info&

进行MD的运行时间一般较长，可以使用程序的并行版本提交集群计算。在主频位2.0GHz的P4单机上，大概需要一个钟头。可以随时查看md1.in文件的程序输出。

4、整系统分子动力学模拟。

这一步中，我们将进行整个系统的分子动力学模拟，而不对某些特定原子位置进行限制。因为知识一个小例子，我们将只进行250ps的MD计算。配置文件md2.in如下：

---

oxytocin: 250ps MD
&cntrl
imin = 0, irest = 1, ntx = 7,
ntb = 2, pres0 = 1.0, ntp = 1,
taup = 2.0,
cut = 10, ntr = 0,
ntc = 2, ntf = 2,
tempi = 300.0, temp0 = 300.0,
ntt = 3, gamma_ln = 1.0,
nstlim = 125000, dt = 0.002,
ntpr = 100, ntwx = 500, ntwr = 1000
/

---

上面一些参数解释如下:
ntb=2：表示分子动力学过程的压力常数。
ntp＝1：表示系统动力学过程各向同性。
taup = 2.0：压力缓解时间，单位为ps。
pres＝1：引用1个单位的压强。
使用以下命令进行MD：
sander –O –i md2.in –o md2.out –p oxy.top –c oxy_md1.rst – r oxy_md2.rst –x oxy_md2.mdcrd –ref oxy_md1.rst –inf md2.info &

模拟时间较长，在P4单机2.0GHz的上需要7.5个小时。
到此，模拟全部完成，接下来要对得到的数据进行分析。主要数据文件.out文件，包含系统能量、温度，压力等等；.mdcrd文件，是分子动力学轨迹文件，可以求系统蛋白的RMSD，回转半径等等。数据分析根据不同的研究目的不同而不同，我们将在后文中进行一些简单的分析。

AMBER分子动力学简例（二）

分子动力学（1）

真空模式

真空模式分子动力学模拟将使用NVT系宗分两步进行，即系统能量最优化和分子动力学过程。

1、系统能量最优化。我们将使用淬火能量最优化解除系统内原子之间的不正常相互作用，这些原子之间的高能量相互作用如果不消除，可能影响后续的分子动力学过程。因为动力学过程是能量梯度变化的，太高的能量壁垒可能让MD局限在某一个能量局部最小化位置中。

Sander程序是分子动力学模拟程序，它的主要功能是能量最优化，动力学模拟和NMR优化计算。我们必须给Sander一个运行的配置文件，使其按照我们的要求进行计算。能量最优化的配置文件如下：

"min_vac.in"

---

oxytocin: initial minimization prior to MD

&cntrl
imin = 1,
maxcyc = 500,
ncyc = 250,
ntb = 0,
igb = 0,
cut = 12
/

---

Sander程序的基本命令参数如下：

sander –O –i in –o out –p prmtop –c inpcrd –r restrt [-ref refc –x mdcrd –v mdvel –e mden –inf mdinfo]

所以我们的命令可以如下：

sander –O –i min_vac.in –o min_vac.out –p oxy_vac.top –c oxy_vac.crd –r oxy_vacmin.rst &

可以使用more命令或者tail命令查看min_vac.out的输出内容，那是能量最优化的记录。

2、分子动力学模拟。

Sander程序的配置文件为：

md_vac.in

---

oxytocin MD in-vacuo, 12 angstrom cut off, 250 ps
&cntrl
imin = 0, ntb = 0,
igb = 0, ntpr = 100, ntwx = 500,
ntt = 3, gamma_ln = 1.0,
tempi = 300.0, temp0 = 300.0,
nstlim = 125000, dt = 0.002,
cut = 12.0
/

---

命令为：

sander –O –i md_vac.in –o md_vac.out –p oxy_vac.top –c oxy_vacmin.rst –r oxy_vacmd.rst –x oxy_vacmd.mdcrd –ref oxy_vacmin.rst –inf mdvac.info &

MD的计算过程一般比较久，真空相对与溶剂中要快以下，250ps的模拟大概在一个主频为2。0GHz的Linux单机上运行5分钟。

以上配置文件中用到很多参数，这些参数这是sander程序参数的一小部分，以下将相应解释。如果要了解sander的其他参数，请阅读AMBER用户指南。

&ctrl和"/"：sander的参数一般要求出现在这两个标识符号之间，参数以及这两个标识符被称做控制模块。

cutoff：以纳米为单位的截矩。即超出截矩范围的非键连接相互作用将不计。

ntr：原子位置能量抑制位，1表示抑制，0表示不抑制。

imin：能量最优化标志位。1表示sander将进行能量最优化，0表示让sander进行分子动力学模拟。

macyc：能量最优化次数。

ncyc：便是经过多少次能量优化以后，能量优化从淬火过程变为梯度变化过程。

ntmin：能量优化方法标志位。0表示前10个能量最优化为淬火过程，然后进行梯度能量优化；1表示ncyc次淬火过程，然后进行梯度能量优化，为默认值；2表示只进行淬火过程。

dx0：表示启动模拟步长。

dxm：最大优化步数。

drms：梯度能量优化标准，默认值为1.0E-4 kcal/mol.A。

更多参数将在后文中解释。

AMBER分子动力学简例（一）

概述

以下是使用AMBER包的简单教程，希望对开始学习分子动力学的同学有用处。申明一下，以下教程原版来自网上，是最最基本的教程，同时也非常实用，有非常好的借鉴意义。

AMBER分子动力学程序包是加州圣弗兰西斯科大学（University of California San Francisco,UCSF)的Peter A Kollman和其同事编的，程序很全，现在已经发展到版本9.0。AMBER功能涵盖种类非常多的生物分子，包括蛋白、核算以及药物小分子。软件详细情况请浏览http://amber.scripps.edu.

以下是AMBER软件包中四个主要的大程序：

---

Leap：用于准备分子系统坐标和参数文件，有两个程序：

           xleap：X－windows版本的leap，带GUI图形界面。

           tleap：文本界面的Leap。

Antechamber：用于生成少见小分子力学参数文件的。有的时候一些小分子Leap程序不认识，

                        需要加载其力学参数，这些力学参数文件就要antechamber生成。

Sander：MD数据产生程序，即MD模拟程序，被称做AMBER的大脑程序。

Ptraj：MD模拟轨迹分析程序。

---

学习项目

本教程研究的题目是脑下垂体荷尔蒙之一的oxytocin，需要X光衍射晶体结构文件1NPO.PDB。该文件包含了该荷尔蒙和其运载蛋白的复合物，可以从蛋白数据库下载。

PDB文件是不包含氢原子的，Leap程序会自动的加上PDB文件缺少的东西。当第一次使用PDB文件的时候，要十分留意文件包含的信息，所以PDB文件缺少的残疾、侧链或者添加的变异都在这个地方记录。可以用文本阅读程序阅读PDB文件头部的信息，即以REMARKS开头的信息文本行。PDB一个重要的信息是SSBOND记录，该记录说明结构中二硫键的位置，这样的信息在使用Leap程序建立分子结构的时候需要。在本教程中，我们将比较oxytocin在真空和溶液中分子动力学的差异，如果没有二硫键，将影响整个结果。整个过程在Linux系统下完成，如对Linux不熟悉，请翻阅有关Linux书籍。

建立项目目录，并进入该目录：

mkdir  myproj

cd myproj

下载1NPO.PDB文件到该目录下，使用文本阅读器阅读该PDB文件。从在文件的开头部分可以得到该文件是一个二聚物的PDB文件，删除文件中A、C、和D链对应的文本行，剩下的就是oxytocin的晶体机构了，保存为oxyt.pdb文件。

使用Swiss PDB viewer分析以下oxyt.pdb文件，可以得知该pdb文件缺少了一个侧链。如果使用Deep View，它会自动给文件加上侧链。重新保存pdb文件为oxyt.pdb文件。

xleap和tleap程序功能是一样的，都是准备分子结构的坐标文件和拓扑文件。xleap启动一个X界面，比较慢。tleap纯文本，要快一些，我们选择tleap：

tleap  -s  -f    leaprc.ff03

其中leaprc.ff03是AMBER的2003力场文件。力场是一个很重要的文件，定义了分子、原子和残疾等的信息，请查阅有关资料。

oxy=loadpdb oxyt.pdb

该命令读入oxyt.pdb文件到oxy变量中。

bond  oxy.1.SG oxy.6.SG

连接oxy变量的第一和第六个残疾的SG原子，即是连接二硫键。可以使用check命令检查oxy是否完好，没有特殊情况的话，最后应该是“OK”，表示分子系统状态是好的，可以保存。

check oxy

接下来就是保存分子系统了：

saveamberparm  oxy  oxy_vac.top  oxy_vac.crd

其中oxy_vac.top和oxy_vac.crd分别是分子系统真空状态下的拓扑和坐标文件。

接下来要给分子系统添加水环境，

solvateoct  oxy  TIP3PBOX   9.0 

该命令让Leap程序使用TIP3PBOX水模型,水环境距离分子为9纳米。也可以用solvatebox命令，但是solvatebox添加的水环境是一个立方体，不是八面体。一般要求水表面到蛋白距离要达到8.5纳米，避免MD过程中蛋白冲出水环境，但是水环境越大计算时间将越长。如果MD过程内含PME，那么水箱长度要大于2倍截矩（cutoff），截矩在MD配置文件中定义。加溶剂之后，要计算系统是否为电中性，使用命令：

charge oxy

如果现实不是电中性，要使用addions添加反性电荷，如Cl－或者Na＋等。

最后保存溶剂环境的系统拓扑结构和坐标文件，并退出Leap程序：

saveamberparm  oxy   oxy.top  oxy.crd

quit

也可以把命令集中到一个，让tleap读取。如将上面的命令集中到以下文件中

"oxy.leaprc"

---

source  leaprc.ff03

oxy=loadpdb oxyt.pdb

bond oxy.1.SG  oxy.6.SG

check oxy

saveamberparm oxy oxy_vac.top oxy_vac.crd

solvateoct oxy

saveamberparm oxy oxy.top oxy.crd

quit

---

将以上两横线之间的命令保存为"oxy.leaprc"，然后使用以下命令一步搞定：

tleap  -s  -f   oxy.leaprc

到此，分子系统准备已经完成，目录中会多了一个leap.log文件，这是Leap程序的记录文件，它包含了leap程序所做的一切，包括自动的和手动命令。

现在已经有了可以进行分子动力学的基本文件了，再编写sander的控制文件，就可以进行模拟了。这些将在后文中介绍。

NAMD分子动力学模拟软件范例

使用NAMD运行分子动力学模拟，需要一下四类文件。
• CHARMM力场文件，可以是CHARMM文件格式或者X－PLOR文件格式。（其他立场也可以，如AMBER或者GROMACS等）
• PSF文件格式的分子结构，可以使用NAMD程序包自带的psfgen程序生成。
• 包含分子系统最初结构的PDB文件，可以从PDB数据库下载。
• NAMD主程序namd2的运行配置文件。

以下是NAMD 2.6用户指南文档中一个简单的例子。

如上所述，要运行这个例子，需要以下内容：
• psfgen程序
• CHARMM力场参数文件 top_all22_prot.inp 和 par_all22_prot.inp（http://www.pharmacy.umaryland.edu/faculty/amackere/force fields.htm）
• BPTI 例子的PDB 文件 6PTI.pdb （http://www.pdb.org/）
•NAMD运行配置文件（后文有述例）

首先，根据PDB文件构建分子结构，生成bpti.psf和bpti.pdb两个文件。以下为在Linux命令

1.准备基本PDB文件(6PTI.pdb已经在当前目录下)
mkdir -p output #建立工作目录
rm -f output/6PTI_protein.pdb output/6PTI_water.pdb #清空工作目录下的文件(多余的)
grep -v ’^HETATM’ 6PTI.pdb > output/6PTI_protein.pdb #从原始PDB文件中提取蛋白分子
grep ’HOH’ 6PTI.pdb > output/6PTI_water.pdb #从原始PDB文件中提取水分子

2.启动psfgen交互程序
psfgen #启动psfgen交互程序
topology top_all22_prot.inp(绝对或相对文件路径) #读入CHARMM立场文件

segment BPTI {
pdb output/6PTI_protein.pdb
} #读入蛋白分子

patch DISU BPTI:5 BPTI:55 #连接二硫键
patch DISU BPTI:14 BPTI:38
patch DISU BPTI:30 BPTI:51

pdbalias atom ILE CD1 CD #纠正ILE残基中原子命名
coordpdb output/6PTI_protein.pdb BPTI #读入原始PDB文件以纠正分子各原子坐标

# 以下是建立分子系统中水分子分子结构

pdbalias residue HOH TIP3 ; # 使用TIP3水分子模型

segment SOLV { #读入水分子
auto none
pdb output/6PTI_water.pdb
}

pdbalias atom HOH O OH2 # 标识水分子的原子
coordpdb output/6PTI_water.pdb SOLV #读入原始PDB纠正各原子坐标

guesscoord #psfgen自动建立分子系统中各原子坐标

writepsf output/bpti.psf #保存分子系统的psf文件
writepdb output/bpti.pdb #保存分子系统的PDB文件
exit #退出psfgen交互程序

到此基本的分子文件准备完毕,可以得到 bpti.psf 和 bpti.pdb两个文件.为了运行namd,需要建立一个namd的配置文件,一下为范例.

###################begin of namd.in########################
# NAMD configuration file for BPTI
# molecular system
structure output/bpti.psf #指定psf文件路径
# force field
paratypecharmm on
parameters par_all22_prot.inp #指定力场文件路径
exclude scaled1-4
1-4scaling 1.0
# approximations
switching on
switchdist 8
cutoff 12
pairlistdist 13.5
margin 0
stepspercycle 20
#integrator
timestep 1.0
#output
outputenergies 10
outputtiming 100
binaryoutput no
# molecular system
coordinates output/bpti.pdb #指定pdb文件路径
#output
outputname output/bpti
dcdfreq 1000
#protocol
temperature 0
reassignFreq 1000
reassignTemp 25
reassignIncr 25
reassignHold 300
#script
minimize 1000
run 20000

#####################end of namd.in#####################
基本参数比较多,请参考namd2.6的用户指南.

到此,已经可以运行namd进行动力学模拟啦,命令如下:

[sen@localhost]$ nohup namd2 namd.in >namd.out &

其中,namd.out为输出文件,包含能量等.

如更多详细,请仔细阅读namd 2.6 用户指南.

NAMD 简介

NAMD Molecular Dynamics Software(NAMD分子动力学计算软件)

自由,开源

开发团队:Illinois大学,Illinois 贝克曼研究所,理论生物物理小组.

NAMD 官方主页: http://www.ks.uiuc.edu/Research/namd

公布论文:

Laxmikant Kale,Robert Skeel,Milind Bhandarkar,Robert Brunner,Attila Gursoy,Neal Krawetz,James Ohillips,Aritomo Shinozaki,Krishnan Varadarajan,and Klaus Schulten.

NAMD2: Greater scalability for parallel molecular dynamics.

J. Comp. Phys.,151:283-312,1999.

NKstars的namd编译

下载NAMD,网址:http://www.ks.uiuc.edu/Research/namd/
包括NAMD源文件,还有tcl,fftw,plugins模块.
#选tcl,fftw模块的时候要选redhat9.0的,南开之星远行的就是这个.

解压NAMD
tar -xzvf NAMD*

在NAMD的目录下又一个charm的tar包,也解压
tar -xvf charm*

#编译namd要用charm++编译,所以先编译charm,然后用charm编译namd.

转换到charm目录,运行目录下build脚本
如下:
./build charm++ mpi-linux gm icc --incdir=/usr/local/mpich/1.2.6..13a/gm-2.1.3aa2nks3/smp/intel32/ssh/include --libdir=/usr/local/mpich/1.2.6..13a/gm-2.1.3aa2nks3/smp/intel32/ssh/lib

#该脚本编译charm++,脚本编译的是mpi-linux架构的charm++,所以编译的时候加入了mpich的头文件路径跟裤路径.

编译完了要修改charm目录下的./src/arch/mpi-linux/conv-mach.sh文件,为CMK_LIBS变量中的"-lmpich"指明路径,又mpich要用到gm库,所以也把gm库的引用在这个变量中申明,如下:
CMK_LIBS='-lckqt -L/usr/local/mpich/1.2.6..13a/gm-2.1.3aa2nks3/smp/intel32/ssh/lib -lpmpich -L/opt/xcat/gm/lib/ -lgm'

#这个文件是charm++编译其他软件的时候的引用文件,如果你想让charm++引用其他的编译器,可以在这个文件中修改.

charm++编译到此结束,可以到charm目录下的tests子目录下测试charm++,如果出现问题,病不能解决,建议利用高一点版本的charm,比如charm5.9.

"cd .."到namd目录,建立tcl,fftw,跟plugins目录,把原先下载的tcl,fftw,plugins模块复制到对应目录下.

修改arch目录下的Linux-i686.fftw,Linux-i686.tcl和All-Unix.plugins文件,修改这些文件中所指的对应库文件的路径.可以查看目录下的Linux-i686-*.arch文件,该文件说明namd的编译要用到的charm++构架,也就是上文所述编译的charm.

运行namd目录下config脚本,引入这些模块,如下:
./config tcl fftw plugins Linux-i686-MPI-icc

#改脚本是一个链接程序,把要建立的namd的所有源文件链接到Linux-i686-MPI-icc目录下.

"cd Linux-i686-MPI-icc"到该目录下,执行make,可以生成namd2程序.

#保佑你不要出错,阿门!!