提取单个水分子势能问题

gongrehk

请问Gromacs跑完后如何提取单个水分子的势能呢?比如说采用SPC/E模型模拟白质溶液体系。
gmx energy 提示可以通过以下步骤计算想要组分的势能：
1). 先建.index文件明确感兴趣的组，然后根据原来.tpr文件建立一个新的.tpr文件。
2). 使用gmx mdrun的-rerun选项指定原轨迹再跑一次模拟
3). 使用gmx energy提取感兴趣的能量项
这个方法的话，我照着试了一下
1) gmx convert-tpr -s md.tpr -n water1.ndx -o mdindex.tpr  (比如提取第一个水分子的3个原子, 系统共含33753个原子)
   该命令提示:Select a group:
   然后我选择:Selected 18: 'r_130'
             Will write subset r_130 of original tpx containing 3 atoms......
2) gmx mdrun -s mdindex.tpr -rerun md.trr -deffnm mdindex
      然后就报错了，提示：Number of atoms in trajectory (33753) does not match the run input file (3)

我还用了在.mdp文件中添加energygrps= r_130 的方法， rerun的时候系统报错提示能量可能未平衡好，很迷惑。
请问有前辈了解过这方面的问题吗?该怎么解决呢?

谢谢大家！

k64_cc

gongrehk 发表于 2021-8-23 18:19
reciprocal-space内的算法我没有具体了解， 不过PME算法不也是计算每个原子的静电势能吗？不然怎么算作用 ...

PME不是pairwise的，不支持能量分解。

gongrehk

sobereva 发表于 2021-8-22 01:44
PME的情况下从原理上就根本没法获得组之间的静电相互作用能，必须用cutoff。你把PME的原理搞懂了自然就知 ...

reciprocal-space内的算法我没有具体了解， 不过PME算法不也是计算每个原子的静电势能吗？不然怎么算作用力呢？既然如此，那么每个分子的势能就应该能获取啊？

sobereva

gongrehk 发表于 2021-8-21 17:18
谢谢老师的回复。
请问老师 把coulombtype设为cutoff是指像计算L-J势一样只计算截断半径内的静电相互作 ...

PME的情况下从原理上就根本没法获得组之间的静电相互作用能，必须用cutoff。你把PME的原理搞懂了自然就知道为什么没法实现

gongrehk

sobereva 发表于 2021-8-21 13:28
你得用cutoff算静电作用，即rerun用的mdp里必须把coulombtype设为cutoff

谢谢老师的回复。
请问老师 把coulombtype设为cutoff是指像计算L-J势一样只计算截断半径内的静电相互作用势能吗
？
MD计算的时候是用的PME啊。我试了一下，用Cut-off得到的Coul-SR:r_130-r_130确实是0，不过这两种方法的结果得到的其他Coul项还是有一点差异的。
还是有些不明白为什么用PME的时候，Coul-SR:r_130-r_130这一项结果不是0。

sobereva

gongrehk 发表于 2021-8-20 18:05
老师@sobereva
按照您的方法，mdrun加上-ntmpi 1 可以运行了。
不过过有个问题，对于水分子r_130来说 ...

你得用cutoff算静电作用，即rerun用的mdp里必须把coulombtype设为cutoff

gongrehk

sobereva 发表于 2021-8-18 22:16
当然可以energygrps= r_130。此时edr会出现r_130 - r_130、r_130 - rest、rest - rest三类作用项，rest代 ...

老师@sobereva
按照您的方法，mdrun加上-ntmpi 1 可以运行了。
不过过有个问题，对于水分子r_130来说，Coul-SR:r_130-r_130这一项为什么不是0呢？我用的OPLS-AA力场，1-2，1-3之前的非键作用不应该是0吗?

@ s0 legend "Coul-SR:r_130-r_130"
@ s1 legend "Coul-14:r_130-r_130"
@ s2 legend "Coul-SR:r_130-rest"
@ s3 legend "Coul-14:r_130-rest"
@ s4 legend "Coul-SR:rest-rest"
@ s5 legend "Coul-14:rest-rest"
    0.000000   -4.046188    0.000000  -79.345924    0.000000  -633651.750000  7718.288574
    0.800000   -4.046196    0.000000  -97.905640    0.000000  -633431.812500  7741.199219
    1.600000   -4.046204    0.000000  -110.546257    0.000000  -633271.687500  7757.114746
    2.400000   -4.046173    0.000000  -88.627541    0.000000  -635106.875000  7632.303711

sobereva

gongrehk 发表于 2021-8-18 17:36
energygrps 是用来计算两个group之间的相互作用能的是吗？所以不能只设定energygrps= r_130。
应该再设定 ...

当然可以energygrps= r_130。此时edr会出现r_130 - r_130、r_130 - rest、rest - rest三类作用项，rest代表体系其余部分
mdrun加上-ntmpi 1再试。当前写-s mdindex.tpr是多余的

gongrehk

energygrps 是用来计算两个group之间的相互作用能的是吗？所以不能只设定energygrps= r_130。
应该再设定一个index文件包含所有其他的原子
然后设定energygrps= r_130 !r_130
我这样重新执行了一遍
1) gmx grompp -f mdindex.mdp -p topol.top -c npt.gro -t npt.cpt -n water1_others.ndx -o mdindex
2) gmx mdrun -s mdindex.tpr -rerun md.trr -deffnm mdindex
还是报错
Fatal error:
15 particles communicated to PME rank 14 are more than 2/3 times the cut-off
out of the domain decomposition cell of their charge group in dimension y.
This usually means that your system is not well equilibrated.

gongrehk

谢谢老师@sobereva
我又试了一遍
首先在原来的md.mdp文件中添加一行 ‘energygrps= r_130’，生成新的mdindex.mdp文件 然后执行以下命令
1) gmx grompp -f mdindex.mdp -p topol.top -c npt.gro -t npt.cpt -n water1.ndx -o mdindex
2) gmx mdrun -s mdindex.tpr -rerun md.trr -deffnm mdindex
第一步没有问题， 第二步就报错了
Fatal error:
15 particles communicated to PME rank 14 are more than 2/3 times the cut-off
out of the domain decomposition cell of their charge group in dimension y.
This usually means that your system is not well equilibrated.
提示能量可能未平衡好， 然后我就不知道怎么解决了。

sobereva

tpr文件必须和trr文件里的原子数相同

你就把水分子单独设一个索引组，energygrps指定它，产生tpr文件，然后结合之前的轨迹文件rerun一下就完了。完全不需要convert-tpr

甭管那个提示