计算相互作用能时正负波动如何理解

zhangfaxue

体系介绍：体系中有ABC分子，想要计算AB分子，AC分子以及AA分子的相互作用能（组分内、组分间）
计算相互作用能的步骤如下：
完整mdp文件：
energygrps     = A B
pbc                = xyz
cutoff-scheme = Verlet
coulombtype   = cut-off
rcoulomb      = 3.0
vdwtype       = cut-off
rvdw          = 3.0
DispCorr      = EnerPre
运行步骤：
gmx trjconv -f traj_comp.xtc -s npt.tpr -n oil.ndx -o AB.xtc
gmx trjconv -f N3030.gro -s npt.tpr -n oil.ndx -o AB.gro
gmx grompp -f rerun.mdp -c AB.gro -p AB.top -o AB.tpr
N3030]$ gmx mdrun -s AB.tpr -rerun AB.xtc -e AB.edr
现在没有报错，使用命令分析能量：gmx energy -f AB.edr -o AB.xvg
选择：
库伦相互作用：Coul-SR:A-B
范德华相互作用：LJ-SR:A-B
其中我的A为带COOH的分子，B为阳离子
他们的库伦相互作用我理解应该是负的表现为吸引力，但是统计整个模拟时间内结果显示正负波动，统计所有时间平均值为正:

请问大家是我理解有问题还是计算处理有问题?

zhangfaxue

sobereva 发表于 2024-5-13 22:25
verlet情况的bug所致（在我来看是bug），之前论坛里我就回答过类似问题

谢谢老师解答

sobereva

zhangfaxue 发表于 2024-5-13 21:29
老师还有一个问题需要请教，我的A为带COOH的分子，B为阳离子，他们的库伦相互作用我理解应该是负的表现为 ...

verlet情况的bug所致（在我来看是bug），之前论坛里我就回答过类似问题

sobereva

zhangfaxue 发表于 2024-5-13 21:12
老师我用cutoff-scheme = group，提示Fatal error:rcoulomb must not be greater than rlist (twin-range ...

把rlist设成和rcoulomb/rvdw相同的值

zhangfaxue

sobereva 发表于 2024-5-10 23:43
用cutoff-scheme = group

老师还有一个问题需要请教，我的A为带COOH的分子，B为阳离子，他们的库伦相互作用我理解应该是负的表现为吸引力，但是为什么我算出来是正，是我renrun.mdp文件设置有问题还是和什么相关？

zhangfaxue

sobereva 发表于 2024-5-10 23:43
用cutoff-scheme = group

老师我用cutoff-scheme = group，提示Fatal error:rcoulomb must not be greater than rlist (twin-range schemes are not supported)
然后逐渐减小rcoulomb和rvdw直到rcoulomb=0.8，rvdw=0.8(都小于1）提示
NOTE 1 【file rerun2.mdp]:The group cutoff scheme is deprecated since GRoMAcs 5.0 and will beremoved in a future release when all interaction forms are supported forthe verlet scheme. The verlet scheme already scales better, and it iscompatible with GPUs and other accelerators
NOTE 2 file rerun2.mdp]:
For a correct single-point energy evaluation with nsteps = 0, usecontinuation =yes to avoid constraining the input coordinates.
ERROR 1 [file rerun2.mdp]:With coulombtype = Cut-off, rcoulomb should be >= rlist unless you use apotential modifier
ERROR 2lfile rerun2.mdp]:With vdwtype = cut-off, rvdw must be >= rlist unless you use a potentialmodifier
NOTE 3 「file rerun2.mdp]:You are using a cut-off for Vdw interactions with NVE, for good energyconservation usevdwtype=shift(possibly with Dispcorr)
NOTE 4 [file rerun2.mdp]:You are using a cut-off for electrostatics with NVE, for good energyconservation use coulombtype=PME-Switch or Reaction-Field-zero
求助老师使用cutoff-scheme = group应该怎么确定rvdw和rcoulomb

sobereva

用cutoff-scheme = group