Gromacs建立非标准残基时修改itp文件的问题

chuxuezhe

各位老师好：
我初学Gromacs，因课题需要，采用Gromacs进行蓝藻蛋白动力学模拟：
主要是参考李老师的文章进行操作（http://blog.sciencenet.cn/home.p ... od=space&uid=548663）
主要步骤：
（1）pdb文件中抽取与非标准残基N端和C端相连的残基（共三个残基）；
（2）采用GV对氢原子补全；
（3）acpype+antechamber建立基于amber力场的拓扑文件；
（4）对产生的以_GMX.itp进行修改，并修改了residuetypes文件，因为没有新原子类型，所以未修改其他力场文件；
问题：
（1）产生的拓扑文件，关于二面角的定义：
;    i      j      k      l      C0         C1         C2         C3         C4         C5     func
     N1    C3    C4    O5    18.20040    0.00000  -18.20040    0.00000    0.00000    0.00000   3   
    ...
但是，所参考的二面角参数仅有2项（图1），使用pdb2gmx时，提示
Fatal error:
Increase MAXSLEN in the grompp code to at least 74, or shorten your definition of bonds like O5    18.20040    0.00000  -18.20040    0.00000    0.00000 0.00000   3 to at most 31
不知道如何处理，有什么可以参考的处理方法？
（2）需要生成hdb文件，文中给出了所有的氢原子的成键方式和书写规则（图2）。但实际的氢原子成键方式特殊，文中并没有介绍，是不是要找结构最相近的代替？
（3）该蛋白质通过二硫键连接了一个小分子，我是要按照蛋白质-配体的方式建拓扑文件（教材中的例子都是未形成共价键的。connect字段表示可以视为连接，是程序按照实际原子距离自动判断么？还是只要定义了相同肽链号，就会连在一起？），还是要重新定义一个新的残基描述这种结构？若需要定义新的残基，是需要把小分子和相连的残基全部定义为一个（据说pdb2gmx判断残基，每个残基都必须是N端和C端相连，否则出错率增加），还是只定义小分子，让程序去自行判断二硫键链接？
问题较多，谢谢各位老师的指导。

sobereva

我没看过那文章。我也不知道你说的“教材中的例子”指什么。如果要按文章作者的做法做，不如直接去问作者。

我不知道你要搞的非标准残基具体是什么样，如果和普通氨基酸差异不大，那么明显不适合用acpype+antechamber，这样搞出来的是基于GAFF力场的参数，而GAFF力场描述氨基酸明显不如Amber力场。基于力场中原本有的各种残基的rtp，举一反三自己写非标准残基的rtp就完了。

自行看力场原文，了解力场中定义的各个原子类型都对应什么情况，gmx自己力场目录下的atomtypes.atp里也有各个原子类型的说明。

“据说pdb2gmx判断残基，每个残基都必须是N端和C端相连，否则出错率增加”这种描述意义不明，应当搞清楚pdb2gmx、rtp文件的基本规则。

氨基酸通过二硫键链接小分子，整体作为一个非标准残基也行；小分子和残基独立定义rtp条目，然后让pdb2gmx自动识别二硫键也行。对于后者，必须自己修改specbond.dat里的判断规则，并且需恰当增加与这个二硫键有关的参数（参考比如G54A7的ffnonbonded.itp的末尾的与二硫键相关的条目）

lonemen

sobereva 发表于 2018-7-4 03:18
我没看过那文章。我也不知道你说的“教材中的例子”指什么。如果要按文章作者的做法做，不如直接去问作者。 ...

嗯嗯，好有用，好详细。谢谢！

chuxuezhe

谢谢Sob老师耐心指导，我明白怎么弄了。非常感谢！

chuxuezhe

老师您好：
查看了相关的资料之后，采用acpype+antechamber建立非标准残基（残基结构与标准残基完全不同，将小分子和与其共价键相连的残基看做一个整体）的rtp文件和hdb文件，非标准参加命名为CYC，结构文件是CYC.zip。通过pdb2gmx建立拓扑结构（采用amber99sb ildn），程序能够成功生成拓扑文件，但是会提示如下警告：
Warning: Long Bond (1234-1291 = 0.413384 nm)
...
Warning: Long Bond (1263-1267 = 0.258767 nm)

查看了全部的警告，发现警告主要出现在以下部分：
（1）主要是新建的非标准残基的键长警告，有些是在bonds中指定成键的，有些是空间相差很远，本就不会成键的原子，
（2）有2个警告是标准残基的（但是去掉了小分子重新建立拓扑文件时，就不会有标准残基的警告），
问题：
（1）在RCSB网站得到的pdb文件中，就发现有些成键的原子初始坐标相差很远，键长很长，请问该如何让处理，
（2）long bond的问题，查不到相关的解决方法，麻烦老师指出错误原因，
（3）原始的pdb文件中有2个残基的缺失，查了整个蛋白质库，都是缺少这两个，不知道如何处理。
附件内容：
蛋白质结构文件，非标准残基结构文件，生成拓扑文件，非标准残基的rtp和hdb文件

谢谢您的指导！

sobereva

chuxuezhe 发表于 2018-7-16 09:54
老师您好：
查看了相关的资料之后，采用acpype+antechamber建立非标准残基（残基结构与标准残基完全不同， ...

最近太忙，没时间仔细看你的文件。如果pdb里本来成键的原子距离就是很长，应该是pdb文件本身精度很差或者错误，应当根据化学直觉恰当把原子放到合适的位置
要么就是弄到acpype处理的时候用的mol2文件连接关系不对，或者自己修改时候弄错了所致

chuxuezhe

谢谢sob老师，我再仔细检查一遍这几个文件。

1971806673

chuxuezhe 发表于 2018-7-5 01:00
谢谢Sob老师耐心指导，我明白怎么弄了。非常感谢！

您好，我模拟的是聚合物，同样构建了非标准残基文件，我在运行 gmx pdb2gmx -f PVDF.pdb -o PVDF.gro 的时候，遇到了和您相似的问题，Fatal error:Increase MAXSLEN in the grompp code to at least 53, or shorten your definition
of bonds like 0.753    2.259    0.000    -3.012    -0.000    0.000 to at most
31，详细的报错我也上传了，请问您当时是怎么解决的呢？可以给我一些解决的思路吗？非常感谢您的回复。

chuxuezhe

1971806673 发表于 2021-8-15 15:57
您好，我模拟的是聚合物，同样构建了非标准残基文件，我在运行 gmx pdb2gmx -f PVDF.pdb -o PVDF.gro 的 ...

你好，我当时构建的文件，其实出现了好几处错误，主要是共价键的识别出错。我觉得需要参考下科学网的非标准残基的构建教程，需要每一处都琢磨明白。严格按照教程就可以成功，但是必须要非常细致。

1971806673

chuxuezhe 发表于 2021-8-16 16:18
你好，我当时构建的文件，其实出现了好几处错误，主要是共价键的识别出错。我觉得需要参考下科学网的非标 ...

非常感谢您的回复，我还有一点疑问，按照教程，刚开始生成的rtp文件中正常二面角有6列参数，而在最后整理之后的rtp却只有2列参数，rtp具有6列参数运行gmx pdb2gmx就会报错Fatal error:
Increase MAXSLEN in the grompp code to at least 74, or shorten your definition of bonds like O5    18.20040    0.00000  -18.20040    0.00000    0.00000 0.00000   3 to at most 31，
我想请教您这个问题是怎么解决的，这个问题已经困扰我好几天了。非常感谢您的帮助，祝好。

chuxuezhe

您好，聚合物我不懂，我就是按照蛋白质非标准残基来说，希望能有点帮助。这个二面角是不是定义非标准残基与标准残基共价键链接的那个。如果是的话，那就得看下定义原子的顺序；相邻的几个分子的参数，也需要一一核对。我当时是就选取了非标准残基和一个与之相连的氨基酸，仅仅对这两个分子做分析，然后排除的问题。此外，您看下初始结构是不是合理（不合理接触或原子距离过大等），看下有没有必要先用半经验方法优化一下。祝你成功。

chenbq18

1971806673 发表于 2021-8-17 09:37
非常感谢您的回复，我还有一点疑问，按照教程，刚开始生成的rtp文件中正常二面角有6列参数，而在最后整理 ...

你好 整理后的rtp文件中的两列数据，是从哪里获得的

0ut0fcontrol

chenbq18 发表于 2021-12-28 16:08
你好 整理后的rtp文件中的两列数据，是从哪里获得的

C0-C5 这6个参数是GROMACS中的Ryckaert-Bellemans dihedral potential，见GROMACS文档 5.5节，p359-361.

RB形式是GROMACS优先支持的格式，不用转化为3个参数periodic形式。

如果一定要从RB形式转为periodic形式，在ParmEd中有讨论，应该能找到转换源码：
https://github.com/ParmEd/ParmEd/issues/714
https://github.com/ParmEd/ParmEd/pull/837/files

苏玖染

这个意思就是说你给的数据太长，可以用一个define的字段来代替，比如：OS_C_CA_CA
然后使用gmx pdb2gmx的时候会在topol.top里面自动用OS_C_CA_CA字段代替你的参数，最后你只用在topol里面#define一下OS_C_CA_CA字段具体的参数值，就行了，注意键参数的函数类型需要删除和修改

具体可以参考TPPMKTOP生成RTP的文件，以及Gromacs自带的力场参数中的RTP文件和ffbonded.itp