计算化学公社

标题: gromacs运行pdb2gmx指令出现问题 [打印本页]

作者
Author: momian 时间: 2015-4-24 13:26
标题: gromacs运行pdb2gmx指令出现问题
本帖最后由 momian 于 2015-4-24 13:26 编辑

各位社友，小弟先简单的讲一下我模拟的思路。我研究的方向是酶催化，具体的说就是用离子液体修饰蛋白，以期改善酶学性质。然后希望用gromacs验证一下实验。由于修饰蛋白的离子液体在gromacs中自带的立场中没有，所以我在DS中将离子液体直接接在要修饰的赖氨酸上保存做成非天然氨基酸，将非天然氨基酸直接突变我要修饰的赖氨酸。非天然氨基酸的立场在ATB在线服务器中生成。根据生成的立场文件去修改对应的*.rtp  *.hdb 文件。修改完成之后，我就开始模拟了。输入的指令如下：pdb2gmx -f PPL-modified.pdb -o PPL-modified.gro -p PPL-modified.top -i PPL-modified.itp -water spc -ignh
之后显示错误信息如下：
Now there are 448 residues with 4586 atoms
Making bonds...

WARNING: atom O1 is missing in residue LLS 96 in the pdb file

WARNING: atom O1 is missing in residue LLS 296 in the pdb file    （注：LLS是我保存的非天然氨基酸的名称，96、296、350、351是我选择修饰的位点）

WARNING: atom O1 is missing in residue LLS 350 in the pdb file

WARNING: atom O1 is missing in residue LLS 351 in the pdb file
Fatal error:
There were 4 missing atoms in molecule Protein_chain_A, if you want to use this incomplete topology anyhow, use the option -missing

下面分别贴上*.pdb和*.rtp文件的相关部分
pdb文件如下：ATOM 961  N LLS A  96    49.537  18.584 124.458  1.00 42.53          N
ATOM 962  CA  LLS A  96    50.035  17.479 123.635  1.00 38.85          C
ATOM 963  C LLS A  96    50.855  16.547 124.518  1.00 38.54          C
ATOM 964  O LLS A  96    50.453  15.415 124.782  1.00 42.25          O
ATOM 965  HN  LLS A  96    49.795  19.510 124.259  1.00  0.00          H
ATOM 966  CG  LLS A  96    50.123  18.461 121.294  1.00 20.00          C
ATOM 967  CB  LLS A  96    50.911  18.013 122.496  1.00 38.87          C
ATOM 968  NZ  LLS A  96    50.802  15.273 119.123  1.00 20.00          N1+
ATOM 969  CE  LLS A  96    51.064  16.535 119.941  1.00 20.00          C
ATOM 970  CD  LLS A  96    49.805  17.275 120.398  1.00 20.00          C
ATOM 971  C1  LLS A  96    50.353  17.765 114.702  1.00 20.00          C
ATOM 972  N2  LLS A  96    49.727  17.520 113.493  1.00 20.00          N
ATOM 973  C6  LLS A  96    49.265  18.616 112.602  1.00 20.00          C
ATOM 974  C7  LLS A  96    47.994  18.314 111.751  1.00 20.00          C
ATOM 975  C8  LLS A  96    48.271  17.360 110.570  1.00 20.00          C
ATOM 976  C9  LLS A  96    46.959  17.131 109.798  1.00 20.00          C
ATOM 977  N1  LLS A  96    50.362  16.686 115.428  1.00 20.00          N1+
ATOM 978  C4  LLS A  96    50.825  16.637 116.906  1.00 20.00          C
ATOM 979  C5  LLS A  96    50.672  15.308 117.736  1.00 20.00          C
ATOM 980  O1  LLS A  96    50.466  14.260 117.107  1.00 20.00          O
ATOM 981  C3  LLS A  96    49.411  16.155 113.464  1.00 20.00          C
ATOM 982  C2  LLS A  96    49.794  15.634 114.630  1.00 20.00          C
ATOM 983  HZ  LLS A  96    49.892  14.855 119.500  1.00 20.00          H
ATOM 984  H1  LLS A  96    49.817  18.552 115.233  1.00 24.39          H
ATOM 985  H3  LLS A  96    48.939  15.619 112.641  1.00 24.39          H
ATOM 986  H2  LLS A  96    49.694  14.590 114.927  1.00 24.39          H

rtp文件如下：
[ LLS ]
[ atoms ]
;  name  type  charge  chargegroup
CA CH1    0.286    1
N NT -0.605    1
H    H    0.391    1
C    C    0.493    1
O    O -0.565    1
CG CH2    0.071    2
CB CH2 -0.071    2
NZ    N -0.548    3
  HZ1    H    0.328    3
CE CH2    0.129    3
CD CH2    0.091    3
C1    C -0.118    4
H1 HC    0.273    4
N2 NR    0.221    4
C6 CH2 -0.017    4
C7 CH2    0.130    4
C8 CH2    0.167    4
C9 CH3 -0.060    4
N1 NR    0.221    4
C4 CH2    0.078    4
C5    C    0.693    4
O1    O -0.588    4
C3    C -0.223    5
H3 HC    0.223    5
C2    C -0.223    6
H2 HC    0.223    6

[ bonds ]
;  ai aj
  CA    N    gb_21
  CA    C    gb_27
  CA CB    gb_27
N    H    gb_2
C    O    gb_5
  CG CB    gb_27
  CG CD    gb_27
  NZ HZ1    gb_2
  NZ CE    gb_21
  NZ C5    gb_10
  CE CD    gb_27
  C1 H1    gb_3
  C1 N2    gb_10
  C1 N1    gb_10
  N2 C6    gb_22
  N2 C3    gb_10
  C6 C7    gb_26
  C7 C8    gb_26
  C8 C9    gb_26
  N1 C4    gb_22
  N1 C2    gb_10
  C4 C5    gb_27
  C5 O1    gb_5
  C3 H3    gb_3
  C3 C2    gb_10
  C2 H3    gb_3
C +N    gb_10

[ angles ]
;  ai aj ak
-C    N    H    ga_32
-C    N CA    ga_31
N CA    C    ga_15
N CA CB    ga_15
C CA CB    ga_15
CA    N    H    ga_11
CA    C    O    ga_35
CB CG CD    ga_27
CA CB CG    ga_15
  HZ1 NZ CE    ga_20
  HZ1 NZ C5    ga_18
CE NZ C5    ga_35
NZ CE CD    ga_27
CG CD CE    ga_27
H1 C1 N2    ga_35
H1 C1 N1    ga_35
N2 C1 N1    ga_55
C1 N2 C6    ga_35
C1 N2 C3    ga_57
C6 N2 C3    ga_35
N2 C6 C7    ga_27
C6 C7 C8    ga_27
C7 C8 C9    ga_27
C1 N1 C4    ga_35
C1 N1 C2    ga_35
C3 N1 C2    ga-35
N1 C4 C5    ga_15
NZ C5 C4    ga_19
NZ C5 O1    ga_33
C4 C5 O1    ga_30
N2 C3 H3    ga_25
N2 C3 C2    ga_7
H3 C3 C2    ga_36
N1 C2 C3    ga_7
N1 C2 H2    ga_25
C3 C2 H2    ga_36
CA    C +N    ga_19
O    C +N    ga_33

   小弟已经反复检查了.rtp文件和.pdb文件，确认原子类型和原子数都是匹配的。错误提示“WARNING: atom O1 is missing in residue LLS 96 in the pdb file”中的atom O1 也都有，已经折腾疯了，希望看到的社友能帮忙给一些思路或者建议，谢谢！

作者
Author: sobereva 时间: 2015-4-24 20:31
还是有不对应的。[atoms]里的HZ1和H在pdb里没有。改后再试试

作者
Author: momian 时间: 2015-4-25 08:49

sobereva 发表于 2015-4-24 20:31
还是有不对应的。[atoms]里的HZ1和H在pdb里没有。改后再试试

Sob，你好，谢谢回复。按照提示将【atom】里面的HZ1和H里与rtp不对应的都改过来了，运行相同的指令后，仍然出现相同的报错。（atom O1 缺少……）我试着按照提示在指令中加入：-missing 强行生成*.top文件,打开后查看后O1原子不在*.top文件里，pdb2gmx确实不认*.pdb里面的O1原子。请问这是原因，或者有其他的解决办法吗？
谢谢回复！

作者
Author: momian 时间: 2015-4-25 09:00

sobereva 发表于 2015-4-24 20:31
还是有不对应的。[atoms]里的HZ1和H在pdb里没有。改后再试试

我使用的gromacs版本是gromacs-4.6.1，刚刚还用gromacs-4.5.4试了一下，完全一样的条件报错也是一样的。

作者
Author: kunkun 时间: 2015-4-25 11:55
把hdb也po上来看看？

作者
Author: momian 时间: 2015-4-25 12:42

kunkun 发表于 2015-4-25 11:55
把hdb也po上来看看？

你好，*.hdb文件的相关部分如下：
LLS    5
1    1    H    N    -C    CA
1    1    H2    C2    N1    C3
1    1    H3    C3    C2    N2
1    1    H1    C1    N1    N2
1    1    HZ1    NZ    CE    C5

谢谢回复！

作者
Author: momian 时间: 2015-4-25 19:19

sobereva 发表于 2015-4-24 20:31
还是有不对应的。[atoms]里的HZ1和H在pdb里没有。改后再试试

sob,你好。我把所有O1统一改成其他名称，指令就通过了，现在做位置限制性模拟的时候，gromacs又给了我两个警告。警告具体如下：
WARNING 1 [file pr.mdp]:
  You are generating velocities so I am assuming you are equilibrating a
  system. You are using Parrinello-Rahman pressure coupling, but this can
  be unstable for equilibration. If your system crashes, try equilibrating
  first with Berendsen pressure coupling. If you are not equilibrating the
  system, you can probably ignore this warning.
WARNING 2 [file pr.mdp]:
  You are using pressure coupling with absolute position restraints, this
  will give artifacts. Use the refcoord_scaling option.

我应该听从警告上的建议，修改.mdp文件，还是忽略警告直接运行呢

作者
Author: sobereva 时间: 2015-4-25 20:52
第一个警告无所谓，只要你模拟结果正常不崩溃就行。

第二个最好还是按它说的加上refcoord_scaling，这样的话，盒子尺寸一变，位置限制的绝对位置也会相应地改变。如果确实不想这样那就别管它。

作者
Author: momian 时间: 2015-4-25 21:43

sobereva 发表于 2015-4-25 20:52
第一个警告无所谓，只要你模拟结果正常不崩溃就行。

第二个最好还是按它说的加上refcoord_scaling，这样 ...

sob，谢谢回复，我按提示加上了refcoord_scaling，然后进行pr模拟，结果妥妥的报错了。我看了一下感觉没有头绪，麻烦再给看看啊！
错误信息为：
3 particles communicated to PME node 3 are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension x.
错误提示为：
starting mdrun 'Protein in water'
250000 steps, 500.0 ps.
step 0
Step 1, time 0.002 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 7.300683, max 120.218826 (between atoms 93 and 94)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length

Step 1, time 0.002 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 5.035449, max 150.348923 (between atoms 96 and 98)
bonds that rotated more than 30 degrees:
atom 1 atom 2  angle  previous, current, constraint length
   77    79 40.9 0.1470 0.2080    0.1470
   93    96 51.9 0.1548 9.5750    0.1530
   79    80 41.5 0.1530 0.2143    0.1530
   79    84  114.6 0.1531 0.2782    0.1530
   31    32 94.1 0.1547  13.1804    0.1530
   31    34 79.9 0.1536 1.8931    0.1530
   32    33 81.2 0.1828 9.4893    0.1830
   33    95 98.0 0.2049  17.2489    0.2040
   34    35  105.5 0.1233 1.1312    0.1230
   34    36  103.3 0.1333 1.1029    0.1330
   36    37  105.6 0.1001 0.1639    0.1000
   36    38  104.1 0.1471 0.0869    0.1470
   84    85  121.3 0.1231 0.2714    0.1230
                        .
                        .
                        .
                        .
                  未列完
能量最小化已经完成了，并且没有出任何问题，进行pr的时候就这样了。
pr.mdp文件内容如下：
;
; User spoel (236)
; Wed Nov  3 17:12:44 1993
; Input file
;
title             =  Yo
cpp                =  /usr/bin/cpp
define             =  -DPOSRES
constraints       =  all-bonds
integrator       =  md
dt                =  0.002 ; ps !
nsteps             =  250000 ; total 500 ps.
nstcomm          =  10
nstxout          =  250 ;collect data every 1 ps
nstvout          =  1000
nstfout          =  0
nstlog             =  10
nstenergy          =  10
nstlist          =  10
ns_type          =  grid
rlist             =  1.2
coulombtype       =  PME
rcoulomb          =  1.2
vdwtype          =  cut-off
rvdw             =  1.2
fourierspacing    =  0.12
fourier_nx       =  0
fourier_ny       =  0
fourier_nz       =  0
pme_order          =  4
ewald_rtol       =  1e-5
optimize_fft       =  yes
; V-rescale temperature coupling is on in two groups
Tcoupl             =  V-rescale
tc-grps       =  Protein non-Protein
tau_t             =  0.1 0.1
ref_t             =  300 300
; Energy monitoring
energygrps       =  Protein non-Protein
; Pressure coupling is on
Pcoupl             =  parrinello-rahman
Pcoupltype       =  isotropic
tau_p             =  0.5
compressibility    =  4.5e-5
ref_p             =  1.0
; Generate velocites is on at 300 K.
gen_vel          =  yes
gen_temp          =  300.0
gen_seed          =  173529

麻烦帮忙看一下，到底是哪里出了问题。再次感谢！

作者
Author: sobereva 时间: 2015-4-25 22:01
这种问题肯定是初始结构不好，或者参数不合适。
先把位置限制去掉，步长改成0.001ps试试
这种问题只能通过反复检查、尝试来解决，从最简单的设置开始，步步为营。

欢迎光临计算化学公社 (http://bbs.keinsci.com/)