gromacs运行pdb2gmx指令出现问题

momian

各位社友，小弟先简单的讲一下我模拟的思路。我研究的方向是酶催化，具体的说就是用离子液体修饰蛋白，以期改善酶学性质。然后希望用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 也都有，已经折腾疯了，希望看到的社友能帮忙给一些思路或者建议，谢谢！

sobereva

还是有不对应的。[atoms]里的HZ1和H在pdb里没有。改后再试试

momian

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

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

momian

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

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

kunkun

把hdb也po上来看看？

momian

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

谢谢回复！

momian

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文件，还是忽略警告直接运行呢

sobereva

第一个警告无所谓，只要你模拟结果正常不崩溃就行。

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

momian

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

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

sobereva

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