GROMACS简单指南

GROMACS官方给出的指南,翻译了一下以供参考,篇幅较长。

简单指南

这儿有一系列简单的指南帮助用户开始模拟。更为详细的实例教程可以在这网站找到:http://www.mdtutorials.com/
原网站:Short How-To guides

初学者

对于那些刚刚开始学习GROMACS或分子动力学模拟的人来说是非常艰难的。强烈建议首先阅读那些为GROMACS提供的各种广泛的文档,以及在感兴趣领域出版的文章。

资源

  • GROMACS参考手册—非常详细的文档,通常也可以作为一个非常好的MD介绍。
  • 流程图—一个蛋白在水盒子中的典型GROMACS分子动力学流程。
  • 分子动力学模拟和GROMACS介绍(幻灯片视频)—力场,积分以及温度和压力的控制(Berk Hess)

添加残基到力场中

添加新的残基

如果你有需要将新的残基添加到已经存在的力场中去以便能够使用pdb2gmx,或者修改存在的残基,那么有几个文件你应该修改。你必须参阅参考手册对于所需格式的描述部分。执行以下步骤:

  • 向你选择的力场中的rtp文件中添加残基.你也可以复制一个已经存在的残基,重命名和适当地修改它,或者你可能需要使用额外的拓扑生成工具并调整其为rtp格式。
  • 如果需要氢原子能够添加到你的残基中,请在相关的hdb文件中创建项目。
  • 如果你正引入新的原子类型,请添加它们进入 atomtypes.atpffnonbonded.itp 文件中。
  • 如果你需要一些新的成键类型,请添加它们到 ffbonded.itp 文件中。
  • 使用正确的指定(Protein,DNA,Ion等)将你的新残基添加到 residuetypes.dat 文件中。
  • 如果你的残基与其他残基涉及特殊连接,请更新 specbond.dat
    注意: 如果你正在模拟一些非自然配体在水中或者它们与正常蛋白质的结合,那么上述的工作远比产生一个独立的包含 [moleculetype] 部分的itp文件要多(例如:通过修改由一些参数化服务器产生的top文件),并将该itp文件的 #include 插入到为系统生成的那个没有非自然配体的top中。

修改力场

修改力场的最佳方法是复制已安装的forcefield目录和 residuetypes.dat 到你的工作路径下:

1
cp -r $GMXLIB/residuetypes.dat $GMXLIB/amber99sb.ff .

然后,按照上面的方法修改这些本地副本。pdb2gmx将同时找到原始版本和修改后的版本,你可以从列表中交互式地选择修改后的版本,或者如果使用pdb2gmx-ff 选项,本地版本将覆盖系统版本。

水溶剂

当使用solvate产生溶剂盒子时,你需要提供一个预先平衡好的盒子,里面装着合适的溶剂,让溶剂在你的溶质周围堆积。然后截断得到你想要的模拟体积大小。当使用3点水模型(例如 SPCSPC/ETIP3P )时你应该指定 -cs spc216.gro,这将使用 the gromacs/share/top路径下的文件。其他的水分子模型(例如 TIP4PTIP5P )也有。检查 /share/top 子目录下的内容。溶剂化以后,你应该确保在期望的温度下平衡至少5-10ps。您需要在top文件中选择正确的水模型,或者使用pdb2gmx-water 选项指定,或者手工适当地编辑top文件。
有关如何使用除纯水以外的溶剂的信息,请参阅非水溶剂化混合溶剂

非水溶剂

在GROMACS中可以使用除水以外的溶剂。唯一的要求是你有一个你所需要的溶剂预先平衡的盒子,和为这个模拟合适的参数。然后使用solvate-cs 选项即可完成溶剂化。
virtualchemistry中,可以找到一系列大约150种不同的平衡液体,这些液体经过验证可以与GROMACS一起使用,也可以用于OPLS/AA和GAFF力场。

制作非水溶剂盒子

选择一个盒子的密度和大小。尺寸不必是您最终的模拟盒的大小—一个1nm的立方体可能就可以了。成单个溶剂分子。计算出一个分子在你选择的密度和大小的盒子里的体积。使用editconf在单个分子周围放置一个同样大小的盒子。然后用editconf将分子移动到中心。然后使用genconf-rot 选项把那个盒子复制成一个大小和密度都合适的大盒子。然后使用NVT和周期性边界条件彻底平衡来去除分子的不正确排列。现在你有了一个可以传递给solvate -cs 选项,它将复制以适应实际模拟体系的大小。

混合溶剂

新用户面临的一个常见问题是如何创建一个混合溶剂(例如,在给定浓度的水中使用尿素或DMSO)的系统。完成这项工作最简单的程序如下:

  • 根据系统的盒形尺寸,确定所需的共溶剂分子数量。
  • 生成一个单分子共溶剂的坐标文件(例如urea.gro)。
  • 使用gmx insert-molecules-ci-nmol 选项将所需的共溶剂分子数目添加到盒中。
  • 使用gmx solvate 或者 gmx insert-molecules把盒子的其余部分装满水(或者其他溶剂分子)
  • 编辑你的top#include 适当的itp文件,并对 [molecules] 部分进行更改,以说明系统中的所有分子类型。

制作二硫键

最简单的方法是使用 specbond.dat 文件和pdb2gmx实现的机制。您可能会发现pdb2gmx-ss yes 非常有用。需要在同一单元的硫原子在pdb2gmx转换为一个 moleculetype,因此,可能需要正确调用pdb2gmx-chainseppdb2gmx -h以查看功能。这就要求两个硫原子之间的距离必须在+公差范围内(通常为10%),这样才能被识别为二硫原子。如果你的硫原子没有这么近,那么你也可以:

  • 编辑 specbond.dat 的内容以允许成键,并且做能量最小化时要非常小心,使键放松到一个合理的长度。
  • 在这些具有较大力常数的硫原子之间运行使用距离约束(且没有二硫键)的初步EM或MD,使它们接近现有 specbond.dat 范围,以便为第二次调用pdb2gmx提供合适的坐标文件。
    否则,手工编辑你的top文件是唯一的选择。

GROMACS运行膜的模拟

运行膜模拟

用户们在运行磷脂双分子层模拟时经常遇到问题,尤其是涉及蛋白的时候。用户在模拟膜蛋白时可能发现这个教程有用。
膜蛋白模拟一般包括以下几个步骤:

  • 1.选择含有蛋白和脂类参数的力场文件
  • 2.将蛋白插入到膜中(例如,成型双层膜使用gmx_membed,或者做粗粒化自组装模拟,然后转换成原子表示)。
  • 3.溶解体系并添加离子使之中和多余电荷,调整最终离子浓度。
  • 4.能量最小化。
  • 5.让膜适应蛋白。一般是对全部蛋白分子的重原子使用限制(1000 kJ/(mol nm2))运行大约5-10 ns的MD。
  • 6.取消限制进行平衡。
  • 7.运行成品MD。

用genbox添加水

当使用solvate在预形成的脂质膜周围添加水分子时,你可能会发现水分子被加入到了膜的空隙中。这里有几种方法可以消除这些问题,包括:

  • 短时间的MD运行,依靠疏水性形象去排除掉水。一般来说这足以实现无水疏水相,因为这些分子通常会被迅速排出而不会破坏整体结构。如果您的设置在一开始就依赖于完全的无水疏水相,您可以尝试以下建议:
  • 使用gmx solvate-radius选项去改变水的排斥半径。
  • 从你的$GMXLIB位置复制vdwradii.dat文件到你的工作目录,然后编辑它,增大脂类原子半径(碳原子一般在0.35到0.5nm)来防止溶剂进入到空隙中。
  • 手手动编辑结构以删除它们(要调整gro文件的原子计数,并考虑拓扑中的一些更修改)。。
  • 使用别人编写的脚本删除它们。

扩充材料

新分子参数化

大多数参数化问题/困难都可以很简单地解决,只要记住以下两条规:

  • 你不应该混合和匹配力场。力场(最好)被设计成自洽的,通常不能很好地与其他力场协同工作。如果你用一个力场来模拟系统的一部分,另一部分用另一个不同的力场来模拟而不是用第一个力场来参数化,您的结果可能会有问题,审阅人员会关注。选择一个力场,就要用那一个力场。
  • 如果您需要开发新的参数,请按照原力场其余部分的推导方式推导它们,这意味着你需要看原始文献。推导力场参数没有单一正确的方法;你需要的是推导出与力场其余部分一致的参数。怎么做取决于你想用哪个力场。例如,利用AMBER力场,推导非标准氨基酸的参数可能需要进行许多不同的量子计算,而推导GROMOS或OPLS参数可能涉及更多 (a) 拟合各种流体和液相性质。 (b) 根据经验/化学知识/类比调整参数。这里可以找到一些关于自动化方法的建议。

在尝试参数化新力场或者现有力场的新分子之前,对GROMACS有一定的仿真经验是明智的,这些都是专家话题,不适合给本科生做研究项目,除非你喜欢昂贵的准随机数生成器。需要非常全面地了解GROMACS参考手册第5章。如果你还没有足够了解,请阅读下面关于外来物种参数化的内容。
另一个建议是:在获取参数时,不要过于随意,否则就会买到高档珠宝。仅仅因为街上有人愿意以10美元的价格卖给你一条钻石项链,并不意味着你应该在那里买一条。同样,从你从未听说过的人的网站上下载你感兴趣的分子的参数也不一定是最好的策略,尤其是如果他们不解释他们是如何得到这些参数的。
预先警告使用PRODRG拓扑而不验证其内容:现已出版,以及一些技巧,以正确地推导参数的GROMOS力场。

外来物种

所以,你想要模拟一个蛋白质/核酸系统,但是它会结合各种其他的金属离子(钌?),或者有一个铁硫团簇对其功能至关重要,或者相似的。但是,(不幸的是?)在你想要使用的力场中没有这些参数可用。你应该怎么做?您可以向GROMACS用户的电子邮件列表发送电子邮件,并获得常见问题解答。
如果你真的坚持在分子动力学中模拟它们,你需要获得它们的参数,或者从文献中,或者通过你自己的参数化。但在这之前,停下来想一想可能是很重要的,因为有时可能存在这样的原子/团簇参数不存在的原因。特别地,这里有几个基本的问题,你可以问你自己,看看是否合理地开发/获得这些标准参数,并在分子动力学中使用它们:

  • 量子效应(即电荷转移)可能很重要吗?(即。,如果你有一个二价金属离子在酶活性位点,并有兴趣研究酶的功能,这可能是一个巨大的问题)。
  • 在我所选择的力场中使用的标准力场参数化技术在这种类型的原子/簇中可能会失败吗?(例如,Hartree-Fock 6-31G*不能充分描述过渡金属)。

如果这两个问题的答案都是“是”,你可能会考虑用经典分子动力学以外的东西来做模拟。
即使这两个问题的答案都是“不”,在尝试自己的参数化之前,您可能也想咨询一下您感兴趣的化合物的专家。此外,在开始这些操作之前,您可能想尝试对一些更直接的东西进行参数化。

平均力势

平均力势(PMF)定义为:在给定系统的所有构型上给出平均力的势。在GROMACS中有几种计算PMF的方法,其中最常见的可能是使用pull代码。使用伞形采样,它允许对统计上不可能的状态进行抽样,获得PMF的步骤如下:

  • 沿着反应坐标生成一系列构象(从SMD模拟、常规MD模拟或任意生成的构象)
  • 使用伞形抽样在抽样窗口内限制这些构象。
  • 利用gmx wham,使用WHAM算法重构PMF曲线。

这里链接了一个更详细的教程,用于伞形抽样。

单点能

计算单个构象的能量有时是有用。GROMACS的最佳方法是使用mdrun-rerun选项,该选项将tpr中的物理模型应用于轨迹或坐标文件中的构象,并提供给mdrun。

1
gmx mdrun -s input.tpr -rerun configuration.pdb

注意,所提供的构象必须与使用grompp生成tpr文件时使用的拓扑相匹配。您提供给grompp的构象除了原子名称之外,都无关紧要。您还可以将此功能用于能量组(看参考手册),或者用多构型轨迹(在本例中,默认情况下mdrun将对每个构象执行邻居搜索,因为它不能假设输入是相似的)。

零步能量最小化在报告能量之前执行一个步骤,而零步MD运行具有(可避免的)并发症,这些并发症与在存在约束的情况下可能重新启动有关,因此不推荐使用这两种程序。

碳纳米管

Robert Johnson的技巧

采纳了Robert Johnson在gmx-users mailing list上的帖子。

  • 要绝对保证拓扑文件中“末端”碳原子是共享一个键。
  • gmx grompp时的mdp输入文件要使用periodic_molecules = yes
  • 即使拓扑文件正确,如果你把碳纳米管放置在错误尺寸的盒子里将使得扭曲发生。因此使用VMD可视化纳米管和它的周期性成像,确保镜像之间的空隙正确。如果空隙太小或太大,在纳米管上将有巨大的压力导致其扭曲或拉伸。
  • 沿着碳纳米管的轴线方向上不要使用压力耦合。事实上,为了调试目的,如果发生某些错误,可能先最好关掉整个压力耦合直到你弄清楚是什么问题。
  • 当使用具有特定力场的x2top时,假设分子的连接性。如果是周期性的,亦或者是非周期性上面有氢原子,那么你的纳米管末端的碳原子最多只能和两个其他分子相结合。
  • 你可以使用x2top工具中的 -pbc 选项产生一个无限长的纳米管,这里x2top将识别为末端C原子共享一个化学键。因此,当你使用grompp时不会得到关于单个碳原子的报错。

Andrea Minoia的指南

GROMACS模拟碳纳米管(也可以参看:http://www.webcitation.org/66u2xJJ3O

原文翻译

使用GROMACS模拟碳纳米管

网上有许多关于GROMACS和CNTs的文献(1,2,3)。即使读了很多,我也没有找到我想要的:一个清楚、简单的模拟单壁碳纳米管和多壁碳纳米管(SWNT,MWNT)指南,周期或者非周期性。
我花了一段时间才找到步骤做这些我想做的事情:构建一个周期性CNT并创建GROMACS的拓扑文件。

构建和准备CNT

第一步很明显是构建CNT,例如采用 YASC buildCstruct。使用这个插件可以构建SWCNT和MSWCNT,包括armchair或zigzag型,有限的或周期性的。从 v1.1版本开始,buildCstruct能够直接保存为grmacs的GRO格式文件。
注意要有一个足够大的晶胞容纳所有的CNT,否则你将不能产生拓扑文件。你能够使用ngmx,molden或者vmd的pbc检查结构。而且,也要注意选择一个足够大的盒子去容纳其他你想和CNT接触的分子,例如聚合物链。

为x2top准备的文件

x2top需要一堆文件来创建拓扑文件。我想使用oplsaa力场。但是我不想搞乱share/gromacs/top目录中的原始ffoplsaa文件,因此我自己创建了我需要的文件。

.n2t itp

此文件用于在拓扑中转换名称。根据它的连接性读取结构中原子的名称(例如成键数,键长和原子结合),根据.n2t文件中定义的尝试去猜合适的原子类型。我的基于oplsaa力场的石墨烯和CNT的.n2t文件如下:

1
2
3
4
5
6
7
8
; Oplsaa-based n2t for carbon-based structures such as CNTs and graphenes
; Andrea Minoia
H HJ 0.00 1.008 1 C 0.109 ;Hydrogen
C CJ 0.00 12.011 3 C 0.142 H 0.109 H 0.109 ;Periferic C
C CJ 0.00 12.011 3 C 0.142 C 0.142 H 0.108 ;Periferic C
C CJ 0.00 12.011 1 C 0.142 ;Internal/periodic C
C CJ 0.00 12.011 2 C 0.142 C 0.142 ;Internal/periodic C
C CJ 0.00 12.011 3 C 0.142 C 0.142 C 0.142 ;Internal/periodic C

我决定把所用的电荷设置为0,但是你可能想不一样。

.itp文件

我更喜欢不让x2top将力场参数放入拓扑文件中,我已经定义了我自己的.itp文件:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
; Oplsaa-based force field for carbon-based structures such as CNTs and graphenes
; Andrea Minoia

[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 3 yes 0.5 0.5
; parameters are taken from the OPLS force field

[ atomtypes ]
; The charges here will be overwritten by those in the rtp file
; name mass charge ptype sigma eps
CJ 6 12.01100 0.000 A 3.55000e-01 2.92880e-01 ;opls_147 naftalene fusion C9
HJ 1 1.00800 0.000 A 2.42000e-01 1.25520e-01 ;opls_146 HA hydrogen benzene. I have set the charges zero

[ bondtypes ]
; i j func b0 kb
CJ CJ 1 0.14000 392459.2 ; TRP,TYR,PHE
CJ HJ 1 0.10800 307105.6 ; PHE, etc.

[ angletypes ]
CJ CJ CJ 1 120.000 527.184 ; PHE(OL)
CJ CJ HJ 1 120.000 292.880 ;
HJ CJ HJ 1 117.000 292.880 ; wlj from HC-CM-HC

[ dihedraltypes ]
CJ CJ CJ CJ 3 30.33400 0.00000 -30.33400 0.00000 0.00000 0.00000 ; aromatic ring
HJ CJ CJ HJ 3 30.33400 0.00000 -30.33400 0.00000 0.00000 0.00000 ; aromatic ring
HJ CJ CJ CJ 3 30.33400 0.00000 -30.33400 0.00000 0.00000 0.00000 ; aromatic ring

目前没有二面角。

.rtp 文件

最后,我需要创建oplsaa的残基拓扑文件,这里还没有定义任何残基。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
; New format introduced in Gromacs 3.1.4.
; Dont use this forcefield with earlier versions.

; Oplsaa-based rtp for carbon-based structures such as CNTs and graphenes
; Andrea Minoia

; NB: OPLS chargegroups are not strictly neutral, since we mainly
; use them to optimize the neighborsearching. For accurate simulations
; you should use PME.

[ bondedtypes ]
; Col 1: Type of bond
; Col 2: Type of angles
; Col 3: Type of proper dihedrals
; Col 4: Type of improper dihedrals
; Col 5: Generate all dihedrals if 1, only heavy atoms of 0.
; Col 6: Number of excluded neighbors for nonbonded interactions
; Col 7: Generate 1,4 interactions between pairs of hydrogens if 1
; Col 8: Remove propers over the same bond as an improper if it is 1
; bonds angles dihedrals impropers all_dihedrals nrexcl HH14 RemoveDih
1 1 3 1 0 3 1 0

FF.dat文件

这里我指定了力场的名字:

1
2
1
ffcnt_oplsaa

使用x2top创建拓扑

一旦所有文件准备好,就可以使用x2top产生拓扑:

1
x2top -f file.gro -o topol.top -ff cnt_oplsaa -name CNT -noparam -pbc

-pbc选项将产生那些周期性边界条件的所有键,角度和二面角(这就是为什么.gro文件中需要指定合适大小的盒子单元)。如果你对纳米管不是周期性的,那就不需要 -pbc

矫正拓扑

最后一步是将用来描述二面体的函数从1更改为3

.mdp文件

在你的mdp文件中打开周期性选项:

1
2
periodic_molecules=yes
pbc = xyz

这就是全部…现在你已经可以模拟周期性碳纳米管或碳基结构。

溶解CNT

为了溶解CNT,我们能够使用GROMACS的genbox工具(gromacs5.0版本以上改为了solvate)采用-cp(纳米管的gro文件)和-cs(溶解的盒子gro文件)。
通常,我首先在NPT下平衡溶剂盒子,然后创建了一个比CNT盒子更大的超胞。这是因为溶质盒子(-cp)将成为最终系统的盒子。
genbox将尝试用溶剂分子填充盒子空隙,甚至在CNT的内部。如果这不可取,我们能写一个创建索引文件的脚本,该文件不包含CNT内部分子中原子的原子索引。然后,使用editconf去产生最终体系:

1
editconf -f starting_system.gro -n index.ndx -o final_system.gro

使用YSAC nosolincnt去除CNT内部的原子

与editconf一起使用的引索文件是为了去除CNT内部的溶剂分子,这能够佷容易通过YSAC nosolincnt编写,它是一个VMD插件。
用VMD读取体系后,打开TCL控制台,执行插件:

1
source /path/of/the/script/nosolincnt.tcl

将提示你插入CNT的片段数和CNT半径(埃米)。这些信息将用于创建一个圆柱形,根据CNT质心X,Z坐标得到的半径。当视图窗口处于激活状态时,你能按 C 去检查选择的结果。如果一切正确,然后按 W 将会写入引索文件而不需要选择原子。
在下面的图片中,分别显示了溶剂化后的碳纳米管、要除去的分子和最终的体系。
CNT体系

复合物实例:溶解CNT+聚合物

我真的不知道有什么兴趣可以单独建立一个溶剂化碳纳米管模型,当然,如果我们加入聚合物,事情会变得更有趣。在下面的流程图中,我发现溶剂化CNT与聚合物链相互作用。
CNT+聚合物模拟流程图

故障排除,提示和技巧

Grompp时提示原子名不匹配

如果您以pdb文件开始生成gro文件,根据您如何生成pdb文件,可能会丢失残基名称。在这种情况下,我经历了grompp的一系列警告:

1
Warning: atom name X in _FILE.top_ and _FILE.gro_ does not match...

我通过在.gro和.top文件中添加一个残基名称来解决这个问题。

VI编辑

学习它!它是您最好的伙伴,特别是在使用gromacs和大型文本文件时。


包括了使用OPLS-AA力场设置一个CNT模拟的详细步骤。一个CNT的结构可以很容易通过buildCstruct(此Python脚本也可以加氢)或者通过TubeGen Online(只需要复制和粘贴PDB输出到文件并命名为cnt.pdb)。
为了能够用GROMACS模拟,你可能要执行以下操作:

  • 新建一个cnt_oplsaa.ff目录文件夹
  • 在此目录下,创建如下文件,以下步骤来源于上面的指南:
    • forcefield.itp来自itp部分文件。
    • atomnames2types.n2t来自n2t部分文件。
    • aminoacids.rtp来自rtp部分文件。
  • 使用传统的力场产生拓扑(cnt_oplsaa.ff目录文件必须与gmx x2top命令执行的目录相同,亦或者它能够在GMXLIB路径中找到),采用gmx x2top-noparam 选项,不使用来自命令行(-kb, -ka, -kd)指定的bond/angle/dihedral等常数,而是依赖力场文件;然而,这就需要下一步(确定二面角函数)
    1
    gmx x2top -f cnt.gro -o cnt.top -ff cnt_oplsaa -name CNT -noparam

二面角的函数类型通过gmx x2top设置为“1”,但是力场文件指定的函数类型为“3”。因此,替换拓扑文件中 [ dihedrals ] 部分的”1“为”3“。一种快速的方式是使用 sed(但是你可能需要调整你的操作系统;也要手动看一下top文件核查一下你更改的二面角函数类型):

1
sed -i~ '/\[ dihedrals \]/,/\[ system \]/s/1 *$/3/' cnt.top

一旦你有了拓扑文件,你就可以设置你的系统了。例如,一个简单的 真空模拟(在em.mdpmd.mdp中使用你自己的参数):
放入一个稍微大的盒子中:

1
gmx editconf -f cnt.gro -o boxed.gro -bt dodecahedron -d 1

真空中的能量最小化:

1
2
gmx grompp -f em.mdp -c boxed.gro -p cnt.top -o em.tpr
gmx mdrun -v -deffnm em

真空中的MD:

1
2
gmx grompp -f md.mdp -c em.gro -p cnt.top -o md.tpr
gmx mdrun -v -deffnm md

查看轨迹:

1
2
3
gmx trjconv -f md.xtc -s md.tpr -o md_centered.xtc -pbc mol -center
gmx trjconv -s md.tpr -f md_centered.xtc -o md_fit.xtc -fit rot+trans
vmd em.gro md_fit.xtc

可视化软件

为了可视化轨迹文件/坐标文件,一些程序是有用:

  • VMD—-一种分子可视化程序,用于显示、动画和分析大型生物分子系统,使用三维图形和内置脚本。读取GROMACS轨迹。
  • PyMOL—-支持动画、高质量渲染、晶体学和其他常见分子图形活动的功能强大的分子查看器。不读取默认配置中的GROMACS轨迹,需要转换为PDB或类似格式。当使用VMD插件编译时,可以加载trr & xtc文件。
  • Rasmol—-派生软件Protein Explorer(如下)可能是更好的选择,但是Chime组件需要windows。Rasmol在Unix上工作得很好。
  • Protein Explorer—-一个RasMol衍生物,是最容易使用和最强大的软件,以查看大分子结构及其与功能的关系。它运行在Windows或Macintosh/PPC电脑上。
  • Chimera—-一个功能齐全的,基于python的可视化程序,具有各种功能,可以在任何平台上使用。当前版本读取GROMACS轨迹。
  • Molscript—-这是一个脚本驱动程序形式的高质量显示分子三维结构的示意图和详细的表示。你可以从Avatar免费获得学术许可。

此外,如果在配置时找到了合适的库,gmx view也会很有用。

拓扑成键 vs 渲染的键

请记住,这些可视化工具只查看您提供的坐标文件(除非您提供gmx view一个tpr文件)。因此它没有使用拓扑文件或tpr文件中描述的拓扑。这些程序中的每一个都对化学键的位置做出自己的猜测,以便进行呈现,所以如果启发式并不总是与您的拓扑结构匹配,不要感到惊讶。

提取轨迹信息

在GROMACS轨迹文件(trrxtctng)中查找信息有几种可用的技术。

  • 使用GROMACS轨迹分析工具。
  • 使用gmx traj编写一个xvg文件,并像上面那样在外部程序中读取该文件。
  • 使用 gromacs/share/template/template.cpp 作为模板编写自己的C代码。
  • 使用gmx dump将shell输出重定向到一个文件,然后在外部程序(如MATLAB、Mathematica或其他电子表格软件)中读取该文件。
  • 资源:GROMACS中MD轨迹分析(幻灯片视频)- (David van der Spoel)。

外部工具执行轨迹分析

近年来,一些外部工具已经足够成熟,可以从多个仿真包中分析不同的轨迹数据集。下面是一个简短的列表,已知能够分析GROMACS轨迹数据的工具。

绘制数据

各种GROMACS分析实用程序可以生成xvg文件。这些是专门为直接在Grace中使用而格式化的文本文件。但是,您可以在所有GROMACS分析程序中使用-xvg none选项运行程序来关闭Grace特定的代码。这就避免了gnuplot和Excel等工具的问题(参见下面)。

注意,Grace使用一些嵌入的反斜杠代码来表示单位中的上标、普通脚本等。面积(Area (nmS2N))等于nm²。

软件

一些软件包,可用于图形数据在一个xvg文件:

  • Grace —- WYSIWYG 2D绘图工具,用于X窗口系统和M*tif。Grace几乎可以运行在任何版本的类unix操作系统上,只要您能够满足它的库依赖关系(Lesstif是Motif的一个有效的免费替代品)。它也适用于其他常见的操作系统。
  • gnuplot—-可移植的命令行驱动的交互式数据和函数绘图工具,适用于UNIX、IBM OS/2、MS Windows、DOS、Macintosh、VMS、Atari和许多其他平台。记住要使用:set datafile commentschars "#@&"去避免gnuplot试图在xvg文件中解释Grace特定的命令,或者在运行分析程序时使用-xvg none选项。对于简单的使用:plot "file.xvg" using 1:2 with lines会得到正确的结果.
  • MS Excel—-将文件扩展名更改为.csv并打开文件(当提示时,选择忽略前20行左右并选择固定宽度的列,如果使用德语MS Excel版本,则必须将十进制分隔符从“,”更改为“.”。,或者使用你最喜欢的*nix工具。
  • Sigma Plot—-一个用于windows的商业工具,其中包含一些有用的分析工具、
  • R—-为统计计算和图形提供了免费的语言和环境,提供了各种各样的统计和图形技术:线性和非线性建模、统计测试、时间序列分析、分类、聚类等。
  • SPSS—-一个商业工具(统计产品和服务解决方案),它也可以绘制和分析数据。

胶束聚类

如果您有一个完全形成的单一聚集体,并希望为该聚集体或该聚集体周围的溶剂生成空间分布函数,gmx spatial工具是必要的.
在计算诸如旋转半径和径向分布函数等性质之前,确保胶束不会在周期边界条件边界上断裂是一个必不可少的步骤。如果没有这个步骤,您的结果将是不正确的(这个错误的一个迹象是,当可视化轨迹看起来很好时,计算值会出现无法解释的巨大波动)。

需要三个步骤:

  • 使用trjconv-pbc cluster获得一个包含单元细胞中所有脂类的单帧。这必须是轨迹的第一帧。以前某个时间点的类似帧将不起作用。
  • 使用grompp根据上述步骤输出的帧创建一个新的tpr文件。
  • 使用新生成的tpr文件使用trjconv-pbc nojump生成所需的轨迹。

更明确地说,同样的步骤是:

1
2
3
gmx trjconv -f a.xtc -o a_cluster.gro -e 0.001 -pbc cluster
gmx grompp -f a.mdp -c a_cluster.gro -o a_cluster.tpr
gmx trjconv -f a.xtc -o a_cluster.xtc -s a_cluster.tpr -pbc nojump

补充

源资源来自:http://www.gromacs.org/Documentation/How-tos

注意:下面的指南有很多的命令已经更改,一些项目也已经太旧。

将3点水模型改为4点水模型

如果你有一个3点的水模型的.gro文件,其中,并且你希望保持水的坐标,但是要更改为一个4点的水模型,请按照以下步骤:

  • cat a.gro |awk '{print $0; if($2=="HW2")printf("%8s MW %4d%8.3f%8.3f%8.3f\n",$1,$3,$4,$5,$6)}' > b.gro
  • 修改文件b.gro第二行中的原子数
  • gmx editconf -f b.gro -o c.gro
  • c.gro运行一个短的能量极小化。

检查点作业

4.0版本开始,GROMACS有内置的自动检查点,见模拟重启
这儿有一系列通用插件为GROMACS 3.x检查点工作,并自动重新提交。

免责声明

虽然作者使用这些脚本并认为它们是高质量的,但并不保证它们的准确性。此外,由于这是一个wiki站点,下面出现的脚本可能已经从原始内容进行了修改。使用这些脚本的风险由您自己承担。

预计时间消耗

如果您是bash脚本编写专家,那么应该不会花费太长时间。否则,计划花一天的时间让它运行起来。请阅读这篇文章。

概述

  • 1.head.sh脚本像在交互式作业中一样从命令行运行。此脚本用于提交一些链接在一起的作业,以便在第N个作业完成后启动第N+1个作业。head.sh提交一连串的sub.$CLUSTER.sh作业,这里的sub.$CLUSTER.sh只是一个sub.sh包装,sub.sh才是脚本的真正内容,不需要定制。
  • 2.每次运行必须自定义head.sh文件顶部的部分。

为什么如此复杂

由于作者在不同的集群上遇到了问题,所以对脚本进行了扩充,以自动处理这些问题。

  • 1.该脚本优雅地处理NFS延迟问题。
  • 2.该脚本优雅地从崩溃中恢复
  • 3.该脚本检查.xtc文件,以确保它们是完整的,没有损坏。
  • 4.该脚本只保留运行中最重要的文件(.edr, .xtc, 初始.gro),为了减少硬盘的使用。

你需要的程序

如果希望使用排序,则需要在用户提交页面上提供G_DESORT包。目前,排序与混排相结合,没有选项不排序就混排的选项(但这将是一个简单的脚本修改)。如果你不想混排或排序,那么你可以在head.sh中设置NEVER_USE_SORT_SHUFFLE为非零值。还要注意head.sh中的DD说明符指定了查找g_desort的目录。关于g_desort的进一步说明:作者使用了它,但是开发人员不支持它,所以使用它的风险由您自己承担。

这些脚本在集群上特定的修改

用户需要为他们的排队系统定制这些脚本。为了设置你的环境,找到所有字符串的情况:

1
case "$CLUSTER" in

并为您的集群添加一个包含适当设置的名称节。例如,在名为bigCluster的cluser上使用sqsub而不是qsub。您需要修改head.sh的以下部分,使原来的:

1
2
3
4
5
case "$CLUSTER" in
*)
subCommand="qsub"
;;
esac

变成:

1
2
3
4
5
6
7
8
case "$CLUSTER" in
bigCluster)
subCommand="sqsub"
;;
*)
subCommand="qsub"
;;
esac

然后在head.sh顶部你要定义:

1
CLUSTER=bigCluster

一旦您完成了这些更改,就不需要进行太多的修改来开始新的运行。

迭代过程是如何发生的

  1. 用户创建两个文件
    1
    2
    echo 0 > finished_test
    echo 0 > finished_next_start_time

其中finished_test文件表示第0次迭代已经完成,finished_next_start_time文件表示下一轮迭代将从0 ps开始。

  1. 用户创建一个名为md0_success的目录,并将其起始.gro文件放在其中。.gro文件的命名很重要。假设你的MYMOL(定义在head.sh的顶部)是:
    1
    MYMOL=proteinInWater

然后,md0_success目录必须包含一个名为proteininwater_md0_deshuffleddesord .gro的文件。

1
2
$ls md0_success
md0_success/proteinInWater_md0_deshuffleddesorted.gro

即使不使用排序或混排,名称的“deshuffleddes”部分也是必需的。

  1. 创建.top.ndx文件,并将它们命名为$MYMOL.top$MYMOL.ndx
  2. 创建您的.mdp文件,并将其命名为$MYMOL.mdp
    4.a) mdp文件的最后一行必须是:EOF
    4.b) 您的.mdp文件必须包含一些标志,这些标志将根据head.sh中的选项使用sed替换。这些选项必须用特殊标志定义。确保你的.mdp文件包含这样的行:

    1
    2
    3
    4
    5
    6
    7
    cpp                 =  CPP
    nsteps = NSTEPS
    tinit = TINIT
    nstxout = NSTEPS
    nstvout = NSTEPS
    nstfout = NSTEPS
    nstxtcout = SAVE_FREQUENCY
  3. cp $MYMOL.mdp $MYMOL_posre.mdp,然后修改新文件,使它可以用于head.sh中定义的N小于$NJUMPS_POSRE的任何段。就我个人而言,作者为第一个片段设置了位置约束的定义,并在两个.mdp文件中分别设置了gen_velunconstrained_start,以便在第一个片段上生成速度,并在其他片段上适当地重新启动。
    注意:这个脚本从未使用tpbconv,因为它是为重新排序而设计的。但是,如果您想要使用tpbconv,那么修改它应该是一件简单的事情。

  4. 所有脚本上执行chmod +x

脚本的位置

A) head.sh.mdp .ndx .top文件以及包含.gro文件的md0_success目录一起在工作目录。
B) sub.shsub.$CLUSTER.sh文件位于/home/`whoami`/scripts/submission/中,因为许多不同的作业都应该使用相同的sub.sh脚本。您将需要修改head.sh的底部,在那里它提交sub.sh作业。

运行脚本

cd到你的工作目录并:

1
nohup ./head.sh &

运行时创建的文件

tmp.sub: 用于捕获作业信息的提交的输出
wait_for_this_pid: 链中最后一个作业的pid
nohup.out: 按正确的链接顺序提交的带有pid的作业列表
DATA/ directory –>这将存储来自N-3和更老版本的.xtc.edr文件,其中N是当前正在运行的尚未完成的作业。注意,N-1N-2目录包含所有文件,但是大多数这些文件在归档到数据目录时都会被删除。如果您想保存更多的文件,那么应该在“finished_test”部分修改sub.sh的源代码。

向已经运行的链提交更多作业

只需要再次运行head.sh。由于这个原因,您不应该在任何时候删除wait_for_this_pid文件。注意,如果删除了该文件,则可以使用nohup来确定链中最后一个作业的pid,并将其放入wait_for_this_pid中。

head.sh脚本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
$cat head.sh
# Automated submission script
# Chris Neale November 2007

#!/bin/bash
PATH=$PATH:.

##########################################################
# Basic setup options:
MYMOL=pagagg # starting name for your .gro, .mdp, .top (etc) files
REPLICA=1 # Just a flag to differentiate runs if you are running repeat identical simulations at once.
MD=/my/directory # working directory
CLUSTER=cluster1 # a tag used to define cluster specific characteristics for ease of porting netween clusters
NUM_TO_SUBMIT=1 # number of chained jobs to submit
MYNP=4 # number of CPUs to occupy
TJUMP=500 # Time (in ps) for a single step
NJUMPS=1000 # Number of steps before dying and allowing a new chained job to start
NFRAMES_IN_XTC=51 # Number of frames that you expect in the .xtc file = TJUMP/(SAVE_FREQUENCY*dt) + 1
SAVE_FREQUENCY=5000 # Save the .xtc every this number of timesteps
NJUMPS_POSRE=1 # The number of initial steps that will be done by position restraints


##########################################################
# More advanced setup options:

RUN_AS_A_TEST=0 # use a test queue
NEVER_USE_SORT_SHUFFLE=0 # set to nonzero to avoid sorting and shuffling always
RUNTIME_LIMIT=1w # flag to apply a runtime limit to your job
APPLY_RUNTIME_LIMIT=0 # set to nonzero to actually apply a runtime limit
DOUBLE_CHECK_DESORT=0 # currently not implemented

##########################################################
# Things below this line do not usually need to be changed

cd ${MD}

if((RUN_AS_A_TEST)); then
# Override the previously setup options
MYNP=2
TJUMP=0.2
NJUMPS=1
NFRAMES_IN_XTC=2
SAVE_FREQUENCY=100
NJUMPS_POSRE=0
testFlag="--test"
else
testFlag=""
fi

case "$CLUSTER" in
*)
if((APPLY_RUNTIME_LIMIT==1)); then
timeFlag="-r $RUNTIME_LIMIT"
else
timeFlag=""
fi
;;
esac

case "$CLUSTER" in
*)
subCommand="qsub"
;;
esac

if((MYNP>1)); then
case "$CLUSTER" in
*)
queueIdent="-q mpi --nompirun"
;;
esac
else
case "$CLUSTER" in
*)
queueIdent="";;
esac
fi

for((i=0;i<NUM_TO_SUBMIT; i++)); do

if [ -e DO_NOT_RUN ]; then
echo "ERROR: file DO_NOT_RUN exists... exiting"
exit
fi

if [ -e wait_for_this_pid ]; then
pid=`cat wait_for_this_pid | awk '{print $1}'`
case "$CLUSTER" in
*)
waitFlag="-w $pid"
;;
esac
else
waitFlag=""
fi

case "$CLUSTER" in
*)
nameFlag="-N ${MYMOL}_${REPLICA}"
mynpFlag="-n ${MYNP}"
subName=".${CLUSTER}"
;;
esac

${subCommand} ${testFlag} ${waitFlag} ${timeFlag} ${queueIdent} ${nameFlag} ${mynpFlag} -o `pwd`/out.${i} -e `pwd`/err.${i} /home/`whoami`/scripts/submission/sub${subName}.sh ${MD} ${MYNP} ${MYMOL} ${TJUMP} ${NJUMPS} ${NFRAMES_IN_XTC} ${SAVE_FREQUENCY} ${NJUMPS_POSRE} ${NEVER_USE_SORT_SHUFFLE} ${CLUSTER} ${DOUBLE_CHECK_DESORT} > tmp.sub 2>&1
case "$CLUSTER" in
*)
pid=`tail -n 1 tmp.sub | awk '{print $4}'`
;;
esac

echo "Submitting ${subCommand} ${testFlag} ${waitFlag} ${timeFlag} ${queueIdent} ${nameFlag} ${mynpFlag} -o `pwd`/out.${i} -e `pwd`/err.${i} /home/`whoami`/scripts/submission/sub${subName}.sh ${MD} ${MYNP} ${MYMOL} ${TJUMP} ${NJUMPS} ${NFRAMES_IN_XTC} ${SAVE_FREQUENCY} ${NJUMPS_POSRE} ${NEVER_USE_SORT_SHUFFLE} ${CLUSTER} ${DOUBLE_CHECK_DESORT} -- received pid = $pid"
echo "$pid" > wait_for_this_pid
sleep 1

done

sub.$CLUSTER.sh脚本

1
2
3
4
5
6
7
8
9
10
$ cat sub.cluster1.sh 
#!/bin/bash

# This is a wrapper script so that you can define PATH, LD_LIBRARY_PATH
# A wrapper is required if you need to define via "#$ -v PATH=", etc...

/home/`whoami`/scripts/submission/sub.sh ${1} ${2} ${3} ${4} ${5} ${6} ${7} ${8} ${9} ${10}

# If you have severe NFS delay, you may need to include this next line
#sleep 60

sub.sh脚本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
$ cat sub.sh 
# Automated submission script
# Chris Neale November 2007

echo -n "TIMING TEST (start): "
date
#!/bin/bash
MD="$1"
MYNP="$2"
MYMOL="$3"
TJUMP="$4"
NJUMPS="$5"
NFRAMES_IN_XTC="$6"
SAVE_FREQUENCY="$7"
NJUMPS_POSRE="$8"
NEVER_USE_SORT_SHUFFLE="$9"
CLUSTER="${10}"
DOUBLE_CHECK_DESORT="${11}"

case "$CLUSTER" in
*)
ED=/tools/gromacs-3.3.1/exec/bin
PED=${ED}
mpiLocation="/tools/openmpi/1.2.1"
mpirunProg="${mpiLocation}/bin/mpirun"
mdrun_mpiProg="mdrun_openmpi_v1.2.1"

###############################################
# For LAM mpi #
# mpiLocation="/tools/lam/lam-7.1.2" #
# mpirunProg="${mpiLocation}/bin/mpirun C" #
# mdrun_mpiProg="mdrun_mpi" #
###############################################

;;
esac

case "$CLUSTER" in
*)
DD=/home/`whoami`/gromacs/template
;;
esac

# The NEW variable allows use of a special writing location (e.g ${TMPDIR} or /scratch/`whoami`)
case "$CLUSTER" in
*)
NEW="."
;;
esac

case "$CLUSTER" in
*)
CPP="cpp"
;;
esac

PATH=$PATH:.

cd ${MD}

TINY_SLEEP=1
SHORT_SLEEP=10
LONG_SLEEP=60
EXTENDED_SLEEP=300

###############################################
# Startup tests:

if [ -e DO_NOT_RUN ]; then
echo "ERROR error: file DO_NOT_RUN exists... exiting"
exit
fi

for((v=0;v<2;v++)); do
num=0;
if [ -e finished_grompp ]; then
let "num=$num+1"
fi
if [ -e finished_mdrun ]; then
let "num=$num+1"
fi
if [ -e finished_desort ]; then
let "num=$num+1"
fi
if [ -e finished_test ]; then
let "num=$num+1"
fi

if((num==0)); then
echo "Unsure how to start the run. Check this out"
echo "$ls -l finished_grompp finished_mdrun finished_desort finished_test"
ls -l finished_grompp finished_mdrun finished_desort finished_test
echo "Perhaps you forgot to set a finished_XXX file upon starting your run?"
echo " - Otherwise there seems to be an error in the script."
else
if((num!=1)); then
echo "Unsure how to start the run. Check this out"
echo "$ls -l finished_grompp finished_mdrun finished_desort finished_test"
ls -l finished_grompp finished_mdrun finished_desort finished_test
echo "Only one of these files should have existed."
fi
fi
if((num==1)); then
break;
fi
echo "Will sleep then try one more time"
sleep {$LONG_SLEEP}
done

if((num!=1)); then
echo "Could not resolve multiple x for finished_x problem. Exiting"
touch ./DO_NOT_RUN
exit
fi

###############################################
# Initializations:

gromppProblemsInARow=0
reverted=0
MAX_CONSECUTIVE_GROMPP_ERRORS=2
MAX_CONSECUTIVE_MDRUN_ERRORS=2

function waitForExistNotEmpty {
# First arg controls usage:
# 0 for existance test
# 1 for not empty test
# -1 for must not exist test
# 12 for not empty test plus require single value to equal third arg
# Second arg is name of file/directory
# Third arg is the expected single value in the file if First arg is = 12
#
# Note: This overly complicated procedure is required for proper usage of a
# cluster where NFS delay can be significant and simple -s tests
# routinely fail to detect the fact that the file is empty as far as
# val=`cat file` is concerned
notEmpty="$1"
case "$notEmpty" in
1*)
eneFlag="-s $2"
;;
-1)
eneFlag="! -e $2"
;;
0)
eneFlag="-e $2"
;;
*)
echo "ERROR error: incorrect argument to waitForExistNotEmpty = $notEmpty"
exit
;;
esac

for((length=1;length<100;length++)); do
if [ ${eneFlag} ]; then
break
fi
sleep ${length}
done
if((length==100)); then
echo "ERROR error: Have slept for 30 minutes while waiting for the file $2 to meet conditions [ ${eneFlag} ] ... is there a problem in your script or is the NFS delay very very large?"
fi

# for all uses other than First Arg = 12 this function is over
if((notEmpty!=12)); then
return
fi
expectedVal="$3"
for((length=1;length<100;length++)); do
currentVal=`cat $2`
if [ -n "$currentVal" ]; then
# make sure variable is non-empty before making the comparison
case "$currentVal" in
$expectedVal)
echo "NOTE: breaking from loop since currentVal($currentVal)=expectedVal($expectedVal)"
break
;;
esac
fi
sleep ${length}
done
if((length==100)); then
echo "ERROR error: Have slept for 30 minutes while waiting for the file $2 to meet conditions `cat $2 = $expectedVal` ... is there a problem in your script or is the NFS delay very very large?"
fi
}

###############################################
# The main loop:

for ((njump=0;njump<NJUMPS;njump++)); do

# GROMPP (WITH SORTING)
if [ -e finished_test ]; then
NIN=`cat finished_test`
TINIT=`cat finished_next_start_time`
let "NOUT=$NIN+1"
DIR=md${NOUT}_running
PREV=md${NIN}_success
if [ ! -e ${PREV} ]; then
echo "There was some problem. Expected ${PREV} to exist, but it does not"
touch ./DO_NOT_RUN
exit
fi
if [ ! -e ${NEW}/${DIR} ]; then
mkdir ${NEW}/${DIR}
fi

nsteps=`echo "$TJUMP/0.002" | bc -l | awk -F '.' '{print $1}'`
if((NOUT<=NJUMPS_POSRE)); then
sed "s/TINIT/${TINIT}/" ${MYMOL}_posre.mdp | sed "s/NSTEPS/${nsteps}/" | sed "s/SAVE_FREQUENCY/${SAVE_FREQUENCY}/" | sed "s/CPP/${CPP}/" > ${NEW}/${DIR}/${MYMOL}_md${NOUT}.mdp
posreFlag="-r md0_success/${MYMOL}_md0_deshuffleddesorted.gro"
else
sed "s/TINIT/${TINIT}/" ${MYMOL}.mdp | sed "s/NSTEPS/${nsteps}/" | sed "s/SAVE_FREQUENCY/${SAVE_FREQUENCY}/" | sed "s/CPP/${CPP}/" > ${NEW}/${DIR}/${MYMOL}_md${NOUT}.mdp
posreFlag=""
fi
waitForExistNotEmpty 1 ${NEW}/${DIR}/${MYMOL}_md${NOUT}.mdp
if [ ! `tail -1 ${NEW}/${DIR}/${MYMOL}_md${NOUT}.mdp | awk '{print $1}'` = ";EOF" ]; then sleep ${TINY_SLEEP}; fi
if [ ! `tail -1 ${NEW}/${DIR}/${MYMOL}_md${NOUT}.mdp | awk '{print $1}'` = ";EOF" ]; then sleep ${SHORT_SLEEP}; fi
if [ ! `tail -1 ${NEW}/${DIR}/${MYMOL}_md${NOUT}.mdp | awk '{print $1}'` = ";EOF" ]; then sleep ${LONG_SLEEP}; fi

if((NIN!=0)); then
if((NOUT<=NJUMPS_POSRE || NEVER_USE_SORT_SHUFFLE==1 || MYNP==1)); then
# Can not do shuffle/sort with posre
${ED}/grompp -np ${MYNP} -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}.mdp ${posreFlag} -c ${PREV}/${MYMOL}_md${NIN}_deshuffleddesorted.gro -t ${PREV}/${MYMOL}_md${NIN}_deshuffleddesorted.trr -p ${MYMOL}.top -n ${MYMOL}.ndx -e ${PREV}/${MYMOL}_md${NIN}.edr -o ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr
rm mdout.mdp &
else
# Shuffle the .trr input file correctly. Assume that it is not currently shuffled
${ED}/grompp -np ${MYNP} -shuffle -sort -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}.mdp ${posreFlag} -c ${PREV}/${MYMOL}_md${NIN}_deshuffleddesorted.gro -p ${MYMOL}.top -n ${MYMOL}.ndx -o ${NEW}/${DIR}/${MYMOL}_md${NOUT}_a.tpr -deshuf ${NEW}/${DIR}/deshuffle_md${NOUT}_a.ndx
rm -f ${NEW}/${DIR}/deshuffle_md${NOUT}_a.ndx mdout.mdp &
echo System | ${ED}/editconf -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}_a.tpr -o ${NEW}/${DIR}/${MYMOL}_md${NOUT}_shuffledsortedInit_a.gro
# g_desort -f original shuffled will unshuffle, therefore g_desort -f shuffled original will REshuffle
${DD}/g_desort -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}_shuffledsortedInit_a.gro ${PREV}/${MYMOL}_md${NIN}_deshuffleddesorted.gro -o ${NEW}/${DIR}/reshuffleresort${MYMOL}_md${NOUT}_a.ndx -n 6
${ED}/trjconv -f ${PREV}/${MYMOL}_md${NIN}_deshuffleddesorted.trr -o ${PREV}/${MYMOL}_md${NIN}_reshuffleresort.trr -n ${NEW}/${DIR}/reshuffleresort${MYMOL}_md${NOUT}_a.ndx
# Create the run input file
${ED}/grompp -np ${MYNP} -shuffle -sort -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}.mdp ${posreFlag} -c ${PREV}/${MYMOL}_md${NIN}_deshuffleddesorted.gro -t ${PREV}/${MYMOL}_md${NIN}_reshuffleresort.trr -p ${MYMOL}.top -n ${MYMOL}.ndx -deshuf ${NEW}/${DIR}/deshuffle_md${NOUT}.ndx -e ${PREV}/${MYMOL}_md${NIN}.edr -o ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr
rm -f ${NEW}/${DIR}/deshuffle_md${NOUT}.ndx mdout.mdp &

# In the future: implement this check. Note that this will require rethinking the _a postfixes
# since the files without the _a postfixes are the ones that I actually use
#if((DOUBLE_CHECK_DESORT!=0)); then

# If a new reshuffle.ndx file differs then the run is invalid.
echo System | ${ED}/editconf -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr -o ${NEW}/${DIR}/${MYMOL}_md${NOUT}_shuffledsortedInit.gro
${DD}/g_desort -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}_shuffledsortedInit.gro ${PREV}/${MYMOL}_md${NIN}_deshuffleddesorted.gro -o ${NEW}/${DIR}/reshuffleresort${MYMOL}_md${NOUT}.ndx -n 6
case "$CLUSTER" in
*)
look=`diff -q ${NEW}/${DIR}/reshuffleresort${MYMOL}_md${NOUT}_a.ndx ${NEW}/${DIR}/reshuffleresort${MYMOL}_md${NOUT}.ndx`
;;
esac
if [ -n "$look" ]; then
echo There was a big problem. ${NEW}/${DIR}/reshuffleresort_md${NOUT}_a.ndx and ${NEW}/${DIR}/reshuffleresort_md${NOUT}.ndx are different.
mv ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr ${NEW}/${DIR}/${MYMOL}_md${NOUT}_notValid.tpr
fi

# End of the new test
#fi

#Create the deshuffle file to properly handle the next run
${DD}/g_desort -f ${PREV}/${MYMOL}_md${NIN}_deshuffleddesorted.gro ${NEW}/${DIR}/${MYMOL}_md${NOUT}_shuffledsortedInit.gro -o ${NEW}/${DIR}/deshuffledesort${MYMOL}_md${NOUT}.ndx -n 6
rm -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}_a.tpr ${NEW}/${DIR}/${MYMOL}_md${NOUT}_shuffledsortedInit_a.gro ${NEW}/${DIR}/reshuffleresort${MYMOL}_md${NOUT}_a.ndx ${PREV}/${MYMOL}_md${NIN}_reshuffleresort.trr ${NEW}/${DIR}/${MYMOL}_md${NOUT}_shuffledsortedInit.gro ${NEW}/${DIR}/reshuffleresort${MYMOL}_md${NOUT}.ndx mdout.mdp &
fi
else
if((NOUT<=NJUMPS_POSRE || NEVER_USE_SORT_SHUFFLE==1 || MYNP==1)); then
# Can not do shuffle/sort with posre
${ED}/grompp -np ${MYNP} -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}.mdp ${posreFlag} -c ${PREV}/${MYMOL}_md${NIN}_deshuffleddesorted.gro -p ${MYMOL}.top -n ${MYMOL}.ndx -o ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr
rm mdout.mdp &
else
${ED}/grompp -np ${MYNP} -shuffle -sort -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}.mdp ${posreFlag} -c ${PREV}/${MYMOL}_md${NIN}_deshuffleddesorted.gro -p ${MYMOL}.top -n ${MYMOL}.ndx -deshuf ${NEW}/${DIR}/deshuffle_md${NOUT}.ndx -o ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr
rm -f ${NEW}/${DIR}/deshuffle_md${NOUT}.ndx mdout.mdp &
echo System | ${ED}/editconf -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr -o ${NEW}/${DIR}/${MYMOL}_md${NOUT}_shuffledsortedInit.gro
${DD}/g_desort -f ${PREV}/${MYMOL}_md${NIN}_deshuffleddesorted.gro ${NEW}/${DIR}/${MYMOL}_md${NOUT}_shuffledsortedInit.gro -o ${NEW}/${DIR}/deshuffledesort${MYMOL}_md${NOUT}.ndx -n 6
rm -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}_shuffledsortedInit.gro &
fi
fi

if [ -e ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr ]; then
echo ${NOUT} > finished_grompp
waitForExistNotEmpty 12 finished_grompp ${NOUT}
rm -f finished_test
gromppProblemsInARow=0
else
# can not revert
let "gromppProblemsInARow=$gromppProblemsInARow+1"
if((gromppProblemsInARow>MAX_CONSECUTIVE_GROMPP_ERRORS)); then
touch ./DO_NOT_RUN
exit
fi
fi

fi

# MDRUN
if [ -e finished_grompp ]; then
NOUT=`cat finished_grompp`
TINIT=`cat finished_next_start_time`
DIR=md${NOUT}_running

# Reversion is important in cases where a crash or time overrun leads to loss of data in ${TMPDIR}
if [ ! -e ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr ]; then
sleep ${SHORT_SLEEP}
if [ ! -e ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr ]; then
rm -f finished_grompp
echo "ERROR error: finished_grompp existed for NOUT=$NOUT but ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr did not exist"
if((reverted==0)); then
echo " reverting to finished_test"
reverted=1
let "NIN=$NOUT-1"
echo "$NIN" > finished_test
if [ ! -s finished_test ]; then sleep ${SHORT_SLEEP}; fi
continue
else
touch ./DO_NOT_RUN
exit
fi
else
reverted=0;
fi
fi

if((MYNP==1)); then
returnValue=`${ED}/mdrun -deffnm ${NEW}/${DIR}/${MYMOL}_md${NOUT}`
else
case "$CLUSTER" in
*)
returnValue=`${mpirunProg} ${PED}/${mdrun_mpiProg} -np ${MYNP} -deffnm ${NEW}/${DIR}/${MYMOL}_md${NOUT}`
;;
esac
fi
if((returnValue!=0)); then
echo "ERROR error: mpirun for mdrun_mpi returned non-zero (${returnValue}). Exiting"
exit
fi
echo ${NOUT} > finished_mdrun
waitForExistNotEmpty 12 finished_mdrun ${NOUT}
rm -f finished_grompp
fi

# DESORT
if [ -e finished_mdrun ]; then
NOUT=`cat finished_mdrun`
TINIT=`cat finished_next_start_time`
DIR=md${NOUT}_running

# Reversion is important in cases where a crash or time overrun leads to loss of data in ${TMPDIR}
if [ ! -e ${NEW}/${DIR}/${MYMOL}_md${NOUT}.xtc ]; then
rm -f finished_mdrun
echo "ERROR error: finished_mdrun existed for NOUT=$NOUT but ${NEW}/${DIR}/${MYMOL}_md${NOUT}.xtc did not exist"
if((reverted==0)); then
echo " reverting to finished_test"
reverted=1
let "NIN=$NOUT-1"
echo "$NIN" > finished_test
waitForExistNotEmpty 12 finished_test ${NIN}
continue
else
touch ./DO_NOT_RUN
exit
fi
else
reverted=0;
fi

if((NOUT<=NJUMPS_POSRE || NEVER_USE_SORT_SHUFFLE==1 || MYNP==1)); then
# Can not do shuffle/sort with posre
mv ${NEW}/${DIR}/${MYMOL}_md${NOUT}.xtc ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.xtc
mv ${NEW}/${DIR}/${MYMOL}_md${NOUT}.trr ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.trr
mv ${NEW}/${DIR}/${MYMOL}_md${NOUT}.gro ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.gro
else
echo System | ${ED}/trjconv -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}.xtc -s ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr -n ${NEW}/${DIR}/deshuffledesort${MYMOL}_md${NOUT}.ndx -o ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.xtc &
echo System | ${ED}/trjconv -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}.trr -s ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr -n ${NEW}/${DIR}/deshuffledesort${MYMOL}_md${NOUT}.ndx -o ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.trr &
echo System | ${ED}/trjconv -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}.gro -s ${NEW}/${DIR}/${MYMOL}_md${NOUT}.tpr -n ${NEW}/${DIR}/deshuffledesort${MYMOL}_md${NOUT}.ndx -o ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.gro &
# Wait for all 3 desorts to finish
wait
fi

echo ${NOUT} > finished_desort
waitForExistNotEmpty 12 finished_desort ${NOUT}
rm -f finished_mdrun
fi

# TEST
if [ -e finished_desort ]; then
TINIT=`cat finished_next_start_time`
runHasNoErrors=1
NOUT=`cat finished_desort`
DIR=md${NOUT}_running

# Reversion is important in cases where a crash or time overrun leads to loss of data in ${TMPDIR}
if [ ! -e ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.xtc ]; then
rm -f finished_desort
echo "ERROR error: finished_desort existed for NOUT=$NOUT but ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.xtc did not exist"
if((reverted==0)); then
echo " reverting to finished_test"
reverted=1
let "NIN=$NOUT-1"
echo "$NIN" > finished_test
waitForExistNotEmpty 12 finished_test ${NIN}
continue
else
touch ./DO_NOT_RUN
exit
fi
else
reverted=0;
fi

${ED}/gmxcheck -f ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.xtc 2> ${NEW}/${DIR}/checkXTC
# Ensure no magic number error
magicNumberError=`grep Error ${NEW}/${DIR}/checkXTC | wc -l`
if((magicNumberError==1));then
mv ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.xtc ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted_magicNumberError.xtc
mv ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.trr ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted_magicNumberError.trr
mv ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.gro ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted_magicNumberError.gro
runHasNoErrors=0
fi
## Ensure expected number of frames is reached -- this is system and mdp file specific
numFramesXTC=`grep "^Step" ${NEW}/${DIR}/checkXTC | awk '{print $2}'`
if((numFramesXTC!=NFRAMES_IN_XTC)); then
mv ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.xtc ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted_incompleteFrames.xtc
mv ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.trr ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted_incompleteFrames.trr
mv ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted.gro ${NEW}/${DIR}/${MYMOL}_md${NOUT}_deshuffleddesorted_incompleteFrames.gro
runHasNoErrors=0
fi
if((runHasNoErrors)); then
mv ${NEW}/${DIR} md${NOUT}_success
waitForExistNotEmpty -1 ${NEW}/${DIR}
while [ -e ${NEW}/${DIR} ]; do
mv ${NEW}/${DIR} md${NOUT}_success
sleep ${LONG_SLEEP}
done

echo ${NOUT} > finished_test
nextTime=`echo "${TINIT}+${TJUMP}" | bc -l`
echo ${nextTime} > finished_next_start_time
waitForExistNotEmpty 12 finished_test ${NOUT}
waitForExistNotEmpty 12 finished_next_start_time ${nextTime}

# Now do some clean up operations
let "NCLEAN=$NOUT-2"
if [ -e md${NCLEAN}_success ]; then
if [ ! -e DATA ]; then
mkdir DATA
fi
if((NCLEAN!=0)); then
# Don't remove or modify the starting directory
mkdir DATA/md${NCLEAN}_success
mv md${NCLEAN}_success/${MYMOL}_md${NCLEAN}_deshuffleddesorted.xtc DATA/md${NCLEAN}_success
waitForExistNotEmpty -1 md${NCLEAN}_success/${MYMOL}_md${NCLEAN}_deshuffleddesorted.xtc
while [ -e md${NCLEAN}_success/${MYMOL}_md${NCLEAN}_deshuffleddesorted.xtc ]; do
mv md${NCLEAN}_success/${MYMOL}_md${NCLEAN}_deshuffleddesorted.xtc DATA/md${NCLEAN}_success
sleep ${LONG_SLEEP}
done
mv md${NCLEAN}_success/${MYMOL}_md${NCLEAN}.edr DATA/md${NCLEAN}_success
waitForExistNotEmpty -1 md${NCLEAN}_success/${MYMOL}_md${NCLEAN}.edr
while [ -e md${NCLEAN}_success/${MYMOL}_md${NCLEAN}.edr ]; do
mv md${NCLEAN}_success/${MYMOL}_md${NCLEAN}.edr DATA/md${NCLEAN}_success
sleep ${LONG_SLEEP}
done
rm -rf md${NCLEAN}_success &
fi
fi
else
# Send the run back to do the mdrun
for((i=1;i<=MAX_CONSECUTIVE_MDRUN_ERRORS;i++)); do
if [ ! -e md${NOUT}_failure${i} ]; then
break
fi
done
mv ${NEW}/${DIR} md${NOUT}_failure${i}
if((i>MAX_CONSECUTIVE_MDRUN_ERRORS)); then
echo "Too many failure for run $NOUT"
touch ./DO_NOT_RUN
exit
fi
# Send it back by setting as if grompp just finished
mkdir ${NEW}/${DIR}
mv md${NOUT}_failure${i}/${MYMOL}_md${NOUT}.tpr md${NOUT}_failure${i}/deshuffledesort${MYMOL}_md${NOUT}.ndx ${NEW}/${DIR}
echo ${NOUT} > finished_grompp
waitForExistNotEmpty 12 finished_grompp ${NOUT}
fi
rm -f finished_desort
fi

wait
done

wait


echo -n "TIMING TEST (end): "
date

恒定pH值模拟

首先,检查你的假设。理想的情况是有一个恒定pH算法来执行MD模拟。然而,传统的显式溶剂MD算法不能做到这一点,因为它们是恒定H+算法。此外,如果你使用MM力场,你将永远不会在你的溶质中得到质子化/去质子化态:自由的H+会四处游荡,而可滴定的位点会停留在它们最初的质子化状态;如果你想克服这个问题,你必须使用非MM哈密顿函数,它明确地允许质子转移。在任何情况下,不管用什么哈密顿量,H+的量永远不会改变。
如果你想做一个显式溶剂恒pH模拟,最简单的方法就是坚持传统的MD,只选择特定数量的H+,这是你的系统的“典型”。您将需要使用比MD模拟中通常使用的水多几个数量级的水。基本上,你可以忽略自由H+浓度在大多数pH值(即使pH=4),它的浓度是数量级以下的其他反离子(如,Na+,Cl-)。此外,如上所述,纯MM哈密顿函数将无法在可滴定位点之间移动质子(例如,与EVB不同),因此实际上需要选择分子中每个可滴定位点的质子化状态。
通过对起始结构执行标准pKa计算,可以很好地猜测起始配置;这些是标准的程序,主要基于连续体静电学计算质子自由能(如使用MEAD, UHBD, DelPhi, APBS)和蒙特卡罗执行质子状态采样(如使用REDTI, PETIT)。遗憾的是,随着MD模拟过程中溶质构象的变化,质子化状态可能会变得不充分,这是由于在许多情况下存在强质子-构象耦合。已经提出了几个解决方案,但大多数都或多或少是启发式的尝试。真正令人满意的解决方案是采用恒pH MD方法,最近几年已经有人提出了一些方法。这已经在GROMACS用户邮件列表中讨论过了。(请注意,Phil Hunenberger的方法存在一个严重的理论问题。)不幸的是,恒定pH方法是最近才出现的,并且还在开发或测试中。希望它们在不久的将来成为标准方法。也许有一天您可以在GROMACS .mdp文件中指定“pH = 7.0”,然后在MD运行期间查看质子状态的变化!在此之前,最好的解决方案可能是上面提到的一个:用标准的pKa计算进行初步的良好猜测,最后,使用MD运行时的快照检查它。
有一个恒定pH值的模拟模型,由Charlie Brooks和其他人开发,你可以在这个模型中模拟质子从一个侧链转移到另一个侧链。到目前为止,它们只在隐式溶剂中起作用(在CHARMM中)。
有其他几种恒pH MD方法,其中一些使用显式溶剂。Antonio Baptista的团队开发了一种基于随机质子化变化的恒pH MD方法(J. Chem. Phys. (2002) 117:4184)。虽然该方法使用波松-玻尔兹曼方法周期性地改变质子化状态,但在显式溶剂下进行了MM/MD模拟。
他们实际上已经使用GROMACS实现了这种随机恒pH MD方法(J. Phys. Chem. B (2006) 110:2927),基本上遵循一种走走停停的方法,使用bashawk脚本将GROMACS与MEAD (通过Don BashfordPB求解)和MCRP(用于蒙特卡罗采样质子状态的内部程序)进行接口。不幸的是,整个过程在某些地方仍然过于混乱和难以处理,这使得它不适合作为GROMACS的贡献提交,至少目前不适合。
Hunenberger(J. Chem. Phys. (2001) 114:9706)提出了另一种基于MM/MD的显式溶剂法,但其理论基础似乎是错误的 (J. Chem. Phys. (2002) 116:7766)。据我所知,目前提出的其他恒pH MD方法均采用隐式溶剂。McCammon(J. Comput. Chem. (2004) 25:2038)和Antosiewicz(Phys. Rev. E (2002) 66:051911)采用隐式溶剂随机方法, Baptista的小组还提出了一种分数电荷法(Proteins (1997)27:523),其中隐式溶剂用于计算速度;所有这些方法都依赖于某种简化的静电导向方法(波松-玻尔兹曼法、广义玻恩法等)来进行质子化计算。Brooks (Proteins (2004) 56:738)提出的方法也使用了隐式溶剂,采用了一种不同但理论上模糊的方法来包含质子效应。
讨论可以离解的水模型见这篇文章
说明:上面的内容来自GROMACS用户邮件列表上的各种电子邮件,但大部分来自Antonio Baptista。

扩散常数

使用GROMACS轨迹,您可以用Einstein) (g_msd)和Green-Kubo (g_velaccg_analyze)方程计算扩散常数。虽然g_msd是不言自明的,但是从版本4.0开始使用g_velacc就有点棘手了。
通过Green-Kubo方程计算粒子的扩散常数,需要执行以下命令:

1
g_velacc -acflen 1001 -nonormalize -mol -n atoms.ndx -s topol.tpr

这里atoms.ndx应该包含组成感兴趣的粒子的原子数。默认情况下,您将得到具有线性速度自相关函数(VACF)的vac.xvg
然后得到扩散常数类型:

1
g_analyze -f vac.xvg -integrate

这就得到了VACF的积分($nm^2/ps$)
最终,如果你从VACF中得到的积分值为$6 \times 10^{-5}nm^2/ps$,这个量乘以$10^{-6}$(同时转换$nm^2到m^2$,$ps到s$),并除以3(根据Green-Kubo方程),最后得到$2 \times 10^{-11} m^2/s$。这是扩散常数的值。

二面角PCA

要为二面角PCA创建索引文件,可以使用mk_angndx,也可以手动创建索引文件。现在您需要g_angle(可能是-oc-or)和g_covar的适当组合。
这儿有关于二面角PCA的一些讨论,并给出了相应的答复

背景

在GROMACS中,PCA的实现首先生成与所选角度匹配的降维轨迹文件,然后生成包含特征向量和特征值的伪轨迹文件。这可能不是最好的方法,但它是有效的。第1步和第2步进行降维,第3步和第4步生成一种坐标文件,其维数足够大,可以帮助gmx covargmx anaeig工作。

逐步的说明

  1. 使用mk_angndx或文本编辑器生成所选二面角中的原子index文件。
  2. 使用步骤1中的索引文件从轨迹中提取角度。例如dangle.ndx,在这样的命令中:
    1
    g_angle -f foo.xtc -s foo.tpr -n dangle.ndx -or dangle.trr -type dihedral

这里将输出dangle.trr,它包含简化到所选角度尺寸的仿真轨迹。

  1. 创建一个索引文件,命名为covar.ndx,其中必然包含一组原子1到整数(2*N/3),其中N是二面角的个数。例如,对于一个5个二面角的肽,covar.ndx的内容应该是:

    1
    2
    [ foo ]
    1 2 3 4
  2. 使用索引文件创建包含整数(2*N/3)原子的.gro文件。盒子大小和原子坐标无关紧要。例如:

    1
    trjconv -s foo.tpr -f dangle.trr -o resized.gro -n covar.ndx -e
  3. 使用类似的命令执行对角化:

    1
    g_covar -f dangle.trr -n covar.ndx -ascii -xpm -nofit -nomwa -noref -nopbc -s resized.gro

这里将输出eigenval.xvg, eigenvec.trr, covar.log, covar.dat, 和 covar.xpm文件。

  1. 要沿着一个特征向量得到PMF,可以使用这样的命令:

    1
    g_anaeig -v eigenvec.trr -f dangle.trr -s resized.gro -first X -last X -proj proj-1

    其中X为特征向量的序号。要可视化结果,请使用:

    1
    xmgrace proj-X.xvg
  2. 要得到沿两个特征向量投影的自由能形貌,可以使用类似这样的命令:

    1
    g_anaeig -v eigenvec.trr -f dangle.trr -noxvgr -s resized.gro -first X -last Y -2d 2dproj_X_Y.xvg

其中XY是特征向量的序号。

  1. 为了将这些数据转换为自由能形貌图,我们将二维景观与网格划分,并计算每个网格内的构象数。这里发布了一个用于此目的的awk脚本g_sham程序(指定-notime选项)还可以用来生成自由能表面的二维图。

二面角限制

在GROMACS 4.6之后的版本中,二面角约束完全在拓扑中指定。以前使用的mdp设置已被删除。
在GROMACS 3.3.1版本中,二面角约束的原子规范是拓扑文件的一部分,力常数等在.mdp文件中指定。
假设在你的peptide.itp中原子的序号是这样的:

1
2
3
4
5
C' (n-1) = 5
N (n) = 7
CA (n) = 9
C' (n) = 15
N (n+1) = 17

然后在.top文件中使用这样的部分,其中180被实际的二面角值(以角度为单位)替换,您想要限制原子的二面角值:

1
2
3
4
5
6
7
8
9
10
#include "peptide.itp"
[ dihedral_restraints ]
; ai aj ak al type label phi dphi kfac power
; phi C'(n-1) - N - CA - C'
5 7 9 15 1 1 180 0 1 2
; psi N - CA - C' - N(n+1)
7 9 15 17 1 1 180 0 1 2

#include "tip3p.itp"
;etc...

确保dihedral_restraints在包含蛋白质拓扑之后。如果您的.top文件中直接包含了蛋白质拓扑结构,那么只需在蛋白质列表之后,但在任何提到不是那个蛋白质分子的东西之前,插入dihedral_restraints即可。
手册第4.3.3节描述了二面角限制的实现。[dihedral_constraints]指令中指定的参数如下:

  • type: 唯一的值1
  • label: 未使用,并已从代码中删除,以供以后的版本使用
  • phi: 手册中4.74和4.75方程中的$\varphi_0$的值
  • dphi:手册中4.75方程中的$\Delta \varphi$的值
  • kfac: 类似于fac在实现距离限制时,力常数(在.mdp文件中指定,见下)乘以的因子。这样,即使只提供一个dihre_fc值,也可以用不同的力常数维护不同的限制
  • power:未使用,并已从代码中删除,以供以后的版本使用

注意:

  1. 在3.3版中,二面角限制要么工作方式不同,要么被破坏;使用GROMACS版本3.3.1(Mobley, Chodera and Dill, J. Chem. Phys. 125, 084902 (2006), page 084902-6, section III D)
  2. 手册中对这种二面角限制在近180度使用时是否稳定还不是很清楚。Chris Neale发现,在包括180度在内的整个二面角范围内,一切似乎都表现得像预期的那样正常。但是,必须避免实际的二面角距离限制二面角接近180度的情况。
    2b). 不是每个人都同意上面的第2部分。看这里。因此,请用户做一些测试来得出自己的结论。
  3. 特别注意单位。请注意在拓扑文件的输入中,角度以度数表示,力常数以kJ mol-1 rad-2表示。

距离限制

通过在拓扑中设置[distance_constraints],可以将相同[molecular uletype]内的两个原子之间的距离限制在规定的距离内。手册第4.3.4节描述了该算法的实现和底层方程。示例[distance_constraints]部分(取自手册)可能如下所示:

1
2
3
4
5
6
[ distance_restraints ]
; ai aj type index type’ low up1 up2 fac
10 16 1 0 1 0.0 0.3 0.4 1.0
10 28 1 1 1 0.0 0.3 0.4 1.0
10 46 1 1 1 0.0 0.3 0.4 1.0
16 22 1 2 1 0.0 0.3 0.4 2.5

上面的指令指定了原子对(ai和aj)之间的距离限制。使用的术语如下:

  • type: 唯一值1
  • index: 具有相同index值的多个限制可以一起控制。上面,第二个和第三个限制一起考虑,而第一个和第四个限制分开的。
  • type’: 值可以是1或2;使用2意味着限制将不是时间或传感器平均(用于限制氢键)
  • low: 手册公式4.76和图4.13中r0的值
  • up1: 手册公式4.76和图4.13中r1的值
  • up2: 手册公式4.76和图4.13中r2的值
  • fac: 乘以力常数(.mdp文件中的disre_fc)的因子。这样,即使只给出一个disre_fc值,不同的限制也可以由不同的力常数来维持

重启

通常

为了精确地重新启动模拟,必须保留系统的所有状态变量。在实际中,这意味着高精度地保存坐标、速度和能量分量。下面的大部分讨论都是关于如何在GROMACS 3.x中重启崩溃模拟。GROMACS 4.x要简单得多,可以先处理。

版本4.x(和以后)

随着检查点的引入,下面给出了3.X的说明部分过时了。它们应该仍然可以工作,但现在,中断长时间的mdrun或从崩溃中恢复的简单情况更容易了。如果模拟崩溃,请利用状态。所写的state.cpt文件,它包含了继续模拟所需的所有信息。为了从模拟停止的地方开始,只需对mdrun使用-cpi-append选项。注意,-append是4.5中的默认值。
在做任何事情之前,先备份您的文件!然后使用:

1
mdrun -s topol.tpr -cpi state.cpt

这个命令将从编写检查点文件时开始继续模拟。除非使用-noappend选项,否则剩余的输出将附加到现有文件(能量、轨迹、日志等)。由于检查点文件包含所有输出文件的校验和,因此这个过程非常简单,不会覆盖或附加到以前的部分输出文件或修改后的输出文件。
在4.0版本中,将新的输出文件附加到旧的输出文件并不是默认值,所以您应该添加-append选项:

1
mdrun -s topol.tpr -cpi state.cpt -append

由于没有校验和,所以需要确保附加到正确的、未修改的文件中。
如果您需要生成一个新的.tpr(例如,因为你必须扩展模拟超出了你计划的步骤,或改变.mdp文件选项)然后你可以按照说明,但不需要供应超过.tpr文件tpbconv只要你旧的检查点文件后续mdrun供应。如果更改了integrator或系综,应该只将检查点文件传递给tpbconv,而不是mdrun,因为状态可能会更改,因此无法添加输出文件。
注意,mdrun将写入state.cpt state_prev.cpt文件。从他们的时间戳可以看出,其中一个大约是在检查点间隔之前编写的(默认情况下是15分钟)。或者您可以使用gmx check查看其中的内容。

版本3.x

在GROMACS版本3.x中,重启模拟有两种方法。一个使用grompp,需要.mdp文件和坐标文件(仅用于原子名称),另一个使用tpbconv,需要.tpr文件。要实现保持模拟集成的重新启动,这两种方法都需要一个完全精确的轨迹(例如.trr文件),当使用压力耦合或Nosé-Hoover温度耦合时,都需要一个.edr文件。注意,当恒定体积系综切换到恒定压力系综时,您不应该提供.edr文件,因为旧的文件中没有相关的术语,而命令行中文件的存在意味着GROMACS需要这些术语。
注意,为了重新启动,需要同时具有坐标和速度(可能还有能量)的轨迹帧,所以在选择.mdp文件中的nstxoutnstvoutnstenergy时要考虑这一点。grompptpbconv都允许您使用一个时间选项来选择将用于重启的所提供的轨迹点。

使用grompp

使用grompp,您需要重新考虑.mdp文件选项的内容,比如tinitinit_stepnsteps。对于一个精确的延续,这些选项应该分别是模拟的第一部分的开始时间、模拟要从哪个步骤开始重新启动以及在这个段中要执行的步骤数。这些选择影响新模拟的起始值,但不影响起始配置的帧的选择,起始配置将是最后一个可用的帧(默认情况下),或者使用-time选项选择。确保在重新启动时设置gen_vel = no,否则将忽略输入轨迹中的速度。在适当的地方使用unconstrained_start。阅读GROMACS手册第7.3节,如果有必要,问问自己为什么在运行第一个模拟之前不这样做。如果可能的话,更简单的方法是使用tpbconv重启,但是如果需要更改系综或输出选项等,那么grompp是唯一的方法。

使用tpbconv

使用tpbconv时,使用-t-e选项来提供模拟结束时的信息,并使用-extend-until选项。阅读手册页!tpbconv将保留原始运行输入文件中的所有参数,以及用于模拟的原始grompp命令行上请求的处理器数量。
注意,-extend选项将扩展剩余的运行时,因此新的模拟将从您使用其他选项选择的位置开始,一直到预期的模拟结束,然后根据您选择的数量进行扩展。
对于那些担心在初始.mdp文件中使用gen_vel覆盖平衡速度的人,grompp在生成时只读取一次该参数,以后再也不会读取。
对于那些担心他们的输入.tpr文件有速度和输入轨迹文件的人,不用担心…tpbconv创建了你跑步的延续,而不是一些奇怪的混合了新位置和旧速度。

崩溃后

你需要从一个.trr坐标系中重新开始,包括位置和速度,如果必要的话,还需要同时从一个.edr坐标系中重新开始。这是否可能取决于您对上述.mdp文件参数的选择。您可以使用gmx check查看哪一帧提供了最新的良好重启。

为崩溃做准备

然您不希望有一个参数,但是选择您的参数以便为您的计算机时间和磁盘空间获得最佳值是有意义的。第一件要确保的事是每次你有完全精确的位置和速度你都有能量。后一种组合通常存储起来相当昂贵,因此如果您需要频繁的位置进行分析,您应该相应地使用nstxtcout。如果崩溃的概率很低,那么每隔几个小时模拟一次全精度帧似乎是存储成本和重新运行模拟直到上次崩溃的成本之间的合理折衷。如果经常崩溃,您可能希望经常编写完全精确的帧,并在事后使用trjconv来增加.trr文件的间隔,并降低长期存储成本。当然,您还应该在集群管理员下点火!
在GROMACS 4.x中,检查点处理这个问题,但是您可能希望在命令行上选择检查点间隔到mdrun

另请参考:

扩展模拟

本质动力学

本质动力学提取蛋白质的相关运动,以了解对蛋白质活性最基本的运动。蛋白质相关运动的分析可以通过g_covarg_anaeig来完成。
资源:

扩展模拟

可以扩展或继续已经完成的模拟,甚至那些已经崩溃的模拟(参见重新启动)。这实际上是一个很好的技术,用于处理较长的模拟,以减少由于计算机崩溃而造成的时间损失。只有从检查点文件或最后一个完全精确的坐标和速度可用的点开始,才能无缝地继续。尽管一些坐标文件格式可能包含速度(例如gro格式),但对于精确重启来说,这些速度不够精确。因此,您必须确保以足够短的时间(每天几次?)编写完全精确的数据,以避免由于崩溃而损失太多的时间。如何实现这一点取决于您正在使用的GROMACS版本。

版本4和以后

已经终止但尚未完成的模拟(例如,由于队列限制、电源故障或使用mdrun-maxh选项)可以继续进行,而不需要使用tpbconv(从5.0版开始称为gmx convert-tpr)。在这种情况下,您可能希望或不希望在mdrun命令中使用-append。否则,可以使用tpbconvmdrun和检查点文件(.cpt)扩展已经完成的模拟。首先,必须在.tpr文件中更改步骤的数量或时间,然后使用mdrun从最后一个检查点继续模拟。这将生成一个与连续运行时“相同”的模拟(但更多讨论请参见重现性)。

1
2
tpbconv -s previous.tpr -extend timetoextendby -o next.tpr
mdrun -s next.tpr -cpi previous.cpt

您可能希望使用mdrun-append选项将新输出追加到旧文件。注意,只有当用户没有修改旧的输出文件时,这才会工作。附加是4.5版本的默认行为。
如果您希望在可管理的部分中运行长时间的模拟时更改默认文件名,那么循环运行以下命令将会有效:

1
2
3
mdrun -deffnm ${name} -cpi ${name} -append
tpbconv -s ${name} -o ${name}_new.tpr -extend ${time}
mv ${name}_new.tpr ${name}.tpr

如果您觉得需要归档检查点并运行输入文件,那么您也可以这样做。如果使用-noappend,那么mdrun将根据您的名称向一系列文件添加数字后缀,就像mdrun -h中描述的那样。
在队列系统中运行时,设置使用grompptpbconv进行总体模拟所需的步骤数,并使用mdrun-maxh选项在队列时间结束之前优雅地停止模拟,这非常有用。使用此过程,您可以简单地通过让mdrun读取检查点文件来继续,而不需要其他工具。但是,如果队列系统允许作业挂起,-maxh机制将不知道挂起所花的时间,您可能会模拟比预期更短的挂起时间。
使用tpbconv-until-nsteps选项也可以延长时间。
模拟可以在没有检查点文件的情况下继续进行,检查点文件将是非二进制相同的,并且会有一些小错误,在大多数情况下,这些错误可以忽略不计。误差产生的原因是轨道和能量文件没有存储所有的状态变量的温控和气压。如果是这种情况,您必须使用下面的3.3.3版本过程。
更改.mdp文件选项
如果您希望/需要更改.mdp文件选项,那么两者都可以:

1
2
grompp -f new.mdp -c old.tpr -o new.tpr -t old.cpt
mdrun -s new.tpr

或者

1
2
grompp -f new.mdp -c old.tpr -o new.tpr
mdrun -s new.tpr -cpi old.cpt

应该工作。在GROMACS 4中,前者是必要的。如果热力学整体改变了。(有人说:“如果你老了。”cpt用于已经完成的运行,然后在grompp之后和mdrun之前使用tpbconv -extend。运行结束由.cpt.tpr上下文中的内容来判断。因此,如果后者被更改,那么运行就没有完成。)

版本3.3.3和以前

要继续进模拟,坐标和速度需要完全精确,这意味着必须使用.trr轨迹文件;xtc格式的轨迹是不够的,因为它们的精度降低了,并且没有速度信息。坐标和速度写入.trr文件的频率由.mdp文件中的nstxoutnstvout设置。该过程由tpbconv按以下方式处理:

1
tpbconv -s previous.tpr -f previous.trr -e previous.edr -o next.tpr -extend timetoextend

然后将生成的.tpr文件返回到mdrun。注意,在输入中必须有一个坐标系,在这个坐标系中,你写出了所有三个位置、速度和能量,这样tpbconv就可以把它们组合成一个新的起点。

精确 vs 二进制相同的延续

如果你有一台精度无限的电脑,或者你用手把时间离散化的运动方程积分,精确的延拓会得到相同的结果。但是由于计算机的精度有限(单精度或双精度),MD是混沌的,如果将一个位(最不重要的位)四舍五入,轨迹就会非常迅速地发散。样的轨迹都是同样有效的,但却截然不同。继续使用检查点文件,使用相同编译器编译的相同代码,并使用相同数量的处理器在相同的计算机体系结构上运行,将导致二进制相同的结果,除非使用计时代码。定时代码由动态负载平衡使用,(默认情况下)由FFTW FFT库使用。
大多数情况下,你并不关心轨迹不是二进制的,当你停止它并继续它,而不是让它在一次运行中运行。最常见的情况是,当您在第567894步发生崩溃时,您希望重新生成它,以跟踪它是否是系统的问题或bug。mdrun有一个选项-reprod(试图)强制执行二进制相同的模拟。注意,在某些情况下,这可能会消耗大量性能。

自由能计算

构建一个Linux集群

为GROMACS和其他程序构建自己的Linux集群(幻灯片视频)(Erik Lindahl)。

多链

为了生成一个拓扑来模拟多个不同的分子(即PDB-speak中的链),我们首先区分两种情况:拓扑相同的和拓扑不相同的。
相同的链
这是一个简单的例子。取一个只有一条链的结构文件,并在上面使用pdb2gmx。(注意,正常情况下,这个结构不需要使用PDB格式。)仔细查看结果.top文件,在底部您将看到如下部分:

1
2
[ molecules ] 
Protein

早些时候,这个分子被命名为“Protein”(或者别的什么)。你可以把计数器从1增加到你需要的值。然后,使用这个.top文件和另一个坐标文件(其链数对应于修改后的.top文件)作为grompp的输入。
您可以使用editconf创建链的平移副本,以便稍后在单个文件中进行连接,或者使用genconf成整个框的复制。
不同的链
这是棘手的。可以使用坐标文件格式来标识不同的链,比如PDB格式,它允许在第22列中分配一个字母(链标识符)。pdb2gmx程序可以识别不同的链标识符,并将它们写成不同的分子(从而分离拓扑)。例如,A链将对应一个名为topol_A.itp的拓扑,B链是topol_B.itp等。注意,在生成的.top文件中,每个链都将被写成一个单独的拓扑结构,即使这些链实际上表示相同的分子(即,均二聚体蛋白)。
在GROMACS 4.5中,pdb2gmx-chainsep参数允许用户控制PDB TER记录是否也应该考虑将一条链拆分为一个新的[moleculetype]部分。
需要分开的链吗?
注意,为了在GROMACS中进行模拟或运行分析,不需要使用单独的链标识符将链分隔成不同的分子。GROMACS并不真正关心[ moleculetype ]部分中包含什么,只是你不能在不同的这个部分之间有键合的相互作用。您可以创建一个索引文件,其中包含根据其编号残数指定链的残数。例如,如果你有一个二聚体蛋白,每个链由200个氨基酸组成,A链可能对应残基1-200,B链 201-400。因此,在make_ndx提示符处,您可以输入:

1
r 1-200

为第一条链,

1
r 201-400

为第二条链。这样,您就可以对每个链分别进行分析,而不用担心链标识符。

多拓扑项

在某些情况下,可能需要生成比所提供的标准形式更为复杂的结合项。这可以通过组合多项来实现,加起来就是最后需要的结合项。例如,当二面角势函数遵循Ryckaert-Bellemans函数,但在GROMACS中有超出标准n=5的附加项时,可以使用这种技术。
[bonds]、[pair]、[angle][dihedrals]部分中的拓扑文件项都是可添加的,每个附加项都添加到前一个节中。在这里复制一个项将使得到的势能项加倍(没有警告!),所以要小心修改。
要将在你 ffbonded.itp中定义成键类型的机制(例如,[bondtypes], [angletypes]和[dihedraltypes])混淆。这些指令中的多个项将覆盖以前的任何实例,导致grompp失败并发出警告。这个规则的例外是CHARMM力场。

简正模式分析

GROMACS简正过程
在简正模式分析(NMA)中,计算是为了研究大分子(通常是蛋白质)潜在的大规模功能运动。计算需要非常彻底的能量最小化和计算Hessian矩阵。
一般来说,对于最终阶段的能量最小化(使用L-BFGS)和实际的简正模式分析,都可以使用以下设置:

  • 双精度的GROMACS
  • 改变库仑和范德瓦尔斯相互作用;截断例如在1.0 nm,从0.8 nm变化(或者从0 nm
  • rlist1.21.3

在创建简正模式运行输入文件时,必须使用grompp-t选项,以便读取完全精确的二进制坐标,而不是三位小数的.gro文件。

资源

pKa计算

高聚物

GROMACS非常适合运行聚合材料的模拟。gmx-users列表中经常被问到的一个问题是,如何为包含多个重复单元的聚合物创建拓扑结构。可能最简单的方法是为所需的力场创建.rtp条目,并使用pdb2gmx创建拓扑。要做到这一点,(至少)需要定义三种残基:

  • 1.一个“开始”的残基,定义了链的“开始”
  • 2.重复(内部)残基
  • 3.一个“结束”的残基,定义链的“结束”

关于如何实现这样的条目的例子可以在这里找到,例如在OPLS-AA力场下的聚乙烯。
还有一个在线服务器允许用户构建某些类型的聚合物,这并不能解决构建拓扑结构的问题,但在生成材料的初始坐标方面肯定是有用的。

位置限制

位置限制算法用于将粒子/原子限制在一个固定的参考位置。因此,原子可以移动,但这样做会消耗大量的能量。
使用位置限制的原因:

  • 避免体系中关键部分的剧烈重排,如在添加溶剂后平衡时限制蛋白质。
  • 用受限制的粒子壳将被研究的区域与由于粒子缺失而缺乏适当相互作用的外部区域隔离(在GROMACS中没有实现)。

位置限制部分可以添加到拓扑文件中,并在需要时使用inclde文件机制激活。如果使用pdb2gmx, -i选项会让它写入非氢原子的posre.itp文件。然后将其包含在拓扑文件中。通过设置.mdp文件中的define = -DPOSRES选项,可以控制是否实际读取限制文件。

位置限制必须物理地位于要约束的[ moleculetype ]。在.top文件中不能有全局位置限制部分。查看相关的错误信息

可以使用genrestr实用程序(GROMACS 3.3.x和早些时候版中的genpr)创建特殊位置约束.itp文件)。例如,如果您希望仅对蛋白质的主链原子施加位置限制,请使用genrestr并选择“Backbone”作为输出组,然后在拓扑中写上#include "backbone_posre.itp"

例子
位置限制要在拓扑文件(.top或者.itp):

  • 原子/粒子的识别要应用于([ moleculetype ]而不是坐标文件)
  • 函数类型
  • 每个维度上的力常数(x, y, z)
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    [ moleculetype]
    This_one

    [ position_restraints ]
    ; ai funct fcx fcy fcz
    1 1 1000 1000 1000 ; restrains to a point
    2 1 1000 0 1000 ; restrains to a line (y-axis)
    2 1 1000 0 0 ; restrains to a plane (y-z-plane)

    [ moleculetype ]
    That_one

    [ position_restraints ]
    ; ai funct fcx fcy fcz
    1 1 1000 1000 1000 ; restrains to a point
    2 1 1000 0 1000 ; restrains to a line (y-axis)
    2 1 1000 0 0 ; restrains to a plane (y-z-plane)
    ; Note that the atom indices ai are relative to the [moleculetype], not the whole coordinate file

资源

  • 参见GROMACS手册第4.3.1节,了解用于在GROMACS内实现位置限制的公式的详细信息。
  • 有关在分子拓扑文件中如何实现位置约束的详细信息,请参阅GROMACS手册第5.6.1节。