GROMACS 力场分布分析（FDA）安装版下载与获取指南

原文链接：https://github.com/HITS-MBM/gromacs-fda

开发者

Mark Abraham <mark.j.abraham@gmail.com>

David van der Spoel <spoel@gromacs.org>

Berk Hess <hess@gromacs.org>

Erik Lindahl <lindahl@gromacs.org>

背景

长期以来，机械工程师一直使用内力和应力分布计算来设计坚固的结构，从高耸的摩天大楼和数公里长的桥梁，到一级方程式赛车和飞机机翼。同样的分析，但应用于原子和分子水平，可以揭示(生物)分子的机械稳定性，变构机制和细胞中的信号传播或机械诱导的过程，如血液凝固。

分子动力学(MD)模拟将牛顿运动方程应用于原子系统。在每个积分时间步，原子通过成对力相互作用。所有作用在原子上的成对力都被加起来，由此产生的力决定了原子在下一个时间步长的位移。这些位移通常用于分析MD模拟。

对单个原子的成对力的分析，我们称之为力分布分析，可以提供一个(生物)分子系统更敏感的观点，无论是在稳定状态还是在状态之间的过渡。这样的分析将能够解释，例如，像在一个原子水平上的牛顿摇篮功能，其中施加在原子上的力并不直接转化为所有球体的位移，而只是外层的球体的位移。类似地，原子可以在没有任何显著的原子位移的情况下传播力，例如通过蛋白质的刚性核。

引文

如果您使用FDA代码，请引用:http://dx.doi.org/10.1186/2046-1682-6-5

安装

GROMACS- fda基于官方GROMACS。实际分子动力学模拟将使用官方GROMACS进行，而力分布分析将使用GROMACS- fda进行。为了避免不同版本的GROMACS和GROMACS- fda之间的不一致，至关重要的是GROMACS版本号相同。请在https://github.com/HITS-MBM/gromacs-fda/releases上找到可用的FDA发布

安装gromacs

GROMACS的默认安装可以通过以下步骤构建:

git clone https://github.com/gromacs/gromacs.git -b <gromacs-version>

cd gromacs

mkdir build

cd build

cmake -DCMAKE_INSTALL_PREFIX=<installation path> ..

make -j <number of cores>

make check

make install

有关安装过程的详细信息，请参阅GROMACS官方文档。

安装GROMACS-FDA

GROMACS-FDA也可以在GitHub上使用，可以通过以下步骤构建:

git clone https://github.com/HITS-MBM/gromacs-fda.git -b

<gromacs-version>-fda<fda-version>cd gromacs-fda

mkdir buildcd build

cmake -DCMAKE_INSTALL_PREFIX=<installation path> -DGMX_BUILD_FDA=ON

-DGMX_DEFAULT_SUFFIX=OFF -DGMX_BINARY_SUFFIX=_fda -DGMX_SIMD=NONE

-DGMX_BUILD_UNITTESTS=ON -DGMX_GPU=OFF ..

make -j <number of cores>

make check

make install

-DGMX_DOUBLE=ON -DGMX_BINARY_SUFFIX=_fda_d

安装VMD可视化插件

安装是通过将归档文件解压缩到VMD插件目录来完成的;通常是plugins/noarch/tcl。此操作将创建两个名为pf_loaduser1.0和pf_draw1.0的目录，其中包含相应的Tcl文件。要测试安装，然后在控制台或Tk控制台类型中启动VMD:

package require pf_loaduser package require pf_draw

规则

残基的操作

对于基于残基的操作，依次进行以下步骤:

1.根据残基对应表建立原子

2.如果一个原子对是为两两力计算所选择的原子对的一部分，则从对应表中获得两个原子的剩余数

3.如果两个剩余数相等(来自同一剩余的两个不同的原子)，则什么也不做，因为计算剩余与自身的相互作用没有意义，代码会带着不同的原子对回到步骤2

4.将原子矢量力加到残基对对应的矢量力上

最后，每个残基对都有一个矢量力，计算公式为:其中i是残基ri的一个原子，j是残基rj的一个原子，ri和rj是不同的。

一旦计算出来，基于残基的成对力就可以像基于原子的那样写入输出文件——作为矢量或标量值——或者可以进一步用于计算每个残基的和/平均值/最小/最大值。

残基(重新)编号

对于GROMACS，在实际的MD运行过程中，残基是次要的——它们只在写出结构文件(例如PDB)时有用。原子和残基序号的编号之间的对应关系是基于结构文件中的信息进行的，其中残基序号编号可能不是从1开始的，也可能不是连续的——为了处理所有这些情况，默认情况下，GROMACS不尝试重新编号残基序号，并回写输入结构文件中分配给原子的任何残基序号。

PF2代码使用这些残基数来识别哪些原子属于哪些残基，并输出基于残基的数据。处理不从1开始或不连续的剩余数很容易-只需填充从1开始的缺失数，并使其输出包含零。但是残基数也可以出现多次，例如，如果结构文件包含几个链(比如蛋白质，配体和溶剂)，并且在每个结构文件中，残基开始以1编号，或者当有许多残基(超过9999，这是PDB文件中的最大数量)并且它们的数量围绕。为了解决这个问题，GROMACS可以对残基进行重新编号，从1开始以递增的数字对残基进行编号。PF2代码可以检测残基序号冲突，默认情况下，自动切换到重新编号残基序号;这种行为可以通过.pfi文件选项来控制。

对残基重新编号既会影响残基数据在内存中的存储方式，也会影响它在输出文件中的写入方式。是否在PF2代码中进行了重新编号的信息(无论是自动的还是通过.pfi文件选项强制执行的)必须传递给VMD插件;不这样做将导致错误的陈述!

成对的力分解为3和4体势

对于只涉及两个原子的势(键，库仑，伦纳德-琼斯)，对偶力与原子力相同，直接从势计算出来。对于涉及两个以上原子的势(角，二面角-正常和不正常，以及交叉键-键和交叉键角)，原子力需要分解以反映每两个原子之间的相互作用。由这种分解产生的成对力的矢量和应该等于原子力。PF2代码中实现的分解详见附录A。

执行

其他mdrun选项

FDA为gmx_fda mdrun添加了几个命令行选项:

Ø-pfi input.pfi – 指定一个输入文件，该文件控制代码应该如何运行

Ø-pfn index.ndx –索引文件是否包含一组或多组原子，这些原子应该用于成对力计算

Ø-pfa output.pfa – 输出文件，包含原子的成对力;内容和格式取决于.pfi文件中的选择

Ø-pfr output.pfr – 输出文件，包含残基的成对力;内容和格式取决于.pfi文件中的选择

Ø-psa output.psa – 输出文件，包含原子的准时应力;内容和格式取决于.pfi文件中的选择

Ø-psr output.psr – 输出文件，包含残基的准时应力;内容和格式取决于.pfi文件中的选择

Ø-vsa output.vsa – 输出文件，包含原子的虚拟应力;内容和格式取决于.pfi文件中的选择

Ø-vma output.vma – 输出文件，包含原子的冯米塞斯维里应力;内容和格式取决于.pfi文件中的选择

Gmx_fda mdrun应该总是与-rerun选项一起运行，要求事先生成初始轨迹;除了上面描述的PF2特定选项外，传递给gmx_fda mdrun的选项应该与生成原始轨迹时使用的选项相同，但有一个例外:应该强制执行单个CPU运行，因为PF2代码(尚未)是并行的。请注意，默认情况下，非mpi版本的GROMACS 4.5以线程方式运行，但仍然被视为并行运行，因此使用'-nt 1'来强制非并行运行:

gmx_fda mdrun -nt 1 -rerun traj.trr -pfi input.pfi -pfn index.ndx -pfa

output.pfa -pfr output.pfr <other options>

输入文件(.pfi)

Onepair-描述如何处理同一对原子I和j之间的成对力。通常，一个力场可以允许几种类型的非键相互作用或几种类型的键相互作用在同一对原子之间。例如，两个距离较远的原子可以通过库仑相互作用和伦纳德-琼斯相互作用相互作用，而两个成键原子可以通过一个键、几个角和几个二面角相互作用相互作用。

这个关键字有2种可能的值:

summed-对同一对原子之间的所有成对力求和;如果原子I和j相互作用，代码只存储它们之间成对作用力的一个值

detailed-在同一对原子之间的所有相同类型的成对力(例如角相互作用)被单独存储;如果原子I和j相互作用，代码存储它们之间至少一个，也可能是多个成对力，对应于它们之间存在的每种相互作用。

与求和选项相比，使用详细选项会导致所需内存量的显著增加。detail应该只在极少数情况下使用，如力场开发/调试，在这种情况下，应该获得由潜在类型分开的成对力的同时输出。

对默认值求和。例如，以下选项集将计算原子i和j之间所有相同类型的相互作用(例如二面角)的总和:

onepair = detailed

type = dihedral

Group1和group2 -命名将用于成对力计算的组;如果不存在，则默认为“蛋白质”。两组是相同的，在内部治疗上没有区别。组名应该存在于使用-pfn命令行选项指定的索引文件中;为此目的，没有必要使用与原始运行的grompp步骤不同的单独索引文件——只有在group1和group2中命名的组将被加载，其他组将被忽略。这两个组定义了哪些原子用于双力计算:只考虑一个原子属于组1而另一个原子属于组2的相互作用。两个基团中同时有一个或多个原子是可以的。在索引文件中只定义一个组，并为group1和group2指定这个组的名称也是可以的——这相当于最初的PF实现[1]，其中只定义了一个组。有两个基团允许更大的灵活性和效率:想象一个配体和蛋白质之间的相互作用的情况-组1可以被定义为包含蛋白质，组2可以被定义为包含配体并且成对相互作用将只计算它们之间的相互作用;在最初的PF实现中，只使用一个基团将意味着包括蛋白质和配体在该组中，并且获得可能无用的蛋白质-蛋白质和配体-配体相互作用。

Group1和group2总是包含原子索引，即使在只请求基于残基的计算时也是如此。指定一个或两个基团中残留的所有原子仍然是用户的责任。如果仅指定残基的某些原子，并且要求进行基于残基的计算，则与该残基相对应的结果将仅包含涉及指定原子的相互作用。代码不会自动填充缺失的原子以获得完整的残基。这样做是为了灵活性:例如，如果只对基于残基的主链相互作用感兴趣，则索引文件应该只包含主链原子，并且应该请求基于残基的计算。

type-指定感兴趣的交互类型。它可以有一个或多个值:

bond angle dihedral polar coulomb lj nb14 bonded nonbonded all

nb14定义为LJ14和Coulomb14相互作用的组合。LJ14和Coulomb14不能单独获得。

bond定义为键、角和二面体的组合。

nonbonded定义为库仑、lj和nb14的组合。

all都被定义为键、极性和非键的组合。

如果未指定类型，则假定为'type = all'。

GROMACS包含几个函数，用于计算键、角和二面角的势和力。使用哪一个取决于力场、.mdp文件中的选择(例如“morse=yes”)和用户对拓扑文件的修改(例如添加额外的谐波电位)。PF2代码根据所涉及的粒子数量来处理它们:所有涉及2个粒子的键相互作用都被认为是键，所有涉及3个粒子的键相互作用都被认为是角，所有涉及4个粒子的键相互作用都被认为是二面角。这使得不可能区分涉及相同数量原子的几种类型的相互作用;例如，在使用OPLS-AA力场的蛋白质拓扑结构中，二面角采用Ryckaert-Bellemans形式建模，而不适当的二面角采用经典的二面角形式建模，这两种相互作用类型都将由二面角来选择。

当请求基于残基的计算时，如果来自一种残基的至少一个原子和来自第二种残基的一个原子涉及该类型的相互作用，则存在基于某种类型的对对残基的相互作用。

atombased 和 residuebased指定如何处理和保存成对的力。这两个关键字是独立的:一个可以设置为“no”，表示没有兴趣存储、处理和保存该类型的数据;默认为“no”。下表显示了可能的值及其影响:

选项

计算

输出数据

输出文件

no

成对的力不存储、处理或保存。

none

.pfa or .pfr

pairwise_forces_vector

成对力被保存为向量。

vector

.pfa or .pfr

pairwise_forces_scalar

成对力保存为有符号标量。

signed scalar

.pfa or .pfr

punctual_stress

对于每个原子/残基，准时应力被保存。只对onepair=sum有效。

scalar

.psa or .psr

virial_stress

对于每个原子，维里应力被保存了下来。仅适用于基于原子的和onepair=sum。

tensor

.vsa

virial_stress_von_mises

对于每个原子，冯米塞斯维里应力被保存。仅适用于基于原子的和onepair=summed。

scalar

.vma

compat_bin (deprecated option)

成对力作为带符号的标量保存到二进制中格式与原来的PF实现兼容。

只对onepair=summed。

signed scalar

.pfa or .pfr

compat_ascii (deprecated option)

成对的力被保存为签名标量，保存为与原始PF实现兼容的文本格式。只对onepair=summed。

signed scalar

.pfa or .pfr

Vector2scalar -控制如何在矢量对力和它们的标量表示之间进行转换。它可以有以下值:投影-标量是由两个原子的位置或两个残基的COM的位置决定的力向量在矢量上的投影的长度。

范数-标量是矢量的范数，符号是根据矢量力和矢量之间的夹角的余弦值确定的，由两个原子的位置或两个残基的COM的位置决定:

Ø如果它是正的，矢量力的方向是相对于原子/残基的-90到90度方向相同=有吸引力=负的

Ø如果它是负的，矢量力的方向是相对于原子/残基的90到270度方向相反=排斥=正的

Ø如果它是零，矢量力的方向与原子/残基正好成90度或270度——这个力既不吸引也不排斥，所以它被设为零

默认值是'norm'。

Residuesrenumber—控制剩余重新编号。可能的值是:自动保持剩余编号，如在结构文件中，如果没有碰撞和做重新编号，如果碰撞存在是-做剩余重新编号独立的剩余编号碰撞的存在否-不做剩余重新编号。如果检测到碰撞，将写入警告，但不会进行重新编号。使用这个选项是危险的，因为残基不能唯一识别!如果存在残基序号碰撞，具有相同残基序号的残基序号的成对力将被加在一起——这可能是不希望的。

只有当residuebased为“no”之外的其他选项时，才会考虑此选项。默认值为“auto”。

time_averages_period—是平均标量数据之间的输入轨迹帧数。这只能与基于原子和/或基于残留的设置为pairwise_forces_scalar、compat_bin或compat_ascii一起使用。默认值是1 -不做平均。对于这个数量的输入帧，输出将包含一个帧;如果输入帧不是这个数字的倍数，则在对其余输入帧进行平均后写入最后一帧。值为0表示应该对输入轨迹中出现的所有帧进行平均;然后输出将包含单个帧。

no_end_zeros – 控制每个原子/残基数据末尾是否存在零。可能的值是:

yes-抑制每个原子/残基数据末尾的所有零

no-不要在每个原子/残基数据的末尾抑制零

默认值为no。当计算由感兴趣的分子和大量水组成的系统的每个原子/残基数据时，此选项非常有用;在这种情况下，对于感兴趣的分子，每个原子/应力只有非零值，后面跟着许多零，这增加了存储空间，但没有带来任何有用的信息。

energy_grp_exclusion -控制能源组排除的使用，以加快md的运行;

yes-不属于FDA组的原子之间不需要的力将不计算

no-所有的力都会被计算出来

默认值是“yes”。该选项仅适用于排除不能正确调试的情况。

normalize_punctual_stress_per_residue –将每个残基的准时应力除以残基中的原子数。可能的值是:

yes-每个残余点应力将被归一化

no-每个残余点应力将不归一化

默认值为“no”。

binary_result_file – 以二进制格式存储FDA结果文件。可能的值是:

yes-二进制格式的FDA结果文件

no-基于文本格式的FDA结果文件

默认值为“no”。

threshold – 是一个实数，用于忽略所有小于阈值的成对力。单位为kJ/mol/nm，默认值为1e-7。

输入文件示例

一个示例。pfi文件，仅用于文本兼容格式的基于残余的输出，考虑到所有交互类型和所有输入帧的平均时间:

; onepair could be detailed or summed

onepair = summed

; group1 and group2 as defined in the -pfn file

; if not defined, defaults to 'Protein'

group1 = Protein

group2 = Ligand

; no interest in atom-based information

atombased = no

; output residue-based information in compatibility format

residuebased = compat_ascii

; interactions type could be one of more of:

; bond angle dihedral polar coulomb lj nb14 bonded nonbonded all

type = all

; scalar summation over all steps time_averages_period = 0

另一个示例。pfi文件，仅用于基于原子的输出，包含所有键合相互作用的每个原子的总和:

; onepair could be detailed or summed

onepair = summed

; group1 and group2 as defined in the -pfn file

; if not defined, defaults to 'Protein'

group1 = Protein

group2 = Ligand

; sum of scalar pairwise forces per atom

atombased = per_atom_sum

; no interest in residue-based information

residuebased = no

; interactions type could be one of more of:

; bond angle dihedral polar coulomb lj nb14 bonded nonbonded all

type = bonded

另一个例子。pfi文件，用于基于原子的输出包含每个原子最大标量成对相互作用和基于残基的输出包含标量成对力:

; onepair could be detailed or summed

onepair = summed

; group1 and group2 as defined in the -pfn file

; if not defined, defaults to 'Protein'

group1 = Protein

group2 = Ligand

; maximum scalar pairwise force per atom

atombased = per_atom_max

; scalar pairwise forces between residues

residuebased = scalar

; interactions type could be one of more of:

; bond angle dihedral polar coulomb lj nb14 bonded nonbonded all

type = all

最后是一个示例。pfi文件，仅用于文本兼容格式的基于残余的输出，考虑到所有交互类型并执行残余重新编号，对所有输入帧进行时间平均:

; onepair could be detailed or summed

onepair = summed

; group1 and group2 as defined in the -pfn file

; if not defined, defaults to 'Protein'

group1 = Protein

group2 = Ligand

; no interest in atom-based information

atombased = no

; output residue-based information in compatibility format

residuebased = compat_ascii

; ask for residues to be renumbered - starting from 1 and continuous

residuesrenumber = yes

; interactions type could be one of more of:

; bond angle dihedral polar coulomb lj nb14 bonded nonbonded all

type = all

; scalar summation over all steps

time_averages_period = 0

输出文件格式

文件格式独立于它们是否包含基于原子的数据或基于残余的数据。Pfa或.pfr文件)。兼容性格式(基于原子的/基于残留的设置为compat_bin或compat_ascii)应该尽可能接近从原始PF实现中获得的格式。PF2和PF交互类型之间使用以下映射:

PF2

PF

bond

iBond

angle

iAngle

dihedral

iDihedral

polar

iPolar

lj

iLJ

coulomb

iCoul

nb14

iDihedral

然而，兼容性格式包含了更多的数据:最初的PF实现近似地认为，对于一个角度或二面角，只有极端原子相互作用;在这些情况下，PF2正确地处理所有原子之间的相互作用，从而导致输出文件中出现更多的相互作用。

这些格式应该仅用于方便将数据读取到当前的R FDA工具中。包含成对力(基于原子/残基的标量或矢量)的输出格式由一组帧组成，每一帧包含一组变量的相互作用。一个框架以包含以下内容的一行开始:

frame <frame_number>

frame <frame_number>是一个从0开始的整数，由单词“frame”之间的空格分隔。这一行后面跟着零行或多行格式:

i j scalar_force interaction_type

对于基于原子/残基的设置为标量或

i j force_x force_y force_z interaction_type

对于基于原子/残基的集合为向量。I和j是整数，表示原子索引(0为基数);Scalar_force是一个浮点数，表示成对力矢量的大小;force_x, force_y和force_z是浮点数，表示成对力矢量的X, Y和Z分量;interaction_type是一个表示交互类型的整数——对于详细的交互，这是一个在include/pf_interactions.h文件中定义的常量，可能会根据PF2版本而改变(现有的常量通常会被保留，新的常量会被添加更高的值);对于求和的相互作用，这是原子i和j之间存在的所有相互作用对应的常数的和。所有浮点数都使用C printf格式说明符“%e”输出。一行中的所有值都用空格分隔。

包含每个原子/每个残基数据的输出格式由每帧一行组成，每行的值由空格分隔。所有值都是浮点数，表示和/平均/最小/最大标量成对交互;它们是使用C printf格式说明符“%e”编写的。所述线上的第一个值对应于所述第一原子/残基，所述线上的第二个值对应于所述第二原子/残基，依此类推;在每一行上，每个原子的值和原子的数目一样多，每个残基的值和残基的数目一样多。如果原子/残基没有任何成对相互作用(例如:

因为它位于距离其他原子的截止点更远的地方，因为它没有指定的类型与type之间的任何交互作用，或者因为它既没有包含在group1中，也没有包含在group2中，所以对应的值是0.0。

对于每残基序号数据，残基序号编号取决于是否进行了残基序号重新编号。如果没有进行残基序号重新编号，则残基序号从0开始，一直到在输入结构文件中找到最大的数字为止。例如，如果输入结构包含编号为5、6、8和9的4个残基序号，则输出将包含10个值:前5个值和第8个值将始终为零。如果进行了残基序号重新编号，则残基序号从0开始，一直到残基序号总数减去1。

二进制文件格式

FDA结果文件也可以以二进制格式写入，以减少和加速磁盘使用。通过将选项' binary_result_file '设置为' yes '来使用二进制格式。所有FDA结果文件(pfa, pfr, psa, psr, vsa, vma)都可以通过使用转换模块从基于文本的格式转换为二进制格式，反之亦然:

gmx_fda fda_convert -i <input-file> -o <output-file>

如果输入文件是基于文本的，输出文件将是二进制的，反之亦然。

分析模块

所有分析模块都集成在GROMACS中，可以通过以下方式执行:

gmx_fda <module> [options]

每个分析模块的文档可以通过以下方式打印:

gmx_fda help <module>

fda_graph模块

fda_graph将FDA强制网络转换为PDB或DIMACS图。如果使用可选文件-ipf-diff，则将取两两作用力的差异。

PDB图形允许使用您选择的程序轻松可视化。

只考虑大于-t的力。默认阈值为零。

网络包含的节点数必须至少与最小值(默认为2)相同。如果使用-big选项，则只会打印节点数最多的网络。将确定每个网络，并为每个网络分配段名称，因此通过段id为其着色将有助于分析(32种不同的颜色)。Bfactor列将用于表示力的值，并作为力大小的函数帮助着色。CONNECT头将用于在节点之间创建绑定。

fda_shortest_path模块

fda_shortest_path计算FDA强制网络的源节点和目的节点之间的最短k条路径。如果使用可选文件-ipf-diff，则将取两两作用力的差异。图形将以pdb格式打印，这允许您选择的程序轻松可视化。

将确定每个路径，并为每个路径分配段名称，因此通过段id为它们着色将有助于分析(32种不同的颜色)。Bfactor列将用于表示力的值，并作为力大小的函数帮助着色。CONNECT头将用于在节点之间创建绑定。

fda_get_stress模块

fda_get_stress通过成对力计算准时应力。如果使用可选文件-ipf-diff，则将取两两作用力的差异。

fda_view_stress模块

fda_view_stress以xpm或pdb-file表示FDA的准时或von Mises病毒应力。xpm文件的x轴表示粒子数，y轴表示帧数。对于xpm文件，可以使用-nbColors选项设置不同颜色的数量。对于pdb文件，Bfactor列将用于表示应力值，并作为应力大小的函数帮助着色。可以使用您选择的程序对pdb文件进行可视化。