gromacs基础知识

admin管理员组

文章数量:1532440

2023年12月24日发(作者：)

（18）gromacs基础知识

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的执行格式如下（假设你的分子坐标文件为）：

pdb2gmx -f -o -p -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模拟的要求高，一般都要有点水。为分子系统添加水环境和离子环境需要其他命令，要慢慢来。

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 -o

或者

editconf -f -o

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

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

-h "查看，但是制作一个重启tpr文件的参数和格式一般如下：

tpbconv -s -f -e -o

其中为原来的tpr文件，为双精度坐标文件（不要用xtc文件，因为精度不够），

为系统能量输出文件，是重启模拟文件。以上的命令得到的是在计算突然终止前一个系统构象的信息。也可以在命令中加上一个"-time "参数来指定从那一个时间重新开始，如一下指定从一纳秒处重新开始模拟：

tpbconv -s -f -e -time 1000 -o

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

tpbconv -s -f -e -extend 1000 -o

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

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

trjcat -f .... -o traj_

eneconv -f ... -o ener_

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

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

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

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

通常用到的genbox参数有一下几个：

-cp ：带盒子参数的分子坐标文件，也就是editconf的输出文件；

-cs ：添加的水分子模型，如spc216、spce、tip3p、tip4p等，关于各个模型的区别，请参考scholar google；

-o ：输出坐标文件，就是添加水分子之后的分子坐标文件，默认是.gro文件，但是也可以输出其他文件格式，如pdb；

-p ：系统拓扑文件，genbox会往里面写入添加水分子的个数，这个不要忘记，不然在进行下一步计算时，会出现坐标文件和拓扑文件原子数不一致的错误；

-ci ：索引文件，genbox可以为分子特定部位添加水环境，这样可以减少原子数，只有研究的分子部位添加水环境，节省计算时间。

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

使用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工具

(20)转自 sobereva 大师姐的blog，

(21)/sobereva/blog/item/

(22)这么好的东西sob 竟然不弄到这里来，看来她不管我们了，。

2009-06-20 19:10

本文主要讨论隐式溶剂模型和GB模型（广义Born模型）的各个方面。

隐式水模型优势主要有下：

减少体系的粒子数从而降低计算开销。

没有溶剂摩擦效应，可以加快构象变化。

不显著、具体地表现溶剂运动而直接表现平均效果，便于计算自由能，不需要动力学很长时间采样。

不需要显式水的预平衡过程。

适用于常PH模拟、REMD等方法。

适合计算能量面，结果比较光滑，真实的水会带来噪音，溶剂的运动、碰撞使主体分子产生许多能量局部极小点。等等......

总体来说，隐式水模型相对于显式水模型，是一个离散化到连续化的过程。溶剂模型的近似关系如下（越左边越真实，越右边近似程度越大）：

显式溶剂模型-->隐式溶剂模型-->明确拆分为静电、非极性贡献-->PB-->GB

分子在溶剂环境下的总能量Etot=溶剂化能deltaG(solv)+溶质在真空条件下的能量Evac。（对于极化力场，溶剂对溶质的极化作用造成的能变也算进在deltaG(solv)里。）

deltaG(solv)=deltaG(elec)+deltaG(nonpolar)=deltaG(elec)+deltaG(vdw)+deltaG(cavity)

注：有时明确包含氢键项。

其中deltaG(elec)就是溶剂环境下的静电作用能与真空下的静电作用能之差，主要包括溶剂与溶质间的静电作用，以及溶剂对溶质原子间静电作用的屏蔽效果。

其中deltaG(nonpolar)就是除此以外，即体系各原子皆无电荷时，溶剂环境与真空环境下能量之差。它又分为deltaG(vdw)和 deltaG(cavity)。deltaG(vdw)是溶质与溶剂的范德华作用，由于范德华作用随距离衰减快，所以只考虑第一溶剂化层的溶剂，此项造成能量降低。deltaG(cavity)表现的是溶质反抗溶剂压力形成洞穴的做功以及溶质周围溶剂（第一溶剂化层为主）的重构造成的熵罚，此项造成能量升高。

因为deltaG(nonpolar)的两项都主要涉及第一溶剂化层，它正比于SASA，故最简单、常用的表达是a*SASA+b，而比例系数a由拟合实验值得到，b一般为0。另有其它方案计算，比如包括溶质体积、溶质直径等信息在内的方程。也有比如deltaG(vdw)通过SAS的积分得到，而 deltaG(cavity)仍用正比于SASA的方案。在GB模型下的MD模拟中，对于原始状态的生物大分子，由于SASA变化很小，原子因非极性作用导致的受力往往可以忽略，若需计算的话则是求它对坐标的导数。由于a*SASA+b的a为正值（一般为7.2cal/mol/Å^2），可见

deltaG(nonpolar)的效果是减小溶质的SASA，此项体现了疏水效应。

deltaG(elec)项一般通过GB或PB方法计算，GB速度快但一般准确度低于PB。目前常用的GB出（原文/filebox/down/fc/90f984299e8583a61e0f33942bd13260），依据是库仑模型计算方程由Still于1990年提

定律和born方程，推导出另一种形式。当粒子相离较远时，方程等价于库仑定律+born方程能量，而相距为0时等价于born方程，而当距离较近时等价于Onsager反应场能量。方程中还可以再引入单价粒子的屏蔽校正项。

GB方程中重要参数是有效born半径。某个原子电荷与溶剂作用对deltaG(elec)的贡献实际上是把分子中其它原子电荷去掉后，溶剂化状态与真空状态能量差。如果用的有效born半径以GB方程算的结果与之一样，就是正确的born半径。或者粗略来说：某原子用有效born半径算得的GB式中的“自身项”的能量应当与真实的分子中此原子与溶剂的静电相互作用能一致。这样，有了正确的born半径，就可以用GB式得到合理的 deltaG(elec)。有效born半径对计算deltaG(elec)有重要影响，模拟过程中分子构象不断改变，依赖于它的有效born半径也不断改变，需要每一步重新计算，若构象变化不明显可每隔一定步数重新计算，这是GB计算量中的主要部分之一。计算有效born半径在不同程序中有不同方法。一般使用CFA（库仑场近似），有效born半径就可以对溶质分子表面内的空间进行积分获得。而分子内空间的范围不便于描述，简单的方法是直接将原子VDW球内区域叠加作为分子内空间，称GB-HCT，但这样就忽略了原子VDW空间之间的缝隙，可以乘以4/3来近似校正。也有人根据原子被埋程度的概念提出了基于经验的简单、快速的函数，其中引入可调参数，在amber里称为GB-OBC方法，结果明显优于GB-HCT，应注意可调参数是通过预先优化得到的，面向不同体系模拟应使用不同参数。amber中还支持更新的GBn方法，结果优于GB-OBC，适合蛋白而不适于核酸体系。

GB方程的优点是形式简单，计算快速，结果较准确，而且有简单的解析导数，可直接计算GB环境下原子在MD过程中的受力。GB可以结合分子力场及量子力学模拟各种分子体系的溶剂化效果。可以用在类似MM/PBSA的结合自由能计算的静电项中，即MM/GBSA。一些模拟方法则只能使用GB，比如 amber中可以实现的结合MC的常PH模拟，原理是根据Metropolis判据在MD过程中动态决定是否某点被质子化或去质子化。对于REMD，因为随着体系粒子数的增加，需要的副本数目飙升（否则能量重叠较差，交换概率低而起不到效果），所以总是结合GB方法来显著降低粒子数量。对于不大的体系，GB比显式溶剂模型有更快的速度。

但GB方法也有缺点。与显式水模型仍有一定差距，将水分子进行了“连续化”的近似，故不能表现与溶剂的强相互作用、特殊相互作用（如水桥）。对于膜体系，溶剂、溶质、膜都有着不同的介电常数，这样复杂的情况介电常数环境下GB也难以表达。对于电荷较多的核酸体系，也不能较好表达多价抗衡粒子与它的相互作用。一些研究也表明用GB模拟多肽会产生与显式溶剂模型不同的构象，与实验结果有一定偏差。计算deltaG(elec)项时，相对于原理严格的PB 也有差距（即便使用最佳的有效born半径），但是好处是GB计算速度快得多。不过PB的结果亦不可作为绝对的金标准，因为一些误差在PB引入的近似中就已经出现了。在计算速度上，显式溶剂体系计算量是O(N^2)，但使用Ewald等方法后，可降至Nlog(N)以下，而GB模型如果在cutoff上不做近似（往往取无限大），对于大体系反倒更慢。

除GB、PB以外也有其它一些方案，最简单的隐式模型是介电常数依赖于作用距离的方法，其中库仑作用的介电常数等于粒子间距离（即静电作用能1/r 变为了1/(r^2)），以体现溶剂的屏蔽效果，此方法精度显然低于GB。近来还有一些新方法，如AGBNP（Analytical Generalized

Born plus Nonpolar），在GB基础上引入另外形式的deltaG(nonpolar)项。ALPB（analytical

linearized Poisson-Boltzmann），其方程类似GB（在无限介电常数下回归为GB方程），故有着基本一致的计算速度，但引入了有效分子静电尺寸参数，模拟水比GB有着更好的精度，结果更接近于PB。这些新模型理论上有着更好的结果，但仍需广泛地应用于模拟来检验。还有隐式与显式溶剂模型相结合的方法，也就是在主体分子外面包一层溶剂分子，往往以弹簧势限制住溶剂避免跑走，而更外面则用隐式溶剂模型，结合了两种方法的一些优点。

目前Amber对隐式溶剂模型支持较好，支持GB、ALPB，也有内建的PB模块pbsa（老版本挂的是第三方的delphi）。而Gromacs在最新版本4.0.5中尚不支持GB/PB，需要外接第三方软件如APBS，但在接下来的版本中即将加入GB。

对于GB模型，适用领域与不适用领域的界限尚为模糊，一般来说，对于必须用GB的地方，或者已证明有明显优势的地方可以使用，但如果不确定，尽量还是用显式溶剂模型，毕竟GB是基于多层近似后的溶剂模型。

g_rms出来的xmgrace图有波谷！

刚完成一个蛋白在水中的simulation，总共有四个xtc文件，分别是, (5ns),

(10ns), (40ns)。然后我利用如下command组合了一个xtc:

trjcat -f -settime -o , 运行后选择的settime是0， C;

然后用g_rms -s -f -n -o ;

最后用xmgrace做出如上图。可以发现在连接处10ns的地方有个波谷，非常奇怪，问过教授，他说是因为做trjcat联合的时候应该从pr energy minimized的文件开始。一时不知道怎么回事？

小弟是初学者，还请各位帮帮忙，究竟是怎么回事？正确的做法又是什么？谢谢！

本文标签： 文件分子溶剂模拟使用