计算化学公社

标题: 甲酸分子的模拟问题 [打印本页]

作者
Author: tingguang 时间: 2019-3-22 09:05
标题: 甲酸分子的模拟问题
各位老师好！

最近希望模拟甲酸作为体系溶剂的情况。首先我准备建立一个装有100个甲酸分子的标准溶剂盒子以方便gromacs进行后续溶剂化。甲酸分子采用网上提供的GAFF力场参数。一直到模拟NVT时没出现问题。当进行NPT时问题不断。由于最初的100个分子的结构用packmol软件生成，因此盒子难免较大，模拟过程中，甲酸分子凝聚在一起，出现真空，因此报错。于是我将凝聚在一起的结构提取出来，重新用editconf命令加上盒子信息。在NPT模拟进行到某个时刻时出现错误“The X-size of the box times triclinic skew factor”。查询网络说是因为盒子塌了。于是根据网上建议，将时间步长改为1fs，压力耦合采用Berendson方法，压力耦合常数tau_p增加至4. 温度耦合方法采用V-rescale方法。仍解决不了问题。采用本站版主的建议，忽略NVT步骤，仍存在问题。是不是因为甲酸分子间氢键相互作用强而发生较强缔合造成的？我模拟其他溶剂分子，譬如乙醇、三氯甲烷都没发现这种情况。

我实在没招了，有同仁经历过类似问题吗？任何建议都表示感谢！

作者
Author: sobereva 时间: 2019-3-22 12:32
甲酸肯定会聚集在一起，但绝对没有理由出现真空，NPT的时候自发就消掉了。要么是初始结构有严重问题，要么是计算方式或mdp有问题，要么力场有问题
做NVT模拟没有丝毫意义，一开始packmol给的结构必然很松散，NVT下必然会产生大量真空区，这种无意义的模拟纯粹是浪费时间

作者
Author: tingguang 时间: 2019-3-22 13:29

sobereva 发表于 2019-3-22 12:32
甲酸肯定会聚集在一起，但绝对没有理由出现真空，NPT的时候自发就消掉了。要么是初始结构有严重问题，要么 ...

已经重复各种计算方案多次了。每当相对密度在100多时就会出现一下错误：“Program gmx mdrun, VERSION 5.1.4
Source code file: /home/chengqianwei/program/gromacs-5.1.4/src/gromacs/domdec/domdec.cpp, line: 3113

Fatal error:
The X-size of the box (2.537634) times the triclinic skew factor (1.000000) is smaller than the number of DD cells (4) times the smallest allowed cell size (0.634346)

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors”。力场参数是网站上获得的GAFF力场参数。我后面采用你的方案，直接进行NPT模拟了。都是同样问题，即当盒子缩小到一定程度时计算崩溃。

作者
Author: tingguang 时间: 2019-3-22 13:30
本帖最后由 tingguang 于 2019-3-22 13:44 编辑

这是我的mdp文件：

integrator             = md
nsteps                = 10000000
dt                   = 0.001
; Output control
nstxout                = 5000
nstvout                = 5000
nstenergy             = 5000
nstlog                = 5000
; Bond parameters
continuation          = yes
constraint_algorithm = lincs
constraints          = h-bonds
lincs_iter             = 1       ; accuracy of LINCS
lincs_order          = 4       ; also related to accuracy
; Nonbonded settings
cutoff-scheme          = Verlet ; Buffered neighbor searching
ns_type                = grid    ; search neighboring grid cells
nstlist                = 10       ; 20 fs, largely irrelevant with Verlet
rcoulomb             = 1.0    ; short-range electrostatic cutoff (in nm)
rvdw                   = 1.0    ; short-range van der Waals cutoff (in nm)
DispCorr             = EnerPres  ; account for cut-off vdW scheme
; Electrostatics
coulombtype          = PME    ; Particle Mesh Ewald for long-range electrostatics
pme_order             = 4       ; cubic interpolation
fourierspacing       = 0.16    ; grid spacing for FFT
; Temperature coupling is on
tcoupl                = V-rescale       ; modified Berendsen thermostat
tc-grps                = LIG ; two coupling groups - more accurate
tau_t                = 0.1             ; time constant, in ps
ref_t                = 300             ; reference temperature, one for each group, in K

; Pressure coupling is on
pcoupl                = berendsen    ; Pressure coupling on in NPT
pcoupltype             = isotropic          ; uniform scaling of box vectors
tau_p                = 4.0                ; time constant, in ps
ref_p                = 1.0                ; reference pressure, in bar
compressibility       = 4.5e-5             ; isothermal compressibility of water, bar^-1
refcoord_scaling       = com
; Periodic boundary conditions
pbc                   = xyz    ; 3-D PBC
; Velocity generation
gen_vel                = no    ; assign velocities from Maxwell distribution
gen_temp             = 300    ; temperature for Maxwell distribution
gen_seed             = -1       ; generate a random seed

作者
Author: tingguang 时间: 2019-3-22 13:36
本帖最后由 tingguang 于 2019-3-22 13:44 编辑

D:\vmdscene.bmp这是模拟崩溃时的图片。我发现模拟的分子并未在Z轴方向上移动，其它方向已经很密实了，而Z轴方向上盒子还很空，这是什么原因了，版主请指点。

作者
Author: sobereva 时间: 2019-3-22 19:14
直接用acpype得到分子的基于GAFF的拓扑文件多好
如果别的地方搞来的拓扑文件里的原子顺序与结构文件里的不符也会导致模拟彻底乱套

作者
Author: tingguang 时间: 2019-3-22 22:53
本帖最后由 tingguang 于 2019-3-22 22:58 编辑

sobereva 发表于 2019-3-22 19:14
直接用acpype得到分子的基于GAFF的拓扑文件多好
如果别的地方搞来的拓扑文件里的原子顺序与结构文件里的不 ...

好的，我试试，谢谢！另外，这个图尽管有点乱，这是周期性边界条件造成的，其实甲酸的分子结构并没有变形。这些分子在盒子的Z轴方向上并没有分散，您觉得是什么原因呢？

作者
Author: tingguang 时间: 2019-3-26 20:45
本帖最后由 tingguang 于 2019-3-26 20:52 编辑

sobereva 发表于 2019-3-22 19:14
直接用acpype得到分子的基于GAFF的拓扑文件多好
如果别的地方搞来的拓扑文件里的原子顺序与结构文件里的不 ...

sobereva您好！根据您的建议，我自己用acpype等工具生成了甲酸分子的GAFF拓扑，并进行分子动力学模拟。然而无论如何改参数，仍旧是同样问题：随着盒子逐步缩小，甲酸分子并未均匀扩散，而是保持凝缩在一起。最后，在x、y、z轴其中一个方向盒子接触到甲酸分子，gromacs随即报错：“the X-size of the box (2.641884) times the triclinic skew factor (1.000000) is smaller than the number of DD cells (4) times the smallest allowed cell size”（见附图和我用的各参数文件）。麻烦老师你指点！

作者
Author: cokoy 时间: 2019-3-26 22:06
你盒子太大了，NPT必然会使盒子减小啊。。而且甲酸之间有氢键肯定也会聚集

作者
Author: tingguang 时间: 2019-3-26 22:13
本帖最后由 tingguang 于 2019-3-26 22:15 编辑

cokoy 发表于 2019-3-26 22:06
你盒子太大了，NPT必然会使盒子减小啊。。而且甲酸之间有氢键肯定也会聚集

盒子大小没什么关系。我尝试过N个大小的盒子

。无论什么盒子，当盒子一边接触到甲酸分子时就报错。以前模拟过甲醇、乙醇和氯仿，都没这种问题。甲酸出现的这个问题真不知道怎么处理！关键是，甲酸分子不均匀扩散，行为类似粘滞流体一般。

作者
Author: sobereva 时间: 2019-3-27 04:26

tingguang 发表于 2019-3-26 22:13
盒子大小没什么关系。我尝试过N个大小的盒子。无论什么盒子，当盒子一边接触到甲酸分子时就报错。以前 ...

你这个结构明显不合理。
盒子那么大，真空区一大片，做NPT必然会导致盒子迅速收缩，还容易因此出现模拟不稳定
如果你要模拟凝聚相，一开始就应该把盒子填满

作者
Author: tingguang 时间: 2019-3-27 07:31

sobereva 发表于 2019-3-27 04:26
你这个结构明显不合理。
盒子那么大，真空区一大片，做NPT必然会导致盒子迅速收缩，还容易因此出现模拟 ...

我调整过盒子的大小。当盒子小时，会出现lincs错误的问题

作者
Author: sobereva 时间: 2019-3-27 10:26

tingguang 发表于 2019-3-27 07:31
我调整过盒子的大小。当盒子小时，会出现lincs错误的问题

这说明你当前的mdp或者拓扑文件本身有问题

作者
Author: tingguang 时间: 2019-3-27 19:23

sobereva 发表于 2019-3-27 10:26
这说明你当前的mdp或者拓扑文件本身有问题

mdp是常用的mdp。拓扑文件我仔细看过，自己做的和网上的差不多，仅有数值上的稍许偏差。

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