计算化学公社

标题: 求助，关于SETTLE算法的一个疑问 [打印本页]

作者
Author: weng 时间: 2021-4-7 14:01
标题: 求助，关于SETTLE算法的一个疑问
各位老师好：    我目前尝试用gromacs进行TIP3P水分子模型的MD模拟。体系中共有1728个水分子。我明确地在mdp文件中设置了constraints = all-bonds，但我发现得到的轨迹当中水分子的O-H键长（图1）和H-H间距（图2）存在明显的波动，波动幅度大致为键长的1%。

(, 下载次数 Times of downloads: 19)
图1：某个水分子的O-H键长随时间的变化（1 frame = 0.2 ns）

(, 下载次数 Times of downloads: 23)
图2：某个水分子的H-H间距随时间的变化（1 frame = 0.2 ns）

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

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

(, 下载次数 Times of downloads: 4)

(, 下载次数 Times of downloads: 2)

(, 下载次数 Times of downloads: 0)

附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

作者
Author: sobereva 时间: 2021-4-8 03:04
根本没必要写constraints = all-bonds，你不指定-DFLEXIBLE直接就是用SETTLS约束水为刚性
all-bonds是对于其它情况靠LINCS约束键长才用的

作者
Author: weng 时间: 2021-4-12 13:32

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

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

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