GMX步骤解释2

1、科研笔记1：Gromacs安装以前跑MD都是用NAMD软件包，今天开始正式使用Gromacs软件包。在安装过程中看了网上很多的安装方法，其实大部分都已经过时，这里专述一文，帮助大家入门。第一步：下载最新版本的Gromacs软件包，比如2014.4.3最新正式版本是gromacs-4.6.5.tar.gz。注意：由于4.6.x及其以上版本和4.5.x及其以下版本使用的安装方式完全不同，4.6.x使用CMake（Cross Platform Make）编译，方便很多。参考网站：/Downloads第二步：进行安装按照以下指令：tar xfz gromacs

2、-4.6.5.tar.gzcd gromacs-4.6.5mkdir buildcd buildcmake . -DGMX _BUILD_OWN_FFTW=ONmakesudo make install参考资料：/Documentation/Installation_Instructions第三步：添加环境变量在我的文件夹里，比如说home/chenqu（举个例子）终端控制台输入：cd ls -a （列举当前目录下所有文件）就能看到文件.bashrc了。vi .bashrc在最后一行添上：export PATH=/usr/local/grom

3、acs/bin:$PATH保存退出wq运行：source .bashrc让PATH路径生效。科研笔记2：如何添加自己的力场参数1GROMACS实例分析1：CO2+H2O体系大部分GROMACS在线指南都是生物体系，但是对于简单的小体系提到的并不多。尽管原理一样，我们还是希望能了解如何自己写简单气体的力场。这里我做一个CO2和H2O混合体系的例子来说明二元混合体系的构建方法，做一个入门指南。目标：计算CO2在水中的扩散系数。第一步：使用Avogadro软件构建一个CO2.pdb。如下所示。如果不会使用Avogadro软件的，可以研究教程http:/avogadro.cc/wiki/Categor

4、y:Tutorials，这里不再赘述。创建完之后如图1所示。                                             

5、  图1：CO2.pdb的脚本注意：不要随意去修改pdb的格式，他的格式每列都是有特殊含义的。文本.pdb的坐标单位是angstrom，而文本.gro的坐标默认单位是nm，千万牢记！第二步：创建CO2的力场。其实这部分是初学者最需要学习的。我喜欢建立一个CO2.itp，这样的好处是：清爽简洁。我选择OPLS-AA力场作为主要力场。由于我们最终的.top文件里面需要#include “forcefield.itp”，而forcefield.itp里面又有#include "ffnonbonded.itp", #include "ffbonded.

6、itp"以及#include "gbsa.itp"，那么在这些力场中，都需要注明自己创建的CO2的各种参数（这些文件都在gromacs-4.6.5/share/top/oplsaa.ff）。首先把ffnonbonded.itp拷贝到自己的文件夹（原gromacs文件夹不要随便改变），然后在最后添加上 atomtypes ，我把CO2的碳定义为CO，而把氧定义为OC。貌似ffbonded.itp和gbsa.itp是不需要修改的，只要写在CO2.itp中就行。修改完如图2所示。       

7、0;                                                 

8、0;              图2：ffnonbonded.itp的部分脚本其次开始撰写CO2.itp。建议使用GROMACS教程4.5.4版本p.131的模板进行写作。CO2的参数需要根据不同的目的而不同。本例子是为了做扩散系数，参考文献如下：Harris, Jonathan G., and Kwong H. Yung. "Carbon dioxide's liquid-vapor coexistence curve and crit

9、ical properties as predicted by a simple molecular model." The Journal of Physical Chemistry 99.31 (1995): 12021-12024. 写完之后的文件请见图3。注意：由于CO2的力场中bonds是固定的，不需要键伸缩力常数。所以我不采用 bonds ，而采用 constraints 。把CO2.itp、forcefield.itp、ffnonbonded.itp、ffbonded.itp、gbsa.itp和CO2.pdb放在一个新文件夹里面。  

10、60;                                                 

11、60;                                          图3：CO2.itp的全部脚本第三步：建立程序的拓扑文件。把所有力场都#include到一

12、个拓扑文件.top里面即可。如下所示，建议大家最好了解下每一个itp的文本到底是什么。有什么含义，这些都是模拟中最重要的东西。                                       

13、60;          图4：topol_spce.top的脚本第四步：添加水盒子。先通过CO2的溶解度，计算一个CO2该溶解在多少个水分子中。其次计算水盒子大小，使得密度尽量控制在0.9-1.1 g/ml。操作指令为：genbox -cp CO2.gro -cs spc216.gro -box 3.83 3.83 4.04 -p topol_spce.top -maxsol 1879其中最后一项控制最多添加水分子的数量，可以根据实际情况添加。添加完别忘记修改里面的原子类型名字和残基类型名字，

14、和自己的拓扑力场参数匹配。        图5：CO2-spce.gro部分脚本第五步：进行能量最小化（Energy Minimization）。在网上搜索，下载一个em.mdp文件，见图5所示。这里从Justin Lemkul的个人教程中下载。下载地址为：/Pages/Personal/justin/gmx-tutorials/lysozyme/Files/ions.这时候em.mdp（参数文件）、CO2-spce.gro（坐标文件）以及topol

15、_spce.top（拓扑文件）三个结合生成一个二进制文件em.tpr。粗浅来看，这是我认为Gromacs比NAMD厉害的地方。运算由于读取二进制文件，所以速度极快。然后就开始运算了！指令是：grompp -f em.mdp -p topol_spce.top -c CO2-spce.gro -o em.tpr（合成一个二进制文件）mdrun -v -deffnm em（开始运算）这时你能发现，程序开始运行了。em.edr是二进制能量文件；em.gro能量最小化后的构型；em.trr二进制轨迹文件；em.log 整个最小化的进程日志。    图6

16、：em.mdp的全部脚本第六步：计算势能变化曲线。能量的信息是通过二进制码存储在em.edr文件中的。可以通过g_energy来计算出势能变化。指令是：g_energy -f em.edr -o potential.xvg根据屏幕提示然后选择，可以画出势能曲线。                            &

17、#160;                               图7：能量最小化势能曲线第七步：进行NPT平衡模拟（Equilibration）。和第五步一样，我们同时也需要npt.mdp（参数文件）、em.gro（坐标文件）以及topol-spce.top（拓扑文件）。最后两个已

18、经有了，我们还是需要从网上下载一个npt.mdp。根据你需要的条件参数进行修改。文件下载地址：/Pages/Personal/justin/gmx-tutorials/lysozyme/Files/npt.mdp.指令照旧：grompp -f npt.mdp -p topol_spce.top -c em.gro -o npt.tprmdrun -v -deffnm npt     图8：npt.mdp的全部脚本第八步：进行模拟运算（Production Run）。和第五步第七部

19、原理一样，我们已经有了npt.gro以及topol-spce.top，所以我们只需要nvt.mdp就可以了。文件下载地址：/Pages/Personal/justin/gmx-tutorials/lysozyme/Files/nvt.mdp.指令照旧：grompp -f nvt.mdp -p topol_spce.top -c npt.gro -o nvt.tprmdrun -v -deffnm nvt第九步：分析数据。我们需要计算CO2的扩散系数。首先，我们需要将周期性边界条件进行unfolding，因为粒子很有可能会跑到盒子

20、外面去，所以运行路程需要累加。指令为：trjconv -s nvt.tpr -f nvt.xtc -o nvt_noPBC.xtc -pbc mol -ur compact其次，做一个index文件。指令是： make_ndx -f nvt.gro -o index.ndx最后，可以利用函数g_msd来进行MSD对time的作图。指令为：g_msd -f nvt_noPBC.xtc -s nvt.tpr -n -o msd.xvg也可以顺便计算扩散系数，具体的建议直接上网去查找g_msd等其他function的使用方法，或者甚至仔细研究里面的程序，gmx功能相当强大，可以计算非常多的

21、东西。因此我使用过一次之后非常着迷。 参考资料：1.       /Pages/Personal/justin/gmx-tutorials/lysozyme/01_pdb2gmx.html2.       http:/verahill.blogspot.nl/2012/01/gromacs-carbon-dioxide-in-water-example.html3. &#

22、160;     http:/www3.mpibpc.mpg.de/groups/grubmueller/start/compbio/p2/科研笔记3：Gromacs对于MD的问题综述问题1：在npt系综运算完后，一定要进行密度的rescaling。在进行气体扩散的初始构型搭建过程中，最容易被人忽略的一点就是当npt系综将体系进行平衡之后，我们就容易拿出最后一帧的数据来进行运算，实际上这样并不合理。npt的最终目的就是让体系更接近于实验值，所以我们也要做准确的模型。参加科研笔记2，在进行完npt模拟后，需要按照如下步骤：先利用指令：g_ener

23、gy -f npt.edr -o density.xvg了解密度的波动，然后选择合适的计算范围。比如CO2在H2O的npt后得到密度波动见图1，我们发现40ps后密度相对比较稳定，于是从40ps到100ps做个平均值。然后利用指令：g_energy -f npt.edr -b 40 -e 100 -o density2.xvg，屏幕中自动会帮你计算平均值，比如密度是986.168 g/L。最后利用指令：editconf -f npt.gro -density 986.168 -o npt2.gro调整密度。npt2.gro就成了最终的构型。   

24、0;     图1 npt系综密度对时间的关系 问题2：遇到一个极其恼人的fatal error叫blowing up该怎么办？如果你遇到一个极其郁闷的致命错误如下：Fatal error:1 particles communicated to PME node 2 are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension x.注：下划线部

25、分，可能会被替换成同类数字或者字母。那么恭喜你，这个问题和其他的fatal error不一样，这个问题非常难解决。论坛上有人建议换版本，显然这不是我们想要的。有人上报为一个bug，这也和我们暂且无关。我的经验如下：首先，非常仔细检查自己的初始构型和itp文件看看有没有什么不妥的。笔者就曾经由于一个二面角的数字位置颠倒而导致浪费了一天时间找问题。千万仔细！其次，问题可能就出在你的热浴选择上，nose hoover是个强coupling的热浴，很容易造成这种类型的blowing up，建议换成berendsen或者v-rescale这类弱coupling热浴。berendsen是不准确的方式，一般

26、不推荐使用，笔者推荐使用v-rescale（比berendsen增加一个随机项）。第三，和你体系未达到平衡有关，建议再次平衡体系，在尝试。笔者有这样的经历，明明两个体系都平衡了，但是换了一个构型就可能把程序跑通了。第四，这个最有意思，只能呵呵。核用的越少，报错的可能性越小。笔者尝试过同样的输入文件，采用mdrun -nt 2 -v -deffnm nvt和mdrun -nt 12 -v -deffnm nvt得到了完全不同的结果，前者速度慢，但是程序容易跑通，后者速度快，但是经常报错。总结起来讲，这个问题非常麻烦，建议大家按照我的方案，反复尝试，也许会有收获。 问题3：关于Nose-

27、hoover的热浴中tau_t的选择问题。Berendsen热浴和v-rescale热浴，time constant一般选择0.1 ps；Langevin热浴也有类似的参数，有个damping coefficient，我们经过尝试在扩散小体系中用1 ps-1，在蛋白质大体系中可以用5 ps-1。那么Nose-hoover热浴呢？以下是我的尝试。我的体系很小1的CO2在500个水中，我选择了0.2 ps, 0.5ps和1ps做了三组平行实验。结果见图2。由于第二组0.5似乎是个相对较好的time constant。         图2 nvt系综nose-hoover控温tau选择依据科研笔记4：如何用VMD将轨迹文件制作动画 很多同学想从dcd(NAMD)或者trr(gromac)文件提取一段轨迹文件做成动画。我们接下来看看该怎么做。这里拿Gromacs举例，NAMD也差不多。第一步： 提取轨迹。往往轨迹文件很大，建议先提取其中的一段，或者把步长放大，提取完整的也行。利用指令：trjconv -f nvt.trr -b 0 -e 500 -o nvt_500.trr