计算化学公社

标题: 我发现gromacs计算能量出现一个bug，居然能计算出某个离子和它自己的coul势不为0 [打印本页]

作者
Author: cavalier 时间: 2020-10-6 20:32
标题: 我发现gromacs计算能量出现一个bug，居然能计算出某个离子和它自己的coul势不为0
本帖最后由 cavalier 于 2020-10-6 20:38 编辑

我在体系里仅加了一个NA离子，并且在能量组中将它称作NA，实际上这么做是为了rerun后能分析能量，总之rerun后发现发现能输出NA-NA之间的Coul相互作用。
绘制后发现这个值是恒定的，但是不为0，想必是一个bug吧。
(, 下载次数 Times of downloads: 27)
注：NA-NA的LJ势是0，其余数据均是正常的。
担心会不会在计算多原子组内的能量时也出现某个原子会和它自己算出非0能量

作者
Author: jimulation 时间: 2020-10-6 22:36
会不会是与镜像之间的相互作用

作者
Author: sobereva 时间: 2020-10-7 00:08
库仑作用随距离衰减很慢。我不知道你怎么算的静电相互作用，应当在贴子里给出mdp。不为零大抵是和镜像的相互作用

作者
Author: cavalier 时间: 2020-10-7 14:24

sobereva 发表于 2020-10-7 00:08
库仑作用随距离衰减很慢。我不知道你怎么算的静电相互作用，应当在贴子里给出mdp。不为零大抵是和镜像的相 ...

感谢jimulation和sob老师的回复，我感觉很有道理。
我用来rerun的mdp如下：
title             =  BPTI in water, 373K
cpp                =  /lib/cpp
integrator       =  md
dt                =  0.001 ; ps !
nsteps             =  150000000
comm-mode             = None
nstcomm          =  0
nstxout          =  0
nstvout          =  0
nstfout          =  0
nstxtcout          =  2500
nstlog             =  2500
nstenergy          =  0
nstlist          =  10
ns_type          =  grid
coulombtype       =  PME
rlist             =  1.4
rcoulomb          =  1.4
rvdw             =  1.4

Tcoupl             =  v-rescale
tc-grps                = SOL CNT CNP NIT NA CB CAP CAN NB
tau_t             = 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
ref_t             = 300 300 300 300 300 300 300 300 300

energygrps       =  SOL CNT CNP NIT NA CB CAP CAN NB
acc_grps          =  CAP CAN
accelerate       =  0.0 0.0 0.7 0.0 0.0 -0.7
freezegrps       = CNT  CNP CAP CAN
freezedim          = Y Y Y Y Y Y Y Y N Y Y N
Pcoupl             = no
gen_vel          =  no

我的盒子是3*3*40nm，很可能是xy方向上比较小导致了这个情况的发生，这么说，我的模拟中NA离子也会受到它的镜像影响吗？那岂不是我的模拟本身就不合理了！
抑或是仅energy命令计算时有这个问题？

作者
Author: k64_cc 时间: 2020-10-7 14:36

cavalier 发表于 2020-10-7 14:24
感谢jimulation和sob老师的回复，我感觉很有道理。
我用来rerun的mdp如下：
title = BP ...

受镜像影响主要体现在长程方面，你可以考虑看一下Ewald Summation。

如果是个常数就合理，不用管。周期性体系里加一个粒子，能量比空盒子高，也挺合理的。

作者
Author: ulosggs 时间: 2020-10-7 20:44
本帖最后由 ulosggs 于 2020-10-7 20:53 编辑

Ewald计算长程作用要求体系是电中性。不然得出来的能量数值是没有意义的。

作者
Author: k64_cc 时间: 2020-10-7 21:15

ulosggs 发表于 2020-10-7 20:44
Ewald计算长程作用要求体系是电中性。不然得出来的能量数值是没有意义的。

实践上倒是会加背景电荷平衡，不过确实这种设置会平添风险。

作者
Author: jimulation 时间: 2020-10-7 21:57

cavalier 发表于 2020-10-7 14:24
感谢jimulation和sob老师的回复，我感觉很有道理。
我用来rerun的mdp如下：
title = BP ...

rerun时coulombtype应设为cut-off，且其值不大于最短边长的一半

欢迎光临计算化学公社 (http://bbs.keinsci.com/)