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

牧生

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

图形如下，看起来没错

3、向盒子中加入 tip4p的水

gmx solvate -cp SDS_box.gro -cs tip4p.gro -p SDS_GMX.top -o SDS_solv.gro

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

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

Lacrimosa

体系崩溃一般是结构不合理或者力场有问题，不要用constraints              = all-bonds，用hbonds，先关掉控压，把控温开启，跑nvt，步长设为0.001，其余参数用默认的即可，不了解的不要随意设置。

Lacrimosa

问题1：你注释的三行在#include "amber99sb-ildn.ff/forcefield.itp"中已经出现过了，不能重复出现

nianbin

你发这么长是想让大佬给你来个基础培训班么

黑择明

请问一下楼主解决这个问题了吗？我前面额外做了控温控压力的平衡，然后进行的模拟，也是出现了一对警告，类似relative constraint deviation after LINCS这种的，楼主是怎么解决的？

牧生

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

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

黑择明

牧生 发表于 2021-1-13 07:37
解决了。如果你的模拟有警告，试试做两次em

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

黑择明

牧生 发表于 2021-1-13 07:37
解决了。如果你的模拟有警告，试试做两次em

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

牧生

能量最小化当然是用steep

出现警告是别的原因，自己多检查一下mdp文件，网上也有别的很多可以参考的文件，最终是可以解决问题的。