首页/文章/ 详情

GROMACS如何使用其他水模型

1月前浏览1861

进行分子动力学模拟,水分子十分重要,除非选择使用连续介质模型(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    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     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    mass  1   OW_opc      1       SOL       OW       1       0        16.00000  2   HW_opc      1       SOL       HW       1       0.67914   1.00800  3   HW_opc      1       SOL       HW       1       0.67914   1.00800  4   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





来源:模拟之家
化学电子理论分子动力学GROMACS
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-09-08
最近编辑:1月前
刘十三613
博士 分子动力学、GROMACS
获赞 133粉丝 95文章 82课程 29
点赞
收藏
作者推荐
未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈