gromacs使用额外水模型

在我们MD中我们有时候很少进行水模型的选择，很多时候我们对于不要求精确的体系使用spc水模型，对于相对要求精确的模型一般使用tip3p模型，对于离子或者小分子的研究有时候使用tip4p的模型，这些在gromacs中都是含有的，但是水模型一直在发展，不同的体系可能对于水模型有较大的变化，具体可以查看水模型的综述页面。例如核酸中用的较多的除了tip4p-eW水模型以外还有OPC水模型，但是该水模型并未在gromacs中自带，所以需要自己构建，以下简单介绍方法:

关于OPC水模型

简化的经典水模型是实际原子模拟中不可缺少的组成部分。然而，尽管经过几十年的深入研究，这些模型还远未完善。我们开发了一种新的方法来构建广泛使用的点电荷水模型，这与主流水模型参数化技术完全不同。与传统方法相比，除了对称性之外，我们不对模型施加任何几何约束。相反，我们优化点电荷的分布以最好地描述水分子的“静电”。我们使用这种新方法开发了4点OPC和3点OPC3刚性水模型，与常用的刚性模型相比，该模型显示出更为精确地重现大部分的性质。

以OPC水模型为例

首先进行下载GROMACS的OPC，topol文件，若没有带拓扑文件，可以使用例如ACPYPE进行转化。
OPC水模型的topol文件完整如下:

[ atomtypes ]
OW            OW      0.0000  0.0000  A   3.16655e-01  8.903586e-01
HW            HW      0.0000  0.0000  A   0.00000e+00  0.00000e+00
MW            MW      0.0000  0.0000  A   0.00000e+00  0.00000e+00

[ moleculetype ]
; molname       nrexcl
SOL             2

[ atoms ]
; id  at type     res nr  res name  at name  cg nr  charge    mass
  1   OW          1       SOL       OW       1       0        16.00000
  2   HW          1       SOL       HW1      1       0.67914   1.00800
  3   HW          1       SOL       HW2      1       0.67914   1.00800
  4   MW          1       SOL       MW       1      -1.35828   0.00000

#ifndef FLEXIBLE

[ settles ]
; i     funct   doh     dhh
1       1       0.08724 0.13712

#else

[ bonds ]
; i     j       funct   length  force.c.
1       2       1       0.08724 502416.0 0.08724        502416.0
1       3       1       0.08724 502416.0 0.08724        502416.0

[ angles ]
; i     j       k       funct   angle   force.c.
2       1       3       1       103.6   628.02  103.6  628.02

#endif


[ virtual_sites3 ]
; Vsite from                    funct   a               b
4       1       2       3       1       0.147722363     0.147722363


[ exclusions ]
1       2       3       4
2       1       3       4
3       1       2       4
4       1       2       3


; The position of the virtual site is computed as follows:
;
;               O
;             
;               V
;         
;       H               H
;
; Ewald OPC:
; const = distance (OV) / [ cos (angle(VOH))    * distance (OH) ]
;         0.01594 nm     / [ cos (51.8 deg)     * 0.0872433 nm    ]
;       then a = b = 0.5 * const = 0.147722363
;
; Vsite pos x4 = x1 + a*(x2-x1) + b*(x3-x1)

我们可以在topol或者ffnobonded.itp中恰当的位置写入

;opc
HW_opc       1       1.008   0.0000  A   0.00000e+00  0.00000e+00
OW_opc       8      16.00    0.0000  A   3.16655e-01  8.903586e-01

为了防止污染我把原拓扑中的原子类型增加了后缀，以免报错。然后将原top中的[ atomtypes ]部分删除，重命名为opc.itp.

由于我们使用的为4点电荷，我是使用的gromacs中tip4p那一套，反正要进行预平衡的，当然也可以自己预平衡一个水，制作一个gro文件。简单介绍一下我用tip4p选择水模型以后，仅需将topol文件中的水拓扑文件替换为如下:

; Include water topology
#include "./amber_na.ff/opc.itp"

值得注意的是opc水模型在md过程中需要进行长程色散校正或者使用Lennard Jones PME,即以下mdp设置二选一（默认的gromacs教程中有1设置）

1. DispCorr = EnerPres (Tested)
2. vdwtype = PME

参考资料:https://bioinformatics.cs.vt.edu/~izadi/