计算化学公社

标题: 求助：关于如何利用pdb2gmx产生聚合物的top文件的问题（.pdb文件相关） [打印本页]

作者
Author: yihanxu 时间: 2019-2-26 15:14
标题: 求助：关于如何利用pdb2gmx产生聚合物的top文件的问题（.pdb文件相关）
本帖最后由 yihanxu 于 2019-2-26 02:19 编辑

老师好，大家好，我试了一下这个PE的例子： https://mailman-1.sys.kth.se/pip ... 9-March/040125.html ，但是链上有12个C的时候就不行了，请问是我的.pdb格式错了（比如哪列没有对齐）吗？而且用gview产生的.pdb在这种情况下感觉不太好用，比如会有多余的列，要手动处理一下，若是文件长的时候就很麻烦，请问怎么样能直接生成好用的.pdb格式文件呢？
谢谢指教。

作者
Author: sobereva 时间: 2019-2-26 16:27
pdb格式对列有严格的要求，看http://www.wwpdb.org/documentati ... /format33/v3.3.html。某些列可以没有（比如元素），但不能错位。你的pdb明显存在错位。利用ultraedit等程序的列模式进行修改很容易。

作者
Author: yihanxu 时间: 2019-2-27 02:28

sobereva 发表于 2019-2-26 02:27
pdb格式对列有严格的要求，看http://www.wwpdb.org/documentation/file-format-content/format33/v3.3.html ...

谢谢老师。用gview画碳链，碳原子被标记为HETATM，是为什么呢？要手动改成ATOM，感觉麻烦呢，能直接生成高分子pdb的时候就是ATOM吗？

作者
Author: sobereva 时间: 2019-2-27 02:44

yihanxu 发表于 2019-2-27 02:28
谢谢老师。用gview画碳链，碳原子被标记为HETATM，是为什么呢？要手动改成ATOM，感觉麻烦呢，能直接生成 ...

没什么麻烦的，用ultraedit列模式五秒钟就搞定了

作者
Author: yihanxu 时间: 2019-2-27 02:58

sobereva 发表于 2019-2-26 12:44
没什么麻烦的，用ultraedit列模式五秒钟就搞定了

谢谢老师

作者
Author: yihanxu 时间: 2019-3-1 09:37
本帖最后由 yihanxu 于 2019-2-28 20:48 编辑

老师好，用gmx pdb2gmx命令后，会让选择力场，用的GAFF并没被列出来，请问这时候应该怎么办呢？选择 6: AMBER99SB-ILDN protein, nucleic AMBER94是好方法吗？谢谢。

作者
Author: sobereva 时间: 2019-3-1 11:55

yihanxu 发表于 2019-3-1 09:37
老师好，用gmx pdb2gmx命令后，会让选择力场，用的GAFF并没被列出来，请问这时候应该怎么办呢？选择 6: AMB ...

本来就不会被列出来，只有把xxx.ff放到了top目录下才会在选择力场时出现xxx。得把机制搞清楚
当前用GAFF的参数，和这里选什么毫无关系
如果你基于GAFF自己写了rtp，应当新建个力场目录，比如叫GAFF.ff

作者
Author: yihanxu 时间: 2019-3-1 13:32
本帖最后由 yihanxu 于 2019-2-28 23:57 编辑

sobereva 发表于 2019-2-28 21:55
本来就不会被列出来，只有把xxx.ff放到了top目录下才会在选择力场时出现xxx。得把机制搞清楚
当前用GAFF ...

谢谢老师。我把gaff.ff文件夹放在top目录里了，gaff.ff文件夹里有两个文件：gaff.hdb和gaff.rtp。但是使用pdb2gmx命令时备选力场里没有gaff，是怎么回事呢？

作者
Author: yihanxu 时间: 2019-3-1 14:42
老师好，GROMACS官网给的聚乙烯的例子中.rtp文件中没有[bondedtypes]字段，但是看GROMACS手册说这是强制的，这是怎么回事呢？哪里可以查到下面各项的定义呢？谢谢。
[ bondedtypes ]
; bonds angles dihedrals impropers all_dihedrals nrexcl HH14 RemoveDih
1 1 9 4 1 3 1 0

作者
Author: sobereva 时间: 2019-3-1 18:17

yihanxu 发表于 2019-3-1 13:32
谢谢老师。我把gaff.ff文件夹放在top目录里了，gaff.ff文件夹里有两个文件：gaff.hdb和gaff.rtp。但是使 ...

照着其它.ff里的内容去写。明显你还没设置forcefield.doc

作者
Author: sobereva 时间: 2019-3-1 18:18

yihanxu 发表于 2019-3-1 14:42
老师好，GROMACS官网给的聚乙烯的例子中.rtp文件中没有字段，但是看GROMACS手册说这是强制的，这是怎么回事 ...

如果rtp里直接定义了参数，则bondedtypes虽然有，但内容可以留空，反正也不需要了。直接定义了参数的情况程序就不会再自动从bondedtypes里试图找匹配的参数

作者
Author: yihanxu 时间: 2019-3-2 15:40
老师好，请问蛋白质和高分子为什么不能用antechamber生成拓扑文件，而要用pdb2gmx命令?
谢谢。

作者
Author: sobereva 时间: 2019-3-2 15:49

yihanxu 发表于 2019-3-2 15:40
老师好，请问蛋白质和高分子为什么不能用antechamber生成拓扑文件，而要用pdb2gmx命令?
谢谢。

本来antechamber就是搞小分子的
巨大的分子本身超过了antechamber程序设计角度能处理的极限，而且其调用SQM进行半经验优化、计算AM1电荷的步骤，对于几千原子体系也根本算不动。
更何况antechamber产生的拓扑文件是基于GAFF力场的，对于蛋白质有明显更适合的蛋白质专用力场。

作者
Author: yihanxu 时间: 2019-3-4 12:58
本帖最后由 yihanxu 于 2019-3-5 11:49 编辑

老师好，生成的高分子的.top文件中没有bond、angle、dihedral等成键参数。为了获得GAFF下的这些参数，我是用antechamber+acpype.py生成了低聚物的.itp文件，然后手动修改低聚物的.itp文件得到的高聚物的成键参数。但是这样耗时有点多，尤其是把低聚物.itp文件中原子序号替换成原子类型时很累眼睛，请问有更好、更快捷的办法吗？谢谢。

作者
Author: tjuptz 时间: 2019-9-9 22:38
请问楼主后来找到更好的方法了吗？我最近也在编rtp，但我的重复单元里分子数有点多，好像从acpype生成的itp改rtp很不容易……

作者
Author: yihanxu 时间: 2019-9-10 05:35
本帖最后由 yihanxu 于 2019-9-9 15:37 编辑

tjuptz 发表于 2019-9-9 08:38
请问楼主后来找到更好的方法了吗？我最近也在编rtp，但我的重复单元里分子数有点多，好像从acpype生成的itp ...

还是用的这种方法。如果重复单元的基团比较多，确实改起来很麻烦，我的重复单元没超过三个环，改起来还行。

作者
Author: tjuptz 时间: 2019-9-10 12:07
本帖最后由 tjuptz 于 2019-9-10 12:54 编辑

sobereva 发表于 2019-3-1 11:55
本来就不会被列出来，只有把xxx.ff放到了top目录下才会在选择力场时出现xxx。得把机制搞清楚
当前用GAFF ...

老师，在GROMACS中，基于acpype获得高分子单元的GAFF参数，并对聚合物用pdb2gmx生成链的top文件，必须新建GAFF.ff而不能基于amber99sb-ildn.ff做吗？因为后面我也想用隐式溶剂。另外，编rtp时需要把itp里的参数都写进去吗？求指点

作者
Author: sobereva 时间: 2019-9-11 09:52

tjuptz 发表于 2019-9-10 12:07
老师，在GROMACS中，基于acpype获得高分子单元的GAFF参数，并对聚合物用pdb2gmx生成链的top文件，必须新 ...

可以用AMBER的.ff。和GAFF的原子类型并不冲突，力场的default部分也都是兼容的。只要把acpype给出的额外参数补进去就行了

rtp里不是非得直接体现出参数。只要把成键关系在里面定义了就够了，grompp的时候自动就会从力场文件里根据原子类型去找对应的参数。

作者
Author: tjuptz 时间: 2019-9-11 13:52
本帖最后由 tjuptz 于 2019-9-11 19:37 编辑

sobereva 发表于 2019-9-11 09:52
可以用AMBER的.ff。和GAFF的原子类型并不冲突，力场的default部分也都是兼容的。只要把acpype给出的额外 ...

我明白了，谢谢老师。

作者
Author: tjuptz 时间: 2019-9-13 15:56
本帖最后由 tjuptz 于 2019-9-15 10:08 编辑

sobereva 发表于 2019-9-11 09:52
可以用AMBER的.ff。和GAFF的原子类型并不冲突，力场的default部分也都是兼容的。只要把acpype给出的额外 ...

再次求助社长。我把低聚体acpype给出的参数补到力场文件里，其中二面角字段的原子名改成原子类型后遇到如下问题：
os  c3  c3  hc 9    0.00 0.00000 0 ;
os  c3  c3  hc 9    0.00 1.04600 1 ;
c3  c3  c3  c3 9    0.00 0.75312 3 ;
c3  c3  c3  c3 9    180.00 0.83680 1 ;
c3  c3  c3  c3 9    180.00 1.04600 2 ;
我查看itp中对应原始原子名字段是如下这样的：
17    18    19    51    9    0.00 0.00000 0 ; O17- C18- C19- H51
17    18    19    51    9    0.00 1.04600 1 ; O17- C18- C19- H51
17    18    19    52    9    0.00 0.00000 0 ; O17- C18- C19- H52
17    18    19    52    9    0.00 1.04600 1 ; O17- C18- C19- H52

30    31    32    76    9    0.00 0.00000 0 ; O30- C31- C32- H76
30    31    32    76    9    0.00 1.04600 1 ; O30- C31- C32- H76
30    31    32    77    9    0.00 0.00000 0 ; O30- C31- C32- H77
30    31    32    77    9    0.00 1.04600 1 ; O30- C31- C32- H77

18    19    20    21    9    0.00 0.75312 3 ; C18- C19- C20- C21
18    19    20    21    9 180.00 0.83680 1 ; C18- C19- C20- C21
18    19    20    21    9 180.00 1.04600 2 ; C18- C19- C20- C21

19    20    21    22    9    0.00 0.75312 3 ; C19- C20- C21- C22
19    20    21    22    9 180.00 0.83680 1 ; C19- C20- C21- C22
19    20    21    22    9 180.00 1.04600 2 ; C19- C20- C21- C22

20    21    22    23    9    0.00 0.75312 3 ; C20- C21- C22- C23
20    21    22    23    9 180.00 0.83680 1 ; C20- C21- C22- C23
20    21    22    23    9 180.00 1.04600 2 ; C20- C21- C22- C23

…………
26    27    28    29    9    0.00 0.75312 3 ; C26- C27- C28- C29
26    27    28    29    9 180.00 0.83680 1 ; C26- C27- C28- C29
26    27    28    29    9 180.00 1.04600 2 ; C26- C27- C28- C29
我的单体里面有长碳链C18-C29，通过acpype获得碳链上的C都是c3原子类型，但这些c3碳原子二面角acpype却给出了多项参数，我想起来，这是多项叠加描述扭转势。这样补到ffbonded.itp中就出现了同样的c3-c3-c3-c3后跟着不同参数的情况；同样的，c3碳原子上的H都是hc原子类型，也遇到了这个问题。我看其他帖子中提到，力常数为零的可以忽略，但对于c3-c3-c3-c3还有三项参数。这时候我该怎么处理呢？把这些参数直接写到rtp里吗，我尝试了一下，参数太多不太现实；而且根据单体的itp参数补还得考虑多聚体成键项改变。

更新：我看了手册，说是二面角函数类型9就可以用多参数，所以我就这样补进ffbonded.itp就可以了吧？

作者
Author: diaok 时间: 2019-9-14 12:29

tjuptz 发表于 2019-9-13 15:56
再次求助社长。我把acpype给出的二面角字段的原子名改成原子类型后遇到如下问题：
os c3 c3 hc 9 ...

写rtp时有没有涉及到片段连接时应该补全的键角二面角之类？
要是不好解决可以直接做成分开的片段，然后用parmed合并成一个拓扑分子，再根据原子序号补上连接的键角之类参数

作者
Author: tjuptz 时间: 2019-9-14 12:39

diaok 发表于 2019-9-14 12:29
写rtp时有没有涉及到片段连接时应该补全的键角二面角之类？
要是不好解决可以直接做成分开的片段，然后 ...

涉及到了，我用共聚体重复单元的itp应该能获得这些参数吧

作者
Author: diaok 时间: 2019-9-14 16:02

tjuptz 发表于 2019-9-14 12:39
涉及到了，我用共聚体重复单元的itp应该能获得这些参数吧

gmx给的那个rtp示例里也没看到键角二面角之类信息，应该是用ffbond里面给补全了
实现gaff应该需要把有关参数都补充完整
It is easier to create .rtp file first, since you only need partial charge and bond information. The dihedral and angle will be automatically generated. Then you run grompp to check for any non-standard bond/angle/dihedral, correct for them (modify the ffbonded.itp).
https://www.researchgate.net/pos ... mulation_in_Gromacs

我一般用的另一种方法，先做好单独的片段，然后直接整合pdb连成完整的聚合物，用antechamber来判断拓扑结构得到各种键连参数，结构需要稍微优化下让自动判断不出错

作者
Author: tjuptz 时间: 2019-9-14 16:54

diaok 发表于 2019-9-14 16:02
gmx给的那个rtp示例里也没看到键角二面角之类信息，应该是用ffbond里面给补全了
实现gaff应该需要把有关 ...

对于很长的聚合物，不能直接用antechamber吧

作者
Author: diaok 时间: 2019-9-15 00:10

tjuptz 发表于 2019-9-14 16:54
对于很长的聚合物，不能直接用antechamber吧

antechamber只限制做单个分子，长了也不要紧，起码我做的800原子可以，只有些说键太多的warning，另外check检查缺失参数部分可以直接拿片段的用，不然要检查特别久。

作者
Author: tjuptz 时间: 2019-9-15 08:10
本帖最后由 tjuptz 于 2019-9-15 10:01 编辑

diaok 发表于 2019-9-15 00:10
antechamber只限制做单个分子，长了也不要紧，起码我做的800原子可以，只有些说键太多的warning，另外che ...

用于800个原子，antechamber做AM1单点不知道我的机子做不做得了；我今天试了下我写的rtp，结果如下：
Now there are 4 residues with 194 atoms
Making bonds...
Number of bonds was 200, now 197
Generating angles, dihedrals and pairs...
Before cleaning: 514 pairs
Before cleaning: 514 dihedrals
Keeping all generated dihedrals
Making cmap torsions...
There are  514 dihedrals, 20 impropers,  364 angles
         502 pairs,    197 bonds and    0 virtual sites
Total mass 1143.680 a.m.u.
Total charge -0.590 e
Writing topology

Writing coordinate file...
            --------- PLEASE NOTE ------------
You have successfully generated a topology from: biru3.pdb.
The Gaff force field is used.

请问Number of bonds was 200, now 197 是提示出有什么问题吗？

作者
Author: diaok 时间: 2019-9-15 10:57

tjuptz 发表于 2019-9-15 08:10
用于800个原子，antechamber做AM1单点不知道我的机子做不做得了；我今天试了下我写的rtp，结果如下：
No ...

antechamber这时候只用来生成拓扑，电荷用片段已经算过了，键的数目要对比好，不匹配肯定不行

作者
Author: tjuptz 时间: 2019-9-30 20:19

yihanxu 发表于 2019-2-27 02:28
谢谢老师。用gview画碳链，碳原子被标记为HETATM，是为什么呢？要手动改成ATOM，感觉麻烦呢，能直接生成 ...

感觉用MS建聚合物的模型更好一些，格式更规范，而且另存时原子顺序不会乱

作者
Author: tjuptz 时间: 2019-11-30 23:01

sobereva 发表于 2019-3-1 11:55
本来就不会被列出来，只有把xxx.ff放到了top目录下才会在选择力场时出现xxx。得把机制搞清楚
当前用GAFF ...

再次借楼请教老师，比如我有很多重复单元要补到ffbond文件中，那对于从不同重复单元获得的同样成键类型（比如 ca h）需要去重吗？对于多次定义的dihedral类型9会叠加，那bond angle呢？是用最后定义的吗？

作者
Author: sobereva 时间: 2019-12-2 02:53

tjuptz 发表于 2019-11-30 23:01
再次借楼请教老师，比如我有很多重复单元要补到ffbond文件中，那对于从不同重复单元获得的同样成键类型（ ...

会用最后一次定义的

作者
Author: tjuptz 时间: 2019-12-2 06:23

sobereva 发表于 2019-12-2 02:53
会用最后一次定义的

谢谢老师

作者
Author: 哪有这么脆弱 时间: 2021-6-10 00:07

sobereva 发表于 2019-3-1 18:17
照着其它.ff里的内容去写。明显你还没设置forcefield.doc

请问自己建立的gaff.ff中ffbonded.itp等文件的参数是copy rtp文件中的参数数值嘛？

作者
Author: sobereva 时间: 2021-6-15 18:39

哪有这么脆弱发表于 2021-6-10 00:07
请问自己建立的gaff.ff中ffbonded.itp等文件的参数是copy rtp文件中的参数数值嘛？

ambertools程序包里有个gaff.dat文件，里面有所有GAFF的力场参数，你可以试图转成gmx的格式

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