请选择 进入手机版 | 继续访问电脑版
第八届北京科音分子动力学与GROMACS培训班将于2021年1月18~21日于北京举办,请点击此链接查看详情,这是系统性学习分子动力学模拟、掌握GROMACS程序使用的最好机会!(再下一届预计于3或4月举办)

计算化学公社

 找回密码
 现在注册!
查看: 308|回复: 8

[GROMACS] 求助,使用gromacs进行SDS水溶液中的模拟过程频繁报错

[复制链接]

125

帖子

0

威望

774

eV
积分
899

Level 4 (黑子)

发表于 2020-9-16 10:23:40 | 显示全部楼层 |阅读模式
本帖最后由 牧生 于 2020-9-16 11:38 编辑

1、首先使用M$画出了十二烷基硫酸钠   阴离子部分 的分子结构,另存为pdb文件,上传到http://bio2byte.be/acpype/,电荷由user指定为-1,获得一个压缩包,使用其中的SDS_GMX的SDS_GMX.gro,SDS_GMX.top,和SDS_GMX.itp三个文件


2、在5×5×5的盒子中加入5个SDS-,
gmx insert-molecules -ci SDS_GMX.gro -nmol 5 -box 5 5 5 -o SDS_box.gro
图形如下,看起来没错
图片1.png


3、向盒子中加入 tip4p的水
gmx solvate -cp SDS_box.gro -cs tip4p.gro -p SDS_GMX.top -o SDS_solv.gro

图片1.png

4、修改top文件  

之前的top文件如下:
; SDS_GMX.top created by acpype (v: 2019-07-10T18:04:00UTC) on Sun Sep 13 07:47:45 2020
[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               2               yes             0.5     0.8333

; Include SDS_GMX.itp topology
#include "SDS_GMX.itp"

[ system ]
SDS in water

[ molecules ]
; Compound        nmols
SDS              1     
SOL              4069

修改为:
; SDS_GMX.top created by acpype (v: 2019-07-10T18:04:00UTC) on Sun Sep 13 07:47:45 2020

;[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
;1               2               yes             0.5     0.8333                     
  (问题1,为何要注释这三行啊)

#include"amber99sb-ildn.ff/forcefield.itp" ; 加上这个力场
; Include SDS_GMX.itp topology
#include "SDS_GMX.itp"
#include"amber99sb-ildn.ff/tip4p.itp" ;加上水类型

[ system ]
SDS in water

[ molecules ]
; Compound        nmols
SDS              5  
         ;改成5个SDS
SOL              4069

5、能量最小化

em.mdp的内容如下:     (问题2:请帮忙看一下em.mdp有无设置错误,不合理的,有没有多的,或缺失的

define          = -DFLEXIBLE
integrator      = steep
emtol           = 250.0
nsteps          = 5000
nstenergy       = 1
energygrps      = System
nstlist         = 10
ns_type         = grid
coulombtype     = cut-off
rlist           = 1.0
rcoulomb        = 1.0
rvdw            = 1.0
constraints     = none
pbc             = xyz
cutoff-scheme= Verlet

命令   gmx grompp -f em.mdp -c SDS_solv.gro -p SDS_GMX.top -o ion.tpr -maxwarn 1     (忽略电荷为-1的错误)

6、添加离子,使得体系不带电
gmx genion -s ion.tpr -o SDS_solv.gro -neutral  -p SDS_GMX.top

选择替换SOL,添加5个Na+

7、修改top文件,将离子的部分力场添加进来

; SDS_GMX.top created by acpype (v: 2019-07-10T18:04:00UTC) on Sun Sep 13 07:47:45 2020


;[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
;1               2               yes             0.5     0.8333

#include "amber99sb-ildn.ff/forcefield.itp"
; Include SDS_GMX.itp topology
#include "SDS_GMX.itp"
#include "amber99sb-ildn.ff/tip4p.itp"
#include "amber99sb-ildn.ff/ions.itp
"        ; 添加 离子的力场      

[ system ]
SDS in water

[ molecules ]
; Compound        nmols
SDS              5     
SOL         4064
NA               5


8、再进行一次能量最小
gmx grompp -f em.mdp -c SDS_solv.gro -p SDS_GMX.top -o em-sol.tpr

gmx mdrun -v -deffnm em-sol


到这里都还算勉强能理解
以下部分就开始稀里糊涂了

9、最后成品模拟
md.mdp

integrator               = md
dt                       = 0.002
nsteps                   = 500000 ; 1 ns

nstxout                 = 500
nstvout                 = 500
nstfout                 = 500
nstenergy               = 500
nstlog                  = 500
;energygrps              = Protein Non-Protein

nstlist                  = 5
ns-type                  = Grid
pbc                      = xyz
rlist                    = 1.0

coulombtype              = PME
pme_order                = 4
fourierspacing           = 0.16
rcoulomb                 = 1.0
vdw-type                 = Cut-off
rvdw                     = 1.0

;Tcoupl                   = v-rescale
;tc-grps                  = Protein  Non-Protein
;tau_t                    = 0.1      0.1
;ref_t                    = 300      300     ;注释这四行,我觉得有错误,但是如果不注释,就不能继续进行,)

DispCorr                 = EnerPres

Pcoupl                   = Parrinello-Rahman
Pcoupltype               = Isotropic
tau_p                    = 2.0
compressibility          = 4.5e-5
ref_p                    = 1.0

gen_vel                  = no

constraints              = all-bonds
continuation             = yes
constraint_algorithm     = lincs
lincs_iter               = 1
lincs_order              = 4


命令为:
gmx grompp -f md.mdp -c SDS_solv.gro -p SDS_GMX.top -o npt-nopr.tpr

报了很多的note   


GROMACS:      gmx grompp, version 2019.3
Executable:   C:\gmx2019.3\bin\gmx.exe
Data prefix:  C:\gmx2019.3
Working dir:  F:\5SDS
Command line:
  gmx grompp -f md.mdp -c SDS_solv.gro -p SDS_GMX.top -o npt-nopr.tpr


NOTE 1 [file md.mdp, line 44]:
  md.mdp did not specify a value for the .mdp option "cutoff-scheme".
  Probably it was first intended for use with GROMACS before 4.6. In 4.6,
  the Verlet scheme was introduced, but the group scheme was still the
  default. The default is now the Verlet scheme, so you will observe
  different behaviour.


NOTE 2 [file md.mdp]:
  With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note
  that with the Verlet scheme, nstlist has no effect on the accuracy of
  your simulation.

Setting the LD random seed to -332505056
Generated 2556 of the 2556 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 0.5
Generated 2556 of the 2556 1-4 parameter combinations
Excluding 3 bonded neighbours molecule type 'SDS'
turning all bonds into constraints...
Excluding 2 bonded neighbours molecule type 'SOL'
turning all bonds into constraints...
Excluding 1 bonded neighbours molecule type 'NA'
turning all bonds into constraints...

NOTE 3 [file unknown]:
  You are using constraints on all bonds, whereas the forcefield has been
  parametrized only with constraints involving hydrogen atoms. We suggest
  using constraints = h-bonds instead, this will also improve performance.


NOTE 4 [file unknown]:
  For energy conservation with LINCS, lincs_iter should be 2 or larger.


Cleaning up constraints and constant bonded interactions with virtual sites
Removing all charge groups because cutoff-scheme=Verlet
Analysing residue names:
There are:     5      Other residues
There are:  4052      Water residues
There are:     5        Ion residues
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
Number of degrees of freedom in T-Coupling group rest is 24749.00

NOTE 5 [file md.mdp]:
  NVE simulation with an initial temperature of zero: will use a Verlet
  buffer of 10%. Check your energy drift!

Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 32x32x32, spacing 0.156 0.156 0.156
Estimate for the relative computational load of the PME mesh part: 0.11
This run will generate roughly 569 Mb of data

There were 5 notes



忽略这些note,再运行
gmx mdrun -deffnm npt-nopr -v

第一步,就出错了

  gmx mdrun -deffnm npt-nopr -v

Compiled SIMD: AVX_256, but for this host/run AVX2_256 might be better (see
log).
Reading file npt-nopr.tpr, VERSION 2019.3 (single precision)
Can not increase nstlist because an NVE ensemble is used
Using 1 MPI thread
Using 8 OpenMP threads

starting mdrun 'SDS in water'
500000 steps,   1000.0 ps.

Step 0, time 0 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.261971, max 2.839417 (between atoms 96 and 121)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
      6     25   90.0    0.1141   0.1675      0.1097
     59     82   90.0    0.1146   0.3110      0.1097
     59     84   90.0    0.1142   0.1041      0.1097
     94    117   90.0    0.1131   0.2503      0.1097
     96    121   90.0    0.1149   0.4211      0.1097
    176    196   90.0    0.1141   0.2024      0.1097
    177    199   41.3    0.1132   0.1095      0.1097
Wrote pdb files with previous and current coordinates
step 0
Step 1, time 0.002 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 21799.666016, max 220749.859375 (between atoms 177 and 199)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
      5     22   40.0    0.1097   0.1068      0.1097
      5     23   40.3    0.1097   0.1061      0.1097
      6      7   46.9    0.1532   0.1609      0.1538
      6     25   89.9    0.1675   2.1332      0.1097
      8     28   90.0    0.1097   2.6021      0.1097
     52     53   90.0    0.1538   0.3086      0.1538
     53     54   90.0    0.1538   0.3315      0.1538
     53     76   90.0    0.1097   0.3326      0.1097
     53     77   90.0    0.1097   2.1089      0.1097
     54     78   34.8    0.1097   0.1324      0.1097
     54     79   36.6    0.1097   0.1370      0.1097
     59     82  155.0    0.3110   4.3932      0.1097
     59     84   89.9    0.1041   0.2050      0.1097
     92     93   67.8    0.1538   0.1849      0.1538
     93     94   90.4    0.1536   0.2252      0.1538
     93    114   77.1    0.1097  27.0166      0.1097
     93    115   90.1    0.1094   0.1766      0.1097
     94     95   33.8    0.1537   0.1826      0.1538
     94    116   31.8    0.1096   0.1270      0.1097
     94    117   53.2    0.2503  25.8629      0.1097
     96    121   90.0    0.4211   0.2488      0.1097
    175    176   44.7    0.1535   0.2061      0.1538
    176    177   89.5    0.1535   0.6377      0.1538
    176    196   97.6    0.2024 24204.2461      0.1097
    176    197   53.8    0.1093   0.1665      0.1097
    177    178   89.9    0.1537   0.4125      0.1538
    177    198   89.7    0.1097   0.5452      0.1097
    177    199   37.5    0.1095 24214.1621      0.1097
    178    179   89.7    0.1538   0.5984      0.1538
    178    200   89.7    0.1097   0.5567      0.1097
    178    201   89.5    0.1097   0.2780      0.1097
    179    180   34.5    0.1538   0.1872      0.1538
    179    202   41.3    0.1097   0.1466      0.1097
    179    203   41.4    0.1097   0.1468      0.1097
Wrote pdb files with previous and current coordinates





请帮忙看一下,哪里的设置,或者操作步骤,或者命令有问题。

131

帖子

0

威望

1107

eV
积分
1238

Level 4 (黑子)

Nerv

发表于 2020-9-16 11:50:56 | 显示全部楼层
体系崩溃一般是结构不合理或者力场有问题,不要用constraints              = all-bonds,用hbonds,先关掉控压,把控温开启,跑nvt,步长设为0.001,其余参数用默认的即可,不了解的不要随意设置。
God's in his heaven,all rights with the world

131

帖子

0

威望

1107

eV
积分
1238

Level 4 (黑子)

Nerv

发表于 2020-9-16 11:52:43 | 显示全部楼层
问题1:你注释的三行在#include "amber99sb-ildn.ff/forcefield.itp"中已经出现过了,不能重复出现
God's in his heaven,all rights with the world

6

帖子

0

威望

538

eV
积分
544

Level 4 (黑子)

发表于 2020-11-27 23:14:47 | 显示全部楼层
你发这么长是想让大佬给你来个基础培训班么

16

帖子

0

威望

37

eV
积分
53

Level 2 能力者

发表于 4 天前 | 显示全部楼层
请问一下楼主解决这个问题了吗?我前面额外做了控温控压力的平衡,然后进行的模拟,也是出现了一对警告,类似relative constraint deviation after LINCS这种的,楼主是怎么解决的?

125

帖子

0

威望

774

eV
积分
899

Level 4 (黑子)

 楼主| 发表于 3 天前 | 显示全部楼层
黑择明 发表于 2021-1-12 20:40
请问一下楼主解决这个问题了吗?我前面额外做了控温控压力的平衡,然后进行的模拟,也是出现了一对警告,类 ...

解决了。如果你的模拟有警告,试试做两次em

16

帖子

0

威望

37

eV
积分
53

Level 2 能力者

发表于 3 天前 | 显示全部楼层
牧生 发表于 2021-1-13 07:37
解决了。如果你的模拟有警告,试试做两次em

是重复做两次能量最小化吗?楼主最小化里integrator用的是steep还是cg?

16

帖子

0

威望

37

eV
积分
53

Level 2 能力者

发表于 3 天前 | 显示全部楼层
牧生 发表于 2021-1-13 07:37
解决了。如果你的模拟有警告,试试做两次em

我看楼主一开始发的问题,先做了两次能量最小化,结果还是出了一堆警告?

125

帖子

0

威望

774

eV
积分
899

Level 4 (黑子)

 楼主| 发表于 3 天前 | 显示全部楼层
能量最小化当然是用steep

出现警告是别的原因,自己多检查一下mdp文件,网上也有别的很多可以参考的文件,最终是可以解决问题的。
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949-1号 )

GMT+8, 2021-1-16 19:57 , Processed in 0.169335 second(s), 27 queries .

快速回复 返回顶部 返回列表