GROMACS中纤维二糖分子做NVT控温过程体系崩溃求助(已解决)

gsbear

在GROMOS53A6Carbo力场下，做了一个分子的纤维二糖加了1013分子的水，之后做了能量极小化，emtol设置为100；接下来跑NVT控温过程，老是出LINCS崩溃，其中都是葡萄糖上各个碳原子所连接羟基的氧原子和氢原子间键旋转超30度的报警，部分出错报警信息如下：
Step 360, time 0.36 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.024789, max 0.070114 (between atoms 30 and 31)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
     30     31   90.0    0.1005   0.1070      0.1000
     30     31   90.0    0.1005   0.1070      0.1000
     30     31   90.0    0.1005   0.1070      0.1000
     30     31   90.0    0.1005   0.1070      0.1000
     30     31   90.0    0.1005   0.1070      0.1000

Step 361, time 0.361 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000001, max 0.000001 (between atoms 2 and 3)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
     30     31   78.8    0.1070   0.1000      0.1000
     30     31   78.8    0.1070   0.1000      0.1000
     30     31   78.8    0.1070   0.1000      0.1000
     30     31   78.8    0.1070   0.1000      0.1000
     30     31   78.8    0.1070   0.1000      0.1000

Step 362, time 0.362 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.021378, max 0.060467 (between atoms 30 and 31)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
     30     31   90.0    0.1000   0.1060      0.1000
     30     31   90.0    0.1000   0.1060      0.1000
     30     31   90.0    0.1000   0.1060      0.1000
     30     31   90.0    0.1000   0.1060      0.1000
     30     31   90.0    0.1000   0.1060      0.1000


------------------------------------------------------------------------------
看了GROMACS手册中关于Blow up的原因和解决方法，也查了网上的解决办法，尝试了如下方法都无法解决：
1、进一步做能量最小化（将之前emtol值从1000改为100）
2、约束算法换LINCS算法为SHAKE算法
3、减小步长从2fs减小到1fs（好像不能在小了吧？）
不知道是否还有其他的解决办法可以处理这个问题，望给位老师不吝赐教，十分感谢！
上传了相关的.gro .mdp .top文件在附件中可供参阅

我是猫猫大王

gsbear 发表于 2020-12-15 09:42
估计你可能是没有把GROMOS 56A6_CARBO的力场文件正确解压到GROMACS的力场目录中，或者是在pdb2gmx的时候 ...

您好，按照您的方法我使用pdb2gmx生成多糖top文件的时候出现了报错，错误“Atom 4 in residue GLC 0 was not found in rtp entry GLC with 14 atoms while sorting atoms.” 请问有什么解决办法呢 谢谢答复！

zdworld

多糖出问题首先把末端糖环4号羟基氢与5号氧原子之间的pair去掉试试。

gsbear

从GMML的GITHUB页面下载压缩包或者自己git以后编译，只是编译需要用Ubuntu的linux系统

如果你

您好，想请问一下您GMML这个软件怎么下载安装呢

gsbear

huhuhutian 发表于 2021-9-10 10:33
想问一下 在GROMACS中怎么处理糖基化蛋白质

目前没有糖蛋白专用力场，估计你得用全原子力场建模

huhuhutian

想问一下 在GROMACS中怎么处理糖基化蛋白质

gsbear

juddtrump 发表于 2020-12-14 19:12
您好！对照你讲的进行了学习，但是用gmx pdb2gmx生成纤维二塘的.top文件报错，错误“Residue '' not found  ...

估计你可能是没有把GROMOS 56A6_CARBO的力场文件正确解压到GROMACS的力场目录中，或者是在pdb2gmx的时候没有选GROMOS 56A6_CARBO的力场。

如果上面两个问题都是正确的，那应该是.pdb文件没有按照力场的规定来划分残基。
GROMOS 56A6_CARBO只支持氧化端、还原端和中间端三种1-4贝塔葡萄糖残基，并且残基中原子命名也一定要跟力场的规范一致，具体的可以参考力场文件中附带的例子那个1-4trimer.pdb，内容如下：

REMARK    GENERATED BY TRJCONV
REMARK    THIS IS A SIMULATION BOX
CRYST1   30.037   30.037   30.037  90.00  90.00  90.00 P 1           1
MODEL        1
ATOM      1  C4  GLC0    1      17.560  10.430  18.280  1.00  0.00            
ATOM      2  O4  GLC0    1      18.450   9.690  19.130  1.00  0.00            
ATOM      3  HO4 GLC0    1      18.130   8.740  19.180  1.00  0.00            
ATOM      4  C3  GLC0    1      16.350  10.930  19.040  1.00  0.00            
ATOM      5  O3  GLC0    1      15.710  10.070  20.000  1.00  0.00            
ATOM      6  HO3 GLC0    1      16.410   9.860  20.680  1.00  0.00            
ATOM      7  C2  GLC0    1      15.400  11.380  17.940  1.00  0.00            
ATOM      8  O2  GLC0    1      14.140  11.650  18.590  1.00  0.00            
ATOM      9  HO2 GLC0    1      13.360  12.010  18.080  1.00  0.00            
ATOM     10  C6  GLC0    1      19.530  11.200  16.870  1.00  0.00            
ATOM     11  O6  GLC0    1      19.940  12.020  15.770  1.00  0.00            
ATOM     12  HO6 GLC0    1      19.340  11.800  15.000  1.00  0.00            
ATOM     13  C5  GLC0    1      18.260  11.600  17.600  1.00  0.00            
ATOM     14  O5  GLC0    1      17.300  12.260  16.760  1.00  0.00            
ATOM     15  C1  GLC0    1      16.100  12.640  17.450  1.00  0.00            
ATOM     16  O1  GLC0    1      15.240  13.340  16.550  1.00  0.00            
ATOM     17  C4  GLC     2      15.730  14.500  15.860  1.00  0.00            
ATOM     18  C3  GLC     2      15.560  14.220  14.370  1.00  0.00            
ATOM     19  O3  GLC     2      16.470  13.170  14.010  1.00  0.00            
ATOM     20  HO3 GLC     2      16.760  12.710  14.850  1.00  0.00            
ATOM     21  C2  GLC     2      15.880  15.460  13.560  1.00  0.00            
ATOM     22  O2  GLC     2      15.860  15.330  12.130  1.00  0.00            
ATOM     23  HO2 GLC     2      16.350  14.480  11.910  1.00  0.00            
ATOM     24  C6  GLC     2      14.750  16.090  17.570  1.00  0.00            
ATOM     25  O6  GLC     2      15.960  16.480  18.230  1.00  0.00            
ATOM     26  HO6 GLC     2      16.150  17.460  18.130  1.00  0.00            
ATOM     27  C5  GLC     2      14.830  15.700  16.100  1.00  0.00            
ATOM     28  O5  GLC     2      15.320  16.880  15.460  1.00  0.00            
ATOM     29  C1  GLC     2      15.010  16.590  14.090  1.00  0.00            
ATOM     30  O1  GLC     2      15.080  17.790  13.310  1.00  0.00            
ATOM     31  C4  GLCN    3      13.900  18.560  13.560  1.00  0.00            
ATOM     32  C3  GLCN    3      14.320  19.930  14.070  1.00  0.00            
ATOM     33  O3  GLCN    3      15.330  19.920  15.090  1.00  0.00            
ATOM     34  HO3 GLCN    3      16.060  19.320  14.780  1.00  0.00            
ATOM     35  C2  GLCN    3      13.050  20.680  14.460  1.00  0.00            
ATOM     36  O2  GLCN    3      13.490  21.830  15.190  1.00  0.00            
ATOM     37  HO2 GLCN    3      14.110  21.520  15.910  1.00  0.00            
ATOM     38  C6  GLCN    3      12.770  17.960  11.320  1.00  0.00            
ATOM     39  O6  GLCN    3      11.860  18.320  10.270  1.00  0.00            
ATOM     40  HO6 GLCN    3      11.050  18.780  10.640  1.00  0.00            
ATOM     41  C5  GLCN    3      13.050  19.010  12.380  1.00  0.00            
ATOM     42  O5  GLCN    3      11.790  19.570  12.780  1.00  0.00            
ATOM     43  C1  GLCN    3      12.180  20.870  13.230  1.00  0.00            
ATOM     44  O1  GLCN    3      10.980  21.640  13.390  1.00  0.00            
ATOM     45  HO1 GLCN    3      11.190  22.610  13.480  1.00  0.00            

juddtrump

您好！对照你讲的进行了学习，但是用gmx pdb2gmx生成纤维二塘的.top文件报错，错误“Residue '' not found in residue topology database”，请问该怎么解决呢？谢谢

gsbear

问题终于解决了，困扰了我大半年，这回五一在家有空又重新梳理了一下，发现了自己之前工作的两个疏漏，修正后终于把这个跑NVT崩溃的问题解决了，现在已经跑完NVT和NPT，开始跑Production MD了。等着第1个纳秒MD结果的档口，上来把解决的过程记录一下，如果也有人遇到相同问题可供参考。
主要的两点：
1、之前把GROMOS 56A6_CARBO力场和GROMOS 56A6_CARBO_R这两个东西混淆起来了，GROMOS 56A6_CARBO是发表再JCC2011年杂志上的，而GROMOS 56A6_CARBO_R是对它的一个修改，发表在JCC2016年杂志上（详细见帖子末尾参考文献）。GROMOS 56A6_CARBO力场文件下载地址是http://www.gromacs.org/@api/deki/files/200/=56a_CARBO4GROMACS.rar
2、GROMOS 56A6_CARBO力场的文件下载后里边有一个python写的程序CONVERT.py这个程序需要放在和topol.top相同目录下运行，程序会给gmx pdb2gmx所产生的拓扑文件打补丁，生成新的一个拓扑文件top_new.top，用这个新的拓扑文件跑就不会NVT崩溃的问题了（Sob老师判断的拓扑文件有问题真准！）

这里要补充说明的问题是力场所带的这个拓扑补丁程序是用python2写的，python3下面运行会出错，之前因为不熟悉python语言，也没有注意到这个问题，直接用pdb2gmx产生的有问题的拓扑文件去跑控温，所以崩溃了。
这个CONVERT.py程序要用python3运行的话要改两个地方：
112行：if语句第一个条件nvW6.has_key(atom1_name+atom2_name)改为atom1_name+atom2_name in nvW6.keys()
117行：nvW6.has_key(atom2_name+atom1_name)改为atom2_name+atom1_name in nvW6.keys()
就可以了，主要原因是python3里边dict类没有了has_key()方法
另外补充一点，这个程序用到了numpy，所以需要安装numpy，不知道怎么安装numpy的直接在命令行运行下面指令（前提你的Linux已经安装了python3）
$pip3 install numpy



[1]        Halvor, S, Hansen, et al. A reoptimized GROMOS force field for hexopyranose-based carbohydrates accounting for the relative free energies of ring conformers, anomers, epimers, hydroxymethyl rotamers, and glycosidic linkage conformers[J]. Journal of Computational Chemistry, 2011.
[2]        Huenenberger, Philippe, H, et al. Revision of the GROMOS 56A6(CARBO) force field: Improving the description of ring-conformational equilibria in hexopyranose-based carbohydrates chains[J]. Journal of Computational Chemistry: Organic, Inorganic, Physical, Biological, 2016.

sobereva

gsbear 发表于 2019-11-4 11:14
我的纤维二糖分子模型是用University of Georgia的CCRC中心开发的一个程序GMML生成的，仔细看了应该没有问 ...

你可以发邮件问问作者怎么回事，肯定比我清楚

gsbear

我的纤维二糖分子模型是用University of Georgia的CCRC中心开发的一个程序GMML生成的，仔细看了应该没有问题。
刚才又把53A6Carbo力场作者发布在GROMACS上的一个修正56A6Carbo_R（顺道上传这个立场文件在附件，下载地址http://www.gromacs.org/@api/deki ... OS_56a6_CARBO_R.zip）里的一个1-4葡萄糖三聚体的分子(1-4trimer.pdb)，拿来跑了一下真空状态的EM和NVT控温，结果还是一样的出现了羟基上氢原子和氧原子键旋转过大的问题。照理说作者上传的模型应该是跑过的，不大可能有问题。因此怀疑是不是我NVT控温过程参数有问题？把.mdp文件贴在后面请老师们帮忙分析一下，谢谢。
title                   = OPLS Lysozyme NVT equilibration
define                  = -DPOSRES  ; position restrain the protein
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 100000    ; 1 * 100000 = 100 ps
dt                      = 0.001     ; 1 fs
; Output control
nstxout                 = 500       ; save coordinates every 1.0 ps
nstvout                 = 500       ; save velocities every 1.0 ps
nstenergy               = 500       ; save energies every 1.0 ps
nstlog                  = 500       ; update log file every 1.0 ps
; Bond parameters
continuation            = no        ; first dynamics run
constraint_algorithm    = Lincs     ; holonomic constraints
constraints             = h-bonds   ; bonds involving H are constrained
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                 = system    ; 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 off
pcoupl                  = no        ; no pressure coupling in NVT
; 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

sobereva

gsbear 发表于 2019-11-4 10:13
尝试了将gen_vel设置为no，可以多跑一段时间，但还是同样会出报警，然后又爆掉。
又重新对纤维二糖分子 ...

糖的拓扑信息或参数有问题

gsbear

sobereva 发表于 2019-11-4 08:16
都1fs了就甭用约束了。等跑一阵子再切换到2fs
gen_vel  = yes结合300K是危险的，初速度太大容易崩。把这个 ...

尝试了将gen_vel设置为no，可以多跑一段时间，但还是同样会出报警，然后又爆掉。
又重新对纤维二糖分子不加水，真空下做了优化，然后跑NVT，同样的葡萄糖环上各个羟基中氢氧原子间角度旋转超过30度的问题还是出现，最终爆掉。

请问Sob老师，这样大体可以判断是糖分子结构的问题了是吗？

少年爱吃地三鲜

sobereva 发表于 2019-11-4 09:42
按照控温的算法温度逐渐升上去

谢谢老师！