GROMACS如何使用其他水模型
进行分子动力学模拟，水分子十分重要，除非选择使用连续介质模型（implictit water model）。水分子模型较多，选择这些模型要结合使用的力场，并参考别人已经的数据。以下简单介绍几种常见的水分子模型，希望对了解它们有点帮助。
按照一般化学常识，水分子由三个原子构成，主要的参数应该有各个原子的质量，电量，氢氧键的长度以及H-O-H的键角。没有错，最简单的水分子模型就是这些参数都固定的刚性水分子模型。如SPC模型和TIP3P模型。这两种模型中，原子质量和电量都在同一个质点上。唯一不同的是TIP3P的H-O-H键角比理论值109.47小，为104.52度。这两种水模型只有氧原子具有范德华作用系数，氢原子的范德华系数为0。
以上两种模型有对应的改进模型，SPC的改进模型为SPC/E，起主要改进其实就是使溶液系统的总能量乘以5.22 kJ/mol。这样可以使SPC溶液属性更加接近实验值。TIP3P在CHARMM力场中的改进是给氢原子一定的范德华系数，这样做的结果的计算根据复杂。
由于真实情况下水分子的电量分布并不是完全在原子上的，如氧原子的一部分负电量就在H-O-H的对角线上，还有两个电子对处在H-O化学键的延长线上。为了得到更加真实的水分子模型，四个粒子以上的模型就被应用到分子动力学模拟中。其中最著名的有TIP4P模型。该模型在三个原子中间，H-O-H化学键的对角线上多了一个不含质量，只带电量的点。很多蛋白质模拟计算中，TIP4P和OPLS力场结合使用都得到很好的效果。
以上提到，水分子的氧原子在H-O化学键延长线上有两个电子对，于是有的人就在这两处添加了两个只带电量的粒子。2000年报道的TIP5P模型，计算结果也很好。还有一些牛人，结合TIP4P和TIP5P，要研制TIP6P，很好很强大。不得不说，并不是模型的所含粒子越多越好。粒子越多，就算付出越大，因为要计算的相互作用更多。以上所提到的各个模型的参数，可以看看这里： http://en.wikipe dia.org/wiki/Water_model
在我们MD中我们有时候很少进行水模型的选择，很多时候我们对于不要求精确的体系使用spc水模型，对于相对要求精确的模型一般使用tip3p模型，对于离子或者小分子的研究有时候使用tip4p的模型，这些在gromacs中都是含有的，但是水模型一直在发展，不同的体系可能对于水模型有较大的变化，具体可以查看水模型的综述页面。例如核酸中用的较多的除了tip4p-eW水模型以外还有OPC水模型，但是该水模型并未在gromacs中自带，所以需要自己构建，以下以OPC水模型为例介绍方法:
简化的经典水模型是实际原子模拟中不可缺少的组成部分。然而，尽管经过几十年的深入研究，这些模型还远未完善。我们开发了一种新的方法来构建广泛使用的点电荷水模型，这与主流水模型参数化技术完全不同。与传统方法相比，除了对称性之外，我们不对模型施加任何几何约束。相反，我们优化点电荷的分布以最好地描述水分子的“静电”。使用这种新方法开发了4点OPC和3点OPC3刚性水模型，与常用的刚性模型相比，该模型显示出更为精确地重现大部分的性质。
OPC水模型的topol文件完整如下:
[ atomtypes ]OW OW 0.0000 0.0000 A 3.16655e-01 8.903586e-01HW HW 0.0000 0.0000 A 0.00000e+00 0.00000e+00MW MW 0.0000 0.0000 A 0.00000e+00 0.00000e+00[ moleculetype ]; molname nrexclSOL 2[ atoms ]; id at type res nr res name at name cg nr charge mass1 OW 1 SOL OW 1 0 16.000002 HW 1 SOL HW1 1 0.67914 1.008003 HW 1 SOL HW2 1 0.67914 1.008004 MW 1 SOL MW 1 -1.35828 0.00000#ifndef FLEXIBLE[ settles ]; i funct doh dhh1 1 0.08724 0.13712#else[ bonds ]; i j funct length force.c.1 2 1 0.08724 502416.0 0.08724 502416.01 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 b4 1 2 3 1 0.147722363 0.147722363[ exclusions ]1 2 3 42 1 3 43 1 2 44 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)
为了防止污染原来的水模型，应当把原拓扑中的原子类型增加一个后缀，如把HW改成HW_opc，OW改成OW_opc。更改过后如下：
[ atomtypes ]OW_opc OW_opc 0.0000 0.0000 A 3.16655e-01 8.903586e-01HW_opc HW_opc 0.0000 0.0000 A 0.00000e+00 0.00000e+00MW MW 0.0000 0.0000 A 0.00000e+00 0.00000e+00[ moleculetype ]; molname nrexclSOL 2[ atoms ]; id at type res nr res name at name cg nr charge mass1 OW_opc 1 SOL OW 1 0 16.000002 HW_opc 1 SOL HW 1 0.67914 1.008003 HW_opc 1 SOL HW 1 0.67914 1.008004 MW 1 SOL MW 1 -1.35828 0.00000#ifndef FLEXIBLE[ settles ]; i funct doh dhh1 1 0.08724 0.13712#else[ bonds ]; i j funct length force.c.1 2 1 0.08724 502416.0 0.08724 502416.01 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 b4 1 2 3 1 0.147722363 0.147722363[ exclusions ]1 2 3 42 1 3 43 1 2 44 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)
之后打开安装目录下需要使用的力场里面的ffnobonded.itp，在[ atomtypes ]位置下写入：
OW_opc OW_opc 0.0000 0.0000 A 3.16655e-01 8.903586e-01HW_opc HW_opc 0.0000 0.0000 A 0.00000e+00 0.00000e+00MW MW 0.0000 0.0000 A 0.00000e+00 0.00000e+00
注意上述为amber力场的格式，即采用的是LJ6-12势函数，若要用到其他力场中去，注意其势函数是否一致，不一致需要统一换算。然后将上诉水的[ atomtypes ]部分删除，保存为opc.itp的文件。然后在总的top文件里面调用此水模型的路径即可。
最后gromacs中tip4p的gro文件，即可。在模拟时候可能会弹出一个警告说gro和top里面原子名字不匹配，直接使用 -maxwarn 选项忽略即可。值得注意的是opc水模型在md过程中需要进行长程色散校正或者使用Lennard Jones PME,即以下mdp设置二选一（默认的gromacs教程中有1设置）
1. DispCorr = EnerPres (Tested)2. vdwtype = PME
