请选择 进入手机版 | 继续访问电脑版
第13届北京科音初级量子化学培训班将于10月5~8日于北京举办,请点击此链接查看详情。这是新人一次性正确、完整学习量子化学计算的最好、最快机会,能少走无数弯路,欢迎参加并相互转告!(已报满)

计算化学公社

 找回密码
 现在注册!
查看: 54|回复: 2

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

[复制链接]

65

帖子

0

威望

531

eV
积分
596

Level 4 (黑子)

发表于 4 天前 | 显示全部楼层 |阅读模式
本帖最后由 牧生 于 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





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

87

帖子

0

威望

874

eV
积分
961

Level 4 (黑子)

Nerv

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

87

帖子

0

威望

874

eV
积分
961

Level 4 (黑子)

Nerv

发表于 4 天前 | 显示全部楼层
问题1:你注释的三行在#include "amber99sb-ildn.ff/forcefield.itp"中已经出现过了,不能重复出现
God's in his heaven,all rights with the world
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

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

GMT+8, 2020-9-20 22:20 , Processed in 0.242031 second(s), 27 queries .

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