粗粒化模拟在平衡过程中遇到问题

Wangyn

短肽分子通过charmm-gui得到粗粒化结构，选用Martini力场，通过packmol搭建溶液体系模型，力场文件用的是通过charmm-gui得到的，经过最小化，平衡过程无法正常进行。在进行平衡之前，输入export GMX_MAXCONSTRWARN=-1，避免因LINCS WARNING 而出错，然而程序却没有正常进行，log文件也只有一个开头，之前有问过，试过人们教我的方法，也有尝试修改mdp文件的参数，但都没有解决，应该是修改的参数作用不大或没有改到关键部分。由于都告诉我是初始结构或拓补文件的问题，所以请问我该怎么进行排查，谢谢。有尝试重新建模，意义不大。

log文件
Constraining the starting coordinates (step 0)

Constraining the coordinates at t0-dt (step 0)
RMS relative constraint deviation after constraining: 1.49e-06
Initial temperature: 300.364 K

Started mdrun on rank 0 Tue Jun  2 21:29:31 2020
           Step           Time         Lambda
              0        0.00000        0.00000

   Energies (kJ/mol)
           Bond       G96Angle  Improper Dih.        LJ (SR)   Coulomb (SR)
    1.34267e+03    1.12112e+03    2.19465e+02   -5.42810e+05   -2.35763e+03
Position Rest.      Potential    Kinetic En.   Total Energy    Temperature
    3.60034e+06    3.05786e+06    1.10291e+05    3.16815e+06    4.67923e+02
Pressure (bar)   Constr. rmsd
   -4.86563e+03    6.19741e-06


top文件：（水盒子边长20nm）
#include "toppar/martini_v2.2.itp"
#include "toppar/martini_v2.0_lipids_all_201506.itp"
#include "toppar/martini_v2.0_ions.itp"
#include "PROA_P.itp"

[ system ]
; name
Martini system

[ molecules ]
; name        number
PROA_P 200
W 16100
CL 200

sobereva

我没用CHARMM-GUI干过这种事，一般来说这种问题都是力场参数不合理、拓扑文件定义有bug、初始结构太烂导致的。mdp不适合当前模拟也可能导致。只能自己慢慢检查

Wangyn

sobereva 发表于 2020-5-21 08:03
我没用CHARMM-GUI干过这种事，一般来说这种问题都是力场参数不合理、拓扑文件定义有bug、初始结构太烂导致 ...

你好，我在用低核数的时候，界面没有变化，但命令是在运行中，用nohup命令，输出文件也没有变化，请问这有可能是什么原因？下面是平衡时的mdp文件，请帮忙看看有哪里需要整改。define                   =  -DPOSRES -DPOSRES_FC=4000
integrator               = md
tinit                    = 0.0
dt                       = 0.020
nsteps                   = 50000

nstlog                   = 1000
nstenergy                = 1000
nstxout-compressed       = 1000
compressed-x-precision   = 100

cutoff-scheme            = Verlet
nstlist                  = 20
ns_type                  = grid
pbc                      = xyz
verlet-buffer-tolerance  = 0.005

epsilon_r                = 15
coulombtype              = reaction-field
rcoulomb                 = 1.1
vdw_type                 = cutoff
vdw-modifier             = Potential-shift-verlet
rvdw                     = 1.1

tcoupl                   = v-rescale
tc-grps                  = protein solute
tau_t                    = 1.0 1.0
ref_t                    = 300 300

; Pressure coupling:
Pcoupl                   = berendsen
Pcoupltype               = isotropic
tau_p                    = 5.0
compressibility          = 4.5e-5
ref_p                    = 1.0

; GENERATE VELOCITIES FOR STARTUP RUN:
gen_vel                  = yes
gen_temp                 = 300
gen_seed                 = 8942579836
refcoord_scaling         = all

最后停止的界面会出现LINCX WARNING,如下
Step 9, time 0.18 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 1234.027907, max 6741.243164 (between atoms 16408 and 16406)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
  18862  18864   90.1   64.2889 240.4977      0.2700
  17319  17318   90.0    0.2699  48.5531      0.2700
  18862  18863   90.1   65.3206 241.1035      0.2700
  17318  17320  114.8    0.2690  23.2625      0.2700
  18863  18864   84.1    1.2228   1.2270      0.2700
  17319  17320   90.0    0.2710  48.0667      0.2700
  16415  16414   93.5    0.2951 334.8731      0.2700
  17152  17151  138.5   15.3910  14.5101      0.2700
  17144  17143   90.5    0.7166   0.3168      0.2700
  16414  16416   85.5    0.3036 233.2608      0.2700
  16407  16406   88.9    0.2702 1340.8149      0.2700
  17143  17142   90.4    0.7143   0.4197      0.2700
  17144  17142   90.5    0.7696   0.3190      0.2700
  16408  16406   89.8    0.2701 1820.4058      0.2700
  16408  16407   91.2    0.2697 820.5641      0.2700
  17147  17146   90.2    0.2731   0.4983      0.2700
  17148  17146   90.2    0.3534   0.4387      0.2700
  16415  16416   90.6    0.3800 371.9580      0.2700

Step 9, time 0.18 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 253.716985, max 2204.279297 (between atoms 18974 and 18976)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
  17151  17150   94.1   42.1391   3.1868      0.2700
  18974  18976   88.6  200.0600 595.4254      0.2700
  17152  17150   90.4   50.2689  13.6350      0.2700
  17148  17147   30.6    0.2697   0.2631      0.2700
  16767  16766   90.0    3.2043   8.6392      0.2700
  16768  16766   90.0    3.2251   7.2913      0.2700
  16768  16767   90.0    2.0290   3.7626      0.2700
  18974  18975   87.2  186.8602 503.4147      0.2700
  18975  18976   81.7   22.7970 115.4593      0.2700
  18859  18858   90.8    8.4864  18.3031      0.2700
  18860  18858   91.4  237.7959  28.4776      0.2700
  18860  18859   93.2  233.0337  21.7875      0.2700
  18863  18864   84.1    1.2228   1.2270      0.2700
  18863  18862   90.1   65.3206 241.1035      0.2700
  18862  18864   90.1   64.2889 240.4977      0.2700
  16428  16426   40.1    0.2719   0.2726      0.2700
  16203  16202   61.4    0.2790   0.2487      0.2700
  16204  16202   77.5    0.2880   0.2612      0.2700
  16204  16203   68.8    0.2696   0.3017      0.2700
Wrote pdb files with previous and current coordinates
Wrote pdb files with previous and current coordinates
Wrote pdb files with previous and current coordinates
Wrote pdb files with previous and current coordinates

sobereva

Wangyn 发表于 2020-5-21 13:40
你好，我在用低核数的时候，界面没有变化，但命令是在运行中，用nohup命令，输出文件也没有变化，请问这 ...

光从这些无从判断，而且不同粗粒化力场适合的mdp参数也不同

我也不知道你说的“界面没有变化”是什么

如果并行线程数少的时候能正常跑、结果也合理，先跑个诸如500 ps让体系弛豫一下，然后再用原本打算用的核数续跑。

Wangyn

sobereva 发表于 2020-5-23 02:35
光从这些无从判断，而且不同粗粒化力场适合的mdp参数也不同

我也不知道你说的“界面没有变化”是什么
...

谢谢老师的答复，我尝试了重新建模，模型是通过packmol搭建的，在跑平衡过程中虽然命令在运行中，但在nohup输出文件中经过一天多仍没有任何进展，仍停留在我上面截取的字段，所以程序并没有正常的进行下去，是不是与出现LINCS WARNING有关，应该怎么做呢

ulosggs

integrator 换成 sd, 步长减到 0.005，控压关掉。如果top文件里有constraints，暂时删掉。
这样跑一段时间弛豫之后，再把参数改回来。

Wangyn

ulosggs 发表于 2020-5-24 23:48
integrator 换成 sd, 步长减到 0.005，控压关掉。如果top文件里有constraints，暂时删掉。
这样跑一段时间 ...

你好，我做了相应更改，但仍有LINCX WARNING,错误如下：
Program gmx, VERSION 5.0.4
Source code file: /home/gromacs/bin/gromacs-5.0.4/src/gromacs/mdlib/constr.c, line: 224

Fatal error:
Too many LINCS warnings (1009)
If you know what you are doing you can adjust the lincs warning threshold in your mdp file
or set the environment variable GMX_MAXCONSTRWARN to -1,
but normally it is better to fix the problem
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors

Wangyn

sobereva 发表于 2020-5-23 02:35
光从这些无从判断，而且不同粗粒化力场适合的mdp参数也不同

我也不知道你说的“界面没有变化”是什么
...

老师您好，我做了一些更改，但程序仍然停滞不前。下面是经过1小时后log文件

There are: 19500 Atoms
Charge group distribution at step 0: 396 54 185 53 265 406 391 205 134 312 144 194 407 179 115 23 99 31 109 237 411 247 241 134 409 97 116 283 171 16 35 225 364 98 352 108 252 138 381 173 87 138 119 270 265 51 193 191 200 208 249 184 304 76 248 37 206 110 262 178 240 455 277 52 160 263 125 13 33 415 237 115 298 92 276 156 115 149 402 268 72 209 195 266 153 213 190 161 189 187 202 439 193 216 84 174 14 50 162 224

Constraining the starting coordinates (step 0)

Constraining the coordinates at t0-dt (step 0)
RMS relative constraint deviation after constraining: 1.70e-04
Initial temperature: 299.78 K

Started mdrun on rank 0 Wed May 27 16:46:43 2020
           Step           Time         Lambda
              0        0.00000        0.00000

   Energies (kJ/mol)
       G96Angle  Improper Dih.        LJ (SR)   Coulomb (SR) Position Rest.
    1.56191e+03    3.66693e+02   -4.98856e+05   -2.71839e+03    3.92034e+06
      Potential    Kinetic En.   Total Energy    Temperature Pressure (bar)
    3.42069e+06    1.04589e+05    3.52528e+06    4.58281e+02   -5.06149e+03
   Constr. rmsd
    3.80335e-04

naoki

可能是初始结构或拓扑文件有问题，先export GMX_MAXCONSTRWARN=-1然后看看跑起来是啥情况

Wangyn

naoki 发表于 2020-5-28 09:55
可能是初始结构或拓扑文件有问题，先export GMX_MAXCONSTRWARN=-1然后看看跑起来是啥情况

你好，我试了一下，可还是有问题，log文件及nohups输出文件始终没有进展，我不知道该怎样解决了
log文件
There are: 19500 Atoms
Charge group distribution at step 0: 393 95 219 129 419 395 361 122 68 316 190 237 337 139 202 42 166 109 128 214 199 204 219 189 263 98 155 197 111 51 69 221 309 101 273 166 285 171 347 212 293 307 180 210 247 103 222 152 232 287 252 243 189 61 317 145 215 90 215 79 183 261 179 139 181 183 118 107 84 377 223 165 300 195 216 131 225 130 253 186 106 215 191 347 207 145 236 155 258 202 90 284 151 188 59 61 90 116 155 248

Constraining the starting coordinates (step 0)

Constraining the coordinates at t0-dt (step 0)
RMS relative constraint deviation after constraining: 4.31e-04
Initial temperature: 300.409 K

Started mdrun on rank 0 Thu May 28 20:18:13 2020
           Step           Time         Lambda
              0        0.00000        0.00000

   Energies (kJ/mol)
       G96Angle  Improper Dih.        LJ (SR)   Coulomb (SR) Position Rest.
    1.50128e+03    2.36689e+02   -5.11048e+05   -2.53198e+03    2.60162e+06
      Potential    Kinetic En.   Total Energy    Temperature Pressure (bar)
    2.08978e+06    9.17413e+04    2.18152e+06    4.01985e+02   -3.43958e+03
   Constr. rmsd
    3.66937e-04

nohup文件
Command line:
  gmx mdrun -nt 100 -pin on -deffnm step4.2_equilibration -v


Back Off! I just backed up step4.2_equilibration.log to ./#step4.2_equilibration.log.1#
Reading file step4.2_equilibration.tpr, VERSION 5.0.4 (single precision)
Using 100 MPI threads
Using 1 OpenMP thread per tMPI thread
Setting the maximum number of constraint warnings to -1

Back Off! I just backed up step4.2_equilibration.xtc to ./#step4.2_equilibration.xtc.1#

Back Off! I just backed up step4.2_equilibration.edr to ./#step4.2_equilibration.edr.1#
starting mdrun 'Martini system'
50000 steps,   1000.0 ps.

step 0
top文件#include "toppar/martini_v2.2.itp"
#include "toppar/martini_v2.0_lipids_all_201506.itp"
#include "toppar/martini_v2.0_ions.itp"
#include "PROA_P.itp"

[ system ]
; name
Martini system

[ molecules ]
; name        number
W 16100
PROA_P 200
CL 200

我试过重新建模，我的体系是通过packmol搭建的，但并没有解决

yuhai

我的martini极化水模型模拟蛋白质也有同样的问题，就是报约束大于90度，蛋白质的两个原子键的真频小于了粗粒化常规步长20fs的5倍（注明如果小于20fs，用5fs是可以不报错运行，但是这不符合粗粒化的常规设置），想请问sob老师极化模型有什么特别吗，在非极化水模型是可以跑的。而且我前面跑另一个蛋白质的时候好像没有类似的错误，请sob老师和各位同学回复解答，万分感谢