能量分析g_energy求助

diaok

我在尝试模拟opls力场下的纤维素
用的是topbuild工具产生的拓扑结构
top文件里面我是直接
#include "oplsaa.ff/forcefield.itp"  
#include "oplsaa.ff/tip3p.itp"
#include "ce8r.itp"
#include "oh.itp"   （我把力场里面的离子都抄到这里了，不过水也有问题，应该不是这里出错  ）

然而分析能量的时候发现除了纤维素自身其他能量项全是0
在md.mdp里面用的
energygrps               = BGC SOL NA OH;  （这里如果不区分应该分析不了吧？）

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Coul-SR:BGC-BGC             -122143        230    1431.08   -1455.24  (kJ/mol)
LJ-SR:BGC-BGC               14208.1         54    395.806    354.901  (kJ/mol)
Coul-14:BGC-BGC             3072.41        1.9    17.3437    13.0047  (kJ/mol)
LJ-14:BGC-BGC               213.496       0.37    11.8793  -0.271114  (kJ/mol)
Coul-SR:BGC-SOL                   0          0          0          0  (kJ/mol)

然后计算的log文件也是
        <====  A V E R A G E S  ====>
        <==  ###############  ======>

        Statistics over 10000001 steps using 100001 frames

   Energies (kJ/mol)
           Bond          Angle Ryckaert-Bell.          LJ-14     Coulomb-14
    3.02468e+02    6.88317e+02    5.77314e+02    2.15790e+02    3.07427e+03
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
    1.42261e+04   -1.22430e+05    3.50924e+02   -1.02995e+05    1.08323e+04
   Total Energy    Temperature Pressure (bar)
   -9.21624e+04    2.67223e+02   -1.40703e+02

          Box-X          Box-Y          Box-Z
    3.34810e+00    3.34810e+00    4.15200e+00

   Total Virial (kJ/mol)
    3.60646e+03    1.09681e+00   -6.75310e-01
    1.07644e+00    3.60586e+03   -1.06565e+00
   -6.62611e-01   -1.02418e+00    4.21244e+03

   Pressure (bar)
    2.56993e+00   -1.83146e-01    3.20050e-01
   -1.68682e-01    2.56961e+00    4.37014e-01
    3.10924e-01    4.07437e-01   -4.27250e+02

  Epot (kJ/mol)        Coul-SR          LJ-SR        Coul-14          LJ-14
        BGC-BGC   -1.22430e+05    1.42261e+04    3.07427e+03    2.15790e+02
        BGC-SOL    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00
         BGC-NA    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00
         BGC-OH    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00
        SOL-SOL    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00
         SOL-NA    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00
         SOL-OH    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00
          NA-NA    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00
          NA-OH    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00
          OH-OH    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00

怎么都想不通是哪里出了问题
求大神指导

diaok

居然找到原因了。。

使用gpu计算就会这样
不用cpu正常  （但是这样gpu的计算可信吗？）

gromacs怎么这么多小bug啊
之前有一次mpi版的gyrate数据不正常，折腾我几个星期，安装了各种版本

sobereva

不能算是bug，只是GPU运行时一些功能会有限制。GPU跑的动力学轨迹是没问题的，但不是所有特征、功能都available。

diaok

sobereva 发表于 2015-11-7 23:08
不能算是bug，只是GPU运行时一些功能会有限制。GPU跑的动力学轨迹是没问题的，但不是所有特征、功能都avail ...

了解了

还好rerun之后能分析

另外我想把无限长的纤维素固定在体系中央，但是使用了
comm-grps                = BGC Water_and_ions
之后纤维素还是到处飘，似乎是受到了npt的影响

lammps中npt可以把盒子正中心放在0处，npt时坐标都是成比例变化，质心不会移动
但是从gromacs轨迹看只定义了盒子大小，所以盒子起点都是0，大部分坐标都大于0
是需要靠COM pulling或者freeze的方法吗？

通过轨迹后处理来实现似乎也很麻烦，可能要逐帧操作

或者我的初始状态就不应该固定在0而应该在正中间？

求指导

sobereva

diaok 发表于 2015-11-8 20:01
了解了

还好rerun之后能分析

给纤维素原子加个坐标的限制势就完了

diaok

sobereva 发表于 2015-11-9 00:41
给纤维素原子加个坐标的限制势就完了

经过验证发现放在盒子中间不会移动    就可以不用限制了

应该还是0位置边界的影响

另外模拟过程中1.5mol/L的NaOH和LiOH溶液都结晶了
这个只能靠降低浓度了吗？
Li、Na和OH都是opls力场里有的，换成其他力场的参数怎么样？

sobereva

diaok 发表于 2015-11-10 15:51
经过验证发现放在盒子中间不会移动    就可以不用限制了

应该还是0位置边界的影响

可以尝试换力场。
有些力场离子参数不合理，结晶倾向过强。

diaok

sobereva 发表于 2015-11-10 16:47
可以尝试换力场。
有些力场离子参数不合理，结晶倾向过强。

不过其他力场中很少有同时含Li和Na的

amber有，然而没发现与纤维素相关的文献

需要找到合适的参数自己来加入力场？

sobereva

diaok 发表于 2015-11-10 19:55
不过其他力场中很少有同时含Li和Na的

amber有，然而没发现与纤维素相关的文献

跟amber兼容的glycam是模拟多糖的，模拟纤维素也是有的

KBFF力场（http://kbff.chem.k-state.edu/）从Li到Cs都有，直接提供了gromacs的格式，可以考虑结合gmx自带力场用。

diaok

sobereva 发表于 2015-11-10 21:02
跟amber兼容的glycam是模拟多糖的，模拟纤维素也是有的

KBFF力场（http://kbff.chem.k-state.edu/） ...

非常感谢sobereva老师的耐心指导

不过我在建模方面还是新手，拓扑文件全是靠网上的工具生成的
现在只会连接已有的重复单元   内部的键角、二面角参数太过复杂还没法修改
很难更换力场

把KBFF力场的离子引入opls中似乎更可行
由于KBFF离子的combination rules存在一个scale factor
所以需要把它和力场中的交叉项全部罗列出来
然而我翻遍了gromacs的力场也没看到ffnonbonded中的两两组合项，这种方法应该可行吧？
或者我需要用table来定义？  

diaok

已经换成kbff的离子了
1.5m的NaOH依然还是结晶。。。

kbff网站上的kbffpv2.ff.tar显示文件不存在
我就按照文献数据自己设定的参数
然而碰到了一个问题
The following combination rules used:
σij=(σiiXσij)^1/2
εij= s(εiiXεij)^1/2
The value of s was set to unity for all interactions except for cation to water oxygen
除了氧其他都用这种组合规则
然而我乘了下发现阳离子与氧用这个组合规则得出的结果和表上面单独列出来的一致

这里的unity是指其他组合规则全用1吗？
如果这样似乎没必要引入这个S因子         （ 经过尝试，好像unity就是1。。。


不过我的体系里面阴离子是OH-，Na-Na、Na-O的参数我都是直接写出来的，没有用到组合率
应该是OH-这边出了问题，我用的是文献中的SPC/E改OH，氧带电-1.4238