求助，关于SETTLE算法的一个疑问

weng

各位老师好：      我目前尝试用gromacs进行TIP3P水分子模型的MD模拟。体系中共有1728个水分子。我明确地在mdp文件中设置了constraints = all-bonds，但我发现得到的轨迹当中水分子的O-H键长（图1）和H-H间距（图2）存在明显的波动，波动幅度大致为键长的1%。


图1：某个水分子的O-H键长随时间的变化（1 frame = 0.2 ns）


图2：某个水分子的H-H间距随时间的变化（1 frame = 0.2 ns）

      我用的是gromacs 5.1.2，4 cpus。但在gtx 1060上运行gromacs 2019也有类似的问题。请问大家有没有碰到过类似的情况？应该怎么解决？

      附件是相关的部分文件。
equ.cpt (123.19 KB, 下载次数 Times of downloads: 0)

2021-4-7 14:00 上传 Uploaded

点击下载
Click to download

topol.top (1.26 KB, 下载次数 Times of downloads: 4)

2021-4-7 13:59 上传 Uploaded

点击下载
Click to download

npt.mdp (3.03 KB, 下载次数 Times of downloads: 2)

2021-4-7 13:59 上传 Uploaded

点击下载
Click to download

equ.gro (349.36 KB, 下载次数 Times of downloads: 0)

2021-4-7 13:58 上传 Uploaded

点击下载
Click to download

附1：我的mdp文件内容如下：
include                  =
; e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)
define                   =

; RUN CONTROL PARAMETERS
integrator               = md
tinit                    = 0
dt                       = 0.002
nsteps                   = 250000
init-step                = 0
comm-mode                = Linear
nstcomm                  = 10
comm-grps                = System

; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout                  = 10000
nstvout                  = 0
nstfout                  = 0
nstlog                   = 10000
nstcalcenergy            = 100
nstenergy                = 0
nstxout-compressed       = 100
compressed-x-precision   = 1000
compressed-x-grps        = system
energygrps               =

; NEIGHBORSEARCHING PARAMETERS
; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
cutoff-scheme            = Verlet
nstlist                  = 10
ns_type                  = grid
pbc                      = xyz
; nblist cut-off        
rlist                    = 1.2
; long-range cut-off for switched potentials

coulombtype              = PME
rcoulomb                 = 1.2
; Method for doing Van der Waals
vdwtype                  = switch
; cut-off lengths      
rvdw-switch              = 1.0
rvdw                     = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                 = EnerPres

; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.14
pme_order                = 4

; Temperature coupling  
tcoupl                   = nose-hoover
tc-grps                  = Water
; Time constant (ps) and reference temperature (K)
tau_t                    = 0.5
ref_t                    = 298
; pressure coupling     
pcoupl                   = Parrinello-Rahman
pcoupltype               = isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau_p                    = 1.0
compressibility          = 4.5e-5
ref_p                    = 1.0
; Scaling of reference coordinates, No, All or COM
refcoord-scaling         = No

; OPTIONS FOR BONDS   
constraints              = all-bonds
constraint-algorithm     = Lincs
continuation             = yes
lincs-order              = 4
lincs-iter               = 1
lincs-warnangle          = 30

; GENERATE VELOCITIES FOR STARTUP RUN
;gen_vel                  = yes
;gen_temp                 = 300
;gen_seed                 = 456820



附2：我使用的水分子参数来自charmm36m的tip3p.itp，用的topol.top如下
[ defaults ]
; nbfunc        comb-rule        gen-pairs        fudgeLJ        fudgeQQ
1        2        yes        1.0        1.0

[ atomtypes ]
;!!!THIS ATOMTYPES DEFINITION ARE ARBITRARY AND STATED HERE TO COMPLY GROMACS REQUIREMENTS
; THEY WILL BE OVERDEFINED BY ATOMS SECTION BELOW!!!
;name  at.num      mass        charge   ptype               c6      c12
OT      8       15.999400       0.000      A    0.315057422683  0.63639
HT      1        1.00800        0.000      A    0.00            0.00000

[ moleculetype ]
; molname           nrexcl
WAT                2

[ atoms ]
; id   at type        res nr        residu name at name                cg nr        charge
1                OT                1                SOL                                 OW                                1                -0.834
2                HT                1                SOL                                HW1                                1                 0.417
3                HT                1                SOL                                HW2                                1                 0.417

#ifdef FLEXIBLE

#ifdef ORIGINAL_TIP3P
[ bonds ]
; i j        funct        length        force.c.
1        2        1        0.09572 502416.0 0.09572        502416.0
1        3        1        0.09572 502416.0 0.09572        502416.0

[ angles ]
; i  j        k        funct        angle        force.c.
2         1        3        1        104.52        628.02        104.52        628.02       

#else
;CHARMM TIP3p
[ bonds ]
; i j        funct        length        force.c.
1        2        1        0.09572 376560.0 0.09572        376560.0
1        3        1        0.09572 376560.0 0.09572        376560.0

[ angles ]
; i  j        k        funct        angle        force.c.
2         1        3        1        104.52        460.24        104.52        460.24       
#endif


#else
[ settles ]
; i   j funct        length
1          1 0.09572 0.15139

[ exclusions ]
1 2 3
2 1 3
3 1 2
#endif

[ system ]
charm36 TIP3P

[ molecules ]
WAT  1728

sobereva

根本没必要写constraints = all-bonds，你不指定-DFLEXIBLE直接就是用SETTLS约束水为刚性
all-bonds是对于其它情况靠LINCS约束键长才用的

weng

sobereva 发表于 2021-4-8 03:04
根本没必要写constraints = all-bonds，你不指定-DFLEXIBLE直接就是用SETTLS约束水为刚性
all-bonds是对于 ...

感谢老师！
我后来弄明白了，原来是存储精度的问题。我设置的xtc存储精度不够，所以四舍五入之后，键长就会有比较大的变化。