【教程】如何在gmx中使用列表势

Lacrimosa

gromacs自2020版本起砍掉了列表势功能，因此本教程只适用于gromacs2020以前的版本。列表势可以有效地扩充gromacs程序支持的力场形式，本文将结合gromacs手册简述一下列表势的用法，本人水平有限，欢迎诸位指正。
gromacs中的非键作用项可以用以下形式表示：

（1）

（1）式右端三项分别对应着:静电势，LJ势中的色散项和互斥项。
（1）中的静电势函数在gromacs中表示如下：





其中的相对介电常数epsilon_r可在mdp文件的epsilon-r参数中进行设置，默认为1。

（1）中的LJ势默认为6-12型，其在gromacs中表示如下：

  

与（1）式对应可知，



当使用列表势时，需要以列表的形式提供 r, f(r), -f'(r), g(r), -g'(r), h(r), -h'(r)七项，在gromacs程序的share/gromacs/top路径下有列表势示例： table6-8/9/10/11/12.xvg，其内容大致如下：
table6-12.xvg

# Table AB1: ndisp=6 nrep=12
# r                           f(r)                        -f'(r)                        g(r)                      -g'(r)                        h(r)                       -h'(r)
0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
2.0000000000e-03   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
4.0000000000e-03   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
6.0000000000e-03   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
8.0000000000e-03   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
1.0000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
1.2000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
1.4000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
1.6000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
1.8000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
2.0000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
2.2000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
2.4000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
2.6000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
2.8000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
3.0000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
3.2000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
3.4000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
3.6000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
3.8000000000e-02   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
4.0000000000e-02   2.5000000000e+01 6.2500000000e+02   -2.4414062500e+08 -3.6621093750e+10   5.9604644775e+16 1.7881393433e+19
4.2000000000e-02   2.3809523810e+01 5.6689342404e+02   -1.8218149332e+08 -2.6025927617e+10   3.3190096508e+16 9.4828847166e+18
4.4000000000e-02   2.2727272727e+01 5.1652892562e+02   -1.3781101808e+08 -1.8792411556e+10   1.8991876704e+16 5.1796027375e+18
4.6000000000e-02   2.1739130435e+01 4.7258979206e+02   -1.0554872947e+08 -1.3767225583e+10   1.1140534293e+16 2.9062263373e+18
4.8000000000e-02   2.0833333333e+01 4.3402777778e+02   -8.1762201338e+07 -1.0220275167e+10   6.6850575676e+15 1.6712643919e+18
5.0000000000e-02   2.0000000000e+01 4.0000000000e+02   -6.4000000000e+07 -7.6800000000e+09   4.0960000000e+15 9.8304000000e+17

...

以上表格可以通过excel或编程等方式轻松实现。需要注意以下几点：1. 表格中r<4.0e-02nm时，其余六项均为0。通常这一设置不会出问题，由于该距离很小，两个原子还未达到这个距离就会因为互斥项分开（如下图所示）；但若初始结构中存在两个原子间距小于4.0e-02 nm时，二者间的能量和力均为0，导致这种异常结构可能不会被察觉到。




2. r应从0开始设置，通常每隔0.002nm写入一行（对于双精度的计算通常以0.0005nm为间隔）。最后一行的r至少要大于rc+1，rc为计算静电作用以及范德华作用的截断半径，其参数可以在mdp文件中进行设置，示例如下：
;===mdp#1===
rlist = 1.0
coulombtype = User ;使用列表势计算静电作用，若不想使用列表势计算静电作用可将此处改为PME
rcoulomb = 1.0 ；短程静电作用的截断距离 (nm)
vdwtype = User ;使用列表势计算范德华作用，若不想使用列表势计算范德华作用可将此处改为cut-off
rvdw = 1.0 ;范德华作用的截断距离 (nm)
;=========如按mdp#1设置，rc=1 nm，最终table.xvg中最后一行的r应大于2nm。

mdp#1的设置适用于对体系内全部粒子均使用同一table进行计算的情况（例如对所有原子均使用table6-10.xvg）。
如需对体系内某些特定组分间使用特定的table则应准备多个table并按table_i_j.xvg的方式进行命名。
例如：系统内包含A,B,C,D四种原子，其中A与所有原子的LJ函数形式为6-10，B，C，D原子间的LJ函数形式为6-12。此时应准备table_A_A.xvg，table_A_B.xvg，table_A_C.xvg，table_A_D.xvg四个文件，其内容与table6-10.xvg相同，还应额外准备一个table.xvg文件以供计算B，C，D间的相互作用使用，其内容应与table6-12.xvg内容相同。在mdp文件中，还要额外设置以下参数：
;===extra===
cutoff-scheme            = group
energygrps      = A B C D ; 定义四个需要计算能量的组，该组中任意两个组间（除energygrp-table中定义的）的相互作用将使用table.xvg进行计算（该计算不支持GPU）
energygrp-table = A A A B A C A D ;对需要用额外table计算组进行定义，该组中的相互作用将使用table_i_j.xvg进行计算

;=======

3. LJ函数中6-8/9/10/11/12型式的table.xvg都已经在gromacs路径下提供了，如有额外需求，可对table中的 f(r), -f'(r), g(r), -g'(r), h(r), -h'(r)项进行修改。
例：若要使用6-13型的LJ势，则要对h(r)以及-h'(r)项进行修改，其函数应如下式，之后按照下式计算table中的数值即可。



======================================================================================================================
以上内容主要涉及到式（1）中的f(r),g(r),h(r)项，其余的参数q,C6,C12则都在top文件中设置。

top文件的设置：top文件起始于[ defaults ]的设置（见以下示例），

;===top_file===
[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               1               yes             0.5     0.8333
;============


其参数的含义如下：

1. nbfunc：设置范德华作用的函数形式，1为LJ势，2为Buckingham势
2. comb-rule：设置LJ势函数中参数的形式以及结合规则；设为1时，参数为C6，C12，结合规则为几何平均；设为2时，参数为sigma，epsilon，结合规则为L-B规则；设为3时，参数为sigma，epsilon，结合规则为几何平均。

以上内容可总结为下表：
#nbfunc + comb-rule    function    combination rules
---------------------------------------------------------------------------------------------------------------------------------
1+1  ：

1+2  ：

1+3：
---------------------------------------------------------------------------------------------------------------------------------

当使用table时，只有组合1+1对应着以C6/C12为参数的LJ函数，因此nbfunc和comb-rule都应设为1。
若你当前的力场参数为sigma/epsilon则可通过下式转换：



3. gen-pairs：决定是否自动产生缺少的1-4范德华参数，与top中的[ pairs ]项有关
4. fudgeLJ/QQ：计算1-4作用时范德华项/静电项要乘的系数

[ defaults ]项设置好后，下面就是在[ atomtypes ]中定义原子类型。其内容大致如下：
[ atomtypes ]
; name       at.num   mass    charge  ptype  V                  W
  HW_2005     1       1.008    0.0000   A     0.00000e+00  0.00000e+00
  OW_2005     8       16.00    0.0000   A     3.15890e-01  7.74898e-01
  MW             0        0.000    0.0000   D    0.00000e+00  0.00000e+00


其参数含义如下：
1. name：原子类型的名称，可自行定义，需避免重复
2. at.num：原子序数，对模拟没有实质影响
3. mass：原子质量
4. charge：原子电荷
原子电荷以及质量会在此处和[ moleculetype ]项中依次赋值两次。
5. ptype：粒子类型，A为原子，D为虚原子
6. V/W为LJ势中用到的两个参数，V对应着sigma或C6；W对应着epsilon或C12；具体对应哪个取决于[ defaults ]中的设置。根据以上[ defaults ]设定，这里应为C6，C12。

如若力场中直接给出了原子i与原子j之间的C6(ij)和C12(ij)，则应使用[ nonbond_params ]项进行参数的指定。其示例如下：
[ nonbond_params ]
;i              j               func  V        W
OW_2005 HW_2005   1       0.1    0.1

其参数含义如下：
1. i/j表示定义原子类型i与j之间的范德华作用参数
2. func 定义参数的类型，1表示LJ势，sigma/epsilon 或 C6/C12具体取决于[ defaults ]中的设置。根据以上设定，这里应为1。
3. V/W 原子类型i与原子类型j之间的sigma/epsilon 或 C6/C12参数，具体取决于[ defaults ]中的设置。根据以上设定，这里应为C6(ij)，C12(ij)。
当[ atomtypes ]中给出了ij原子的C6(ii)，C12(ii)，C6(jj)，C12(jj)同时设置了[ nonbond_params ]中的C6(ij)和C12(ij)时，[ nonbond_params ]中的参数优先级更高。








后续更新内容：

1.WCA potential

牧生

可是我没有权限看这一篇

lyj714

咦。话说回来，gromacs用户table用的不多也是有原因的，首先新版本不支持table，只能啃旧版本gmx（<=2019)。另外计算只能用cpu，速度拉胯，这和gmx速度快相违背。所以就导致用的人不多，极少数用的也只在测试开发力场参数这块用，大体系用这东西速度完全无法忍受。要是开发者能在这方面改善就好了，可惜鱼和熊掌不可兼得，他们舍弃了，而追求速度去了。

Lacrimosa

lyj714 发表于 2023-3-14 17:40
咦。话说回来，gromacs用户table用的不多也是有原因的，首先新版本不支持table，只能啃旧版本gmx（

确实，GPU不支持这点其实让table的存在显得很尴尬

小瓶不会做计算

请问大佬有生成table文件的脚本吗，求分享

Lacrimosa

小瓶不会做计算 发表于 2025-5-23 11:12
请问大佬有生成table文件的脚本吗，求分享

在excel里写好公式就可以呀