计算化学公社

标题: 求助，使用packmol建立了H2O+CTAB的体系以后，下一步如何做才能顺利进行计算 [打印本页]

作者
Author: 牧生 时间: 2020-8-12 22:24
标题: 求助，使用packmol建立了H2O+CTAB的体系以后，下一步如何做才能顺利进行计算
本帖最后由牧生于 2020-8-12 22:37 编辑

最近在做表面活性剂自组装，花了好大力气才装好了带cuda的gromacs。。目前刚开始做，就遇到了难题。

用MS建立并另存了CTAB+，Br-和水分子的pdb文件，然后用packmol粗略的建立了一个模型（先不管是否很合理），命名为ctab+h2o.pdb

tolerance 2.0
add_box_sides 1.2
output ctab+h2o.pdb
structure water.pdb
  number 1000
  inside cube 0. 0. 0. 50.
end structure
structure ctab+.pdb
  number 50
  inside cube 0. 0. 0. 50.
end structure
structure br-.pdb
  number 50
  inside cube 0. 0. 0. 50.
end structure

end structure

自我感觉似乎也没错 (, 下载次数 Times of downloads: 27)

下一步操作：读取pdb文件，生成gro和top文件，使用水模型为tip3p

gmx pdb2gmx -f ctab+h2o.pdb -o ctab+h2o.gro -p ctab+h2o.top -water tip3p

选择力场15 OPLS-AA以后，就报错

Using the Oplsaa force field in directory oplsaa.ff
going to rename oplsaa.ff/aminoacids.r2b
Opening force field file C:\gmx2019.3/share/gromacs/top/oplsaa.ff/aminoacids.r2b
Reading ctab+h2o.pdb...
Read 'Built with Packmol', 6150 atoms
Analyzing pdb file
Splitting chemical chains based on TER records or chain id changing.
There are 3 chains and 0 blocks of water and 1100 residues with 6150 atoms

  chain  #res #atoms
  1 'A'  1000 3000
  2 'B' 50 3100
  3 'C' 50    50

WARNING: there were 3000 atoms with zero occupancy and 0 atoms with
      occupancy unequal to one (out of 6150 atoms). Check your pdb file.       #pdb肯定是有问题的，但是作为新人小白，又不太看得出来哪里有问题，或者如何解决

Opening force field file C:\gmx2019.3/share/gromacs/top/oplsaa.ff/atomtypes.atp
Atomtype 814
Reading residue database... (Oplsaa)
Opening force field file C:\gmx2019.3/share/gromacs/top/oplsaa.ff/aminoacids.rtp
Residue 51
Sorting it all out...
Opening force field file C:\gmx2019.3/share/gromacs/top/oplsaa.ff/aminoacids.hdb
Opening force field file C:\gmx2019.3/share/gromacs/top/oplsaa.ff/aminoacids.n.tdb
Opening force field file C:\gmx2019.3/share/gromacs/top/oplsaa.ff/aminoacids.c.tdb

Back Off! I just backed up ctab+h2o.top to ./#ctab+h2o.top.1#
Processing chain 1 'A' (3000 atoms, 1000 residues)

Warning: No residues in chain starting at 1 identified as Protein/RNA/DNA.
This makes it impossible to link them into a molecule, which could either be
correct or a catastrophic error. Please check your structure, and add all
necessary residue names to residuetypes.dat if this was not correct.

Warning: No residues in chain starting at 2 identified as Protein/RNA/DNA.
This makes it impossible to link them into a molecule, which could either be
correct or a catastrophic error. Please check your structure, and add all
necessary residue names to residuetypes.dat if this was not correct.

Warning: No residues in chain starting at 3 identified as Protein/RNA/DNA.
This makes it impossible to link them into a molecule, which could either be
correct or a catastrophic error. Please check your structure, and add all
necessary residue names to residuetypes.dat if this was not correct.

Warning: No residues in chain starting at 4 identified as Protein/RNA/DNA.
This makes it impossible to link them into a molecule, which could either be
correct or a catastrophic error. Please check your structure, and add all
necessary residue names to residuetypes.dat if this was not correct.

Warning: No residues in chain starting at 5 identified as Protein/RNA/DNA.
This makes it impossible to link them into a molecule, which could either be
correct or a catastrophic error. Please check your structure, and add all
necessary residue names to residuetypes.dat if this was not correct.

Disabling further warnings about unidentified residues at start of chain.
Problem with chain definition, or missing terminal residues.
This chain does not appear to contain a recognized chain molecule.
If this is incorrect, you can edit residuetypes.dat to modify the behavior.
8 out of 8 lines of specbond.dat converted successfully

-------------------------------------------------------
Program:    gmx pdb2gmx, version 2019.3
Source file: src\gromacs\gmxpreprocess\resall.cpp (line 618)

Fatal error:
Residue '' not found in residue topology database                #我如何解决这个问题啊？

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors            #看了手册的说明，这是一个新手常犯的错误，需要残基，但是手册上说的解决办法比较官方，没有实例的解释说明，我理解起来还是比较吃力。
-------------------------------------------------------

得到一个空的top文件

;
;       File 'ctab+h2o.top' was generated
;       By user: Administrator (-1)
;       On host: PCOS-2020NOAPHI
;       At date: Wed Aug 12 22:12:38 2020
;
;       This is a standalone topology file
;
;       Created by:
;                         :-) GROMACS - gmx pdb2gmx, 2019.3 (-:
;
;       Executable: C:\gmx2019.3\bin\gmx.exe
;       Data prefix:  C:\gmx2019.3
;       Working dir:  E:\
;       Command line:
;       gmx pdb2gmx -f ctab+h2o.pdb -o ctab+h2o.gro -p ctab+h2o.top -water tip3p
;       Force field was read from the standard GROMACS share directory.
;

; Include forcefield parameters
#include "charmm27.ff/forcefield.itp"

于是怀疑带电的部分是否有问题，重新建立了只有水的pdb，但是运行gmx pdb2gmx还是一样的报错。

附件为我使用的三个文件
(, 下载次数 Times of downloads: 3)

(, 下载次数 Times of downloads: 6)

(, 下载次数 Times of downloads: 11)

请诸位智者帮帮忙。

作者
Author: sobereva 时间: 2020-8-13 07:47
这种体系本来就不适合用pdb2gmx，除非你把所有分子信息都定义到rtp文件里

这种体系应当用acpype之类程序分别产生各种分子的拓扑信息，然后再合并到一起去
几种生成有机分子GROMACS拓扑文件的工具
http://sobereva.com/266（http://bbs.keinsci.com/thread-428-1-1.html）
Br-在ions.itp里、水在spce.itp里已经有，不需要上述步骤。

作者
Author: 牧生 时间: 2020-8-13 09:45
本帖最后由牧生于 2020-8-13 11:37 编辑

我也试过使用LigParGen得到了CTAB+的pdb，top, itp等文件，如附件中压缩包内的信息，
(, 下载次数 Times of downloads: 22)

但是，下一步，我该怎么继续进行呢？

  "产生.itp文件然后include到主.top文件是最好的作法。"这个操作是如何进行的？？

即使认为得到的pdb，gro，top等生成文件就是正确的，那么，我进行能量最小化，

gmx grompp -f em.mdp -c MOL_95525A-PBC.gro -p MOL_95525A.top -o emctab.tpr

其中，em.pdb内容如下

define       = -DFLEXIBLE
integrator    = steep
emtol          = 500.0
emstep       = 0.01
nsteps       = 1000
nstenergy    = 1
energygrps    = System
nstlist       = 1
ns_type       = grid
coulombtype    = PME
rlist          = 1.0
rcoulomb       = 1.0
vdwtype       = cut-off
rvdw          = 1.0
constraints    = none
pbc          = xyz

得到的结果如下：

  gmx grompp -f em.mdp -c MOL_95525A-PBC.gro -p MOL_95525A.top -o emctab.tpr -maxwarn 2

NOTE 1 [file em.mdp, line 22]:
  em.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 em.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 -1039400983

-------------------------------------------------------
Program:    gmx grompp, version 2019.3
Source file: src\gromacs\gmxpreprocess\grompp.cpp (line 545)

Fatal error:
No molecules were defined in the system    #这个错误，我该如何解决

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------

所谓“然后再合并到一起去”，由于没有实际操作教程流程，我一时间还无法理解。

作者
Author: wyfgg 时间: 2022-5-24 18:44

牧生发表于 2020-8-13 09:45
我也试过使用LigParGen得到了CTAB+的pdb，top, itp等文件，如附件中压缩包内的信息，

楼主老师，我现在就卡在您当年这一步，用packmol建立了mix.pdb，合并了top和itp，运行em.mdp提示指定的原子类型无效

'[ atomtypes ]'
Invalid order for directive atomtypes

复制代码

作者
Author: 牧生 时间: 2022-5-24 20:03

wyfgg 发表于 2022-5-24 18:44
楼主老师，我现在就卡在您当年这一步，用packmol建立了mix.pdb，合并了top和itp，运行em.mdp提示指定的原 ...

http://bbs.keinsci.com/thread-19761-1-1.html
看这个贴子第二楼，第四楼

作者
Author: wyfgg 时间: 2022-5-24 20:35
本帖最后由 wyfgg 于 2022-5-24 20:49 编辑

牧生发表于 2022-5-24 20:03
http://bbs.keinsci.com/thread-19761-1-1.html
看这个贴子第二楼，第四楼

谢谢牧生老师，已经解决了，现在踩到了您当年top和pdb原子类型不匹配这一步，正在研究您的帖子

使用现成的甲烷水合物跑MD，在能量最小化这一步，原子名就对不上，无法继续 - 分子模拟 (Molecular Modeling) - 计算化学公社 - 第2页
http://bbs.keinsci.com/thread-24586-2-1.html

不过我感觉我整体的原子顺序是没问题的，原子类型的顺序对比之后都是对的上的，是不是因为top文件包含原子类型，pdb文件只有元素类型导致的？这样的话可不可以直接 -maxwarn
更新：
看到了Sob老师的帖子6#，直接忽略试试

请问top文件和gro文件的atom name不匹配该怎么办？ - 分子模拟 (Molecular Modeling) - 计算化学公社
http://bbs.keinsci.com/thread-28342-1-1.html

作者
Author: 牧生 时间: 2022-5-24 20:54

wyfgg 发表于 2022-5-24 20:35
谢谢牧生老师，已经解决了，现在踩到了您当年top和pdb原子类型不匹配这一步，正在研究您的帖子
不过我感 ...

只要你确保确实能对应起来，可以忽略。但也有一个问题，warning末尾有提到more than 20 XXXXXXX,表明还有更多的原子不能匹配，有可能实际上是对的，也有可能不对。你先试试忽略这个警告以后跑一阵，再放进vmd看看有没有异常。

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