为什么计算精氨酸和水的rdf总是小于1

frige

本人计算初学者，按照文献想计算精氨酸上酰胺键的O和水中O的rdf图，模拟结果中rdf总是小于1，没有显示氢键的峰

一步步寻找问题四五天了，不知道是哪里出错，请大家帮忙看看，万分感激！

rdf图对比： 左自己模拟结果 右文献结果（10.1021/acs.cgd.1c01468）



L-精氨酸结构式以及文献中静电势--表明最可能形成氢键的位置是酰胺键上的O：



1 模拟流程

创建盒子
gmx insert-molecules -ci L-citrulline.gro -o box.gro -nmol 8 -box 4 4 4
gmx solvate -cp box.gro -cs spc216.gro -o box2.gro -p topol.top
制作top文件
能量最小化
gmx grompp -f em.mdp -c box2.gro -p topol.top -o minim.tpr -maxwarn 2
gmx mdrun -v -s minim.tpr
NVT平衡
gmx grompp -f nvt.mdp -p topol.top -c confout.gro -o nvt-eq.tpr -maxwarn 2
gmx mdrun -v -s nvt-eq.tpr
NPT平衡
gmx grompp –f npt.mdp –p topol.top –c confout.gro –o npt-eq.tpr -maxwarn 2
gmx mdrun –v –s npt-eq.tpr
成品输出
gmx grompp –f npt-run.mdp –p topol.top –c confout.gro –o npt-run.tpr -maxwarn 2
gmx mdrun –v –s npt-run.tpr
分析结果
gmx make_ndx -f confout.gro -o index.ndx
gmx rdf -f traj.trr -s npt-run.tpr -n index.ndx -o rdf.xvg -ref 6 -sel 7


(NVT NPT NPT-run 均步长2fs constraints = hbonds 跑的时间是0.1 0.1 0.2 ns)


2 index索引设置
ref 6  设置为酰胺上的氧

[ CONH-O ]
   3   28   53   78  103  128  153  178
ref 7  设置为水中的氢
[ SOL_&_OW ]
201  204  207  210  213  216  ...

3 topol文件
; Created by Sobtop (http://sobereva.com/soft/sobtop) Version 1.0(dev3.1) on 2022-07-25

#include "amber99sb.ff/forcefield.itp"
#include "amber99sb.ff/tip3p.itp"
#include "L-citrulline.itp"

[ system ]
L-citrulline in water

[ molecules ]
; Molecule      nmols
L-citrulline     8
SOL              2086


后附整个md文件

ps1：虽然精氨酸rdf不对，但是溶液中的水分子的O-O的rdf是对着的，就感觉很奇怪


ps2：自己也换DMF或丙酮的水溶液跑过羰基氧与水分子氧的rdf，也不对劲，图基本和图1一样，g(r)没有高于1的峰

sobereva

把问题简化，先只统计一个分子上的氧和OW的rdf。也统计羧基的氧看看是什么情况

另外，当前模拟并不合理，在水中，此体系的羧基应该是质子解离状态，氨基应当是质子化的状态