参考1:GROMACS非标准残基教程1：修改力场与增加残基
参考2:GROMACS非标准残基教程2：芋螺毒素小肽实例
参考3使用AmberTools+ACPYPE+Gaussian创建小分子GAFF力场的拓扑文件
前言
先看上面的参考资料，下面的介绍仅为参考资料的补充。
pdb2gmx在生成拓扑结构时按照以下顺序：r2b->hdb->rtp->tdb
aminoacids.c.tdb和aminoacids.n.tdb是对端基的处理，在amber力场中均为空。（大概因此amber力场无-ter）
下面以amber99sb-ildn力场为例介绍端基残基的生成：
实现自定义氨基酸残基要处理的文件主要有4个：
- top/residuetypes.dat
- top/amber99sb-ildn_m.ff/aminoacids.r2b
- top/amber99sb-ildn_m.ff/xxx.hdb
-
top/amber99sb-ildn_m.ff/xxx.rtp
- xxx为自定义残基的名字，为了不污染原有的aminoacids.hdb和aminoacids.rtp而将其独立出来。
- 如果不编写residuetypes.dat文件，则会提示“The residues in the chain DCF1–DA53 do not have a consistent type. The first residue has type ‘Other’, while residue DA3 is of type ‘DNA’. Either there is a mistake in your chain, or it includes nonstandard residue names that have not yet been added to the residuetypes.dat file in the GROMACS library directory. If there are other molecules such as ligands, they should not have the same chain ID as the adjacent protein chain since it’s a separate molecule.”
步骤
生成top文件
- 片段结构的获取：若为端基残基，则用另一中性氨基酸配成酰胺键后用H将外加氨基酸的C端加氢（或N端减氢），实现中性化封端。若为中间氨基酸，则两端用其它氨基酸分别配成酰胺键后用H将外加氨基酸C端加氢（或N端减氢），实现中性化封端。（若为DNA则需注意磷酸带1个负电，否则电子数和自旋不匹配）
- 生成mol2文件，现发现VMD生成的mol2文件不能用，gaussian view生成的可以。
- 利用acpype生成含有自定义残基的片段结构的拓扑文件，生成的同时确定好带电情况，因为两端均做中性化，所以带电量为自定义残基带电量。检查生成的gro文件结构合理性。acpype要用sf版，sf版中二面角默认类型为9或4，可通过-z改变；github版中为3或1。（见参考3）
xxx.rtp文件的生成
- 拓扑文件中只保留下列部分的内容后另存为rtp格式: [ atoms ], [ bonds ], [ angles ], [ dihedrals ] ; propers, [ dihedrals ] ; impropers。
- 调整电荷。如果采用AM1-BCC电荷, 简单的处理方法是将相邻残基的净电荷加到相应的连接原子上。RESP电荷可采用约束拟合，看参考资料如RESP拟合静电势电荷的原理以及在Multiwfn中的计算。
- acpype生成的拓扑文件的原子类型均为小写，将其改为大写，其中C3改写为CT，HN改写为H。
- 效仿自带库中的氨基酸修改原子名称，可参考氨基酸在PDB文件中的原子命名规则
- rtp文件主要内容如下，其中XXX为自定义残基名，处理方式见参考资料，文件内容样式可参考力场自带的aminoacids.rtp但内容比其要多。
[ bondedtypes ]
;
1 1 9 4 1 3 1 0
[ XXX ]
[ atoms ]
...
[ bonds ]
...
[ angles ]
...
[ dihedrals ] ; propers
...
[ impropers ]
...
xxx.hdb文件的生成
参照参考资料和力场自带的aminoacids.hdb的内容自行编写，原子名称一定要与前面写的xxx.rtp文件一致。注意：如果一个原子上要加多个H时会按照xxx.hdb里的H原子名+序号的顺序生成，比如要加3个名为“HA”的氢原子，则会在位点生成名为“HA1”“HA2”“HA3”的3个H原子，如果这3个原子名称在xxx.rtp文件中没有定义则会报错，所以一定要在xxx.rtp中做好相应的定义。
生成文件的放置位置
将生成的xxx.rtp和xxx.hdb添加到自定义力场top/amber99sb-ildn_m.ff/中，在top/residuetypes.dat中添加自定义残基的名称。
检查
找一个包含自定义残基的pdb文件并将其中的自定义残基名与xxx.rtp调成一致。
用 gmx pdb2gmx -f xxx.pdb -ignh 检查能否正常运行。
Q&A
- 之前试过用gaussian view画出来的pdb结构跑不通（会包其他残基的错误）但从网上下载的pdb可以跑通，还不清楚原因。
- Warning:”Residue 1 named MET of a molecule in the input file was mapped to an entry in the topology database, but the atom H used in an interaction of type angle in that entry is not found in the input file. Perhaps your atom and/or residue naming needs to be fixed.”是个废话警告，参见：链接1或链接2。此情况出现，可能因为有多个同种H虽在[atom]里做了区分但在后面的[bond][angle]等里未区分，将其区别表示后即可解决，不解决也可能得到正确结果，因为参数一样。
- WARNING: “Duplicate line found in or between hackblock and rtp entries.” 说明rtp有问题，比如没有设定键长和力常数。参考链接
小工具
- 残基原子取代修饰工具DIYtool_residue_modification
- itp文件转rtp文件工具DIYtool_top2rtp
附录
自己写过的几个自定义残基：
中性的C端LEU:clec.rtp
中性的C端LEU:clec.hdb
中性的N端VAL:nvan.rtp
中性的N端VAL:nvan.hdb
中性的C端LYN:clyn.rtp
中性的C端LYN:clyn.hdb
中性的N端PHE:nphn.rtp
中性的N端PHE:nphn.hdb
手机版“神探玺洛克”请扫码
