本帖最后由 JohnCase 于 2025-12-18 16:38 编辑
此教程用pymol编辑分子，multiwfn计算resp电荷，用sobtop产生.rtp文件，以及手写用于加氢的.hdb文件，可以让gmx pdb2gmx直接生成（稍加修改即可使用的）topol.top文件。
此流程摸索了一天，然后在supernova的指导下半天通关，https://zhuanlan.zhihu.com/p/87402839
确实繁琐，但是走完此流程会加深对gromacs力场文件的格式的理解。
这里以绿荧光蛋白为例，PBDID 5exb，64号发光核残基5SQ在表达时由三个氨基酸缩合而成，前一个残基是61号，后一个是65号。
一、准备mol2文件
第一件要操作的事就是改原子名，加氢，改氢原子名
原子名对于gromacs识别肽键的键连顺序至关重要，否则gromacs不知道前一个氨基酸的C要连接本残基哪一个N。 根据Amber力场的标准原子命名惯例，设涉及肽键的CNOH四个原子就是分别命名为‘C’, ‘N’, ‘H’, ‘O’
我们点击选择(pymol左上角Residues选择模式改成Atoms)这个5SQ残基的N端，以及C端的C和O，终端输入 label sele, name
发现他们的名不是标准命名， 终端输入（直接复制粘贴，敲击回车） - select C_tmp, sele and element C
- alter C_tmp, name='C'
- select N_tmp, sele and element N
- alter N_tmp, name='N'
- select O_tmp, sele and element O
- alter O_tmp, name='O'
- label C_tmp or N_tmp or O_tmp, name
此时显示修改成标准原子名了 save 5exb_monomer.pdb, all 保存修改原子名后的pdb文件 （顺便打开文本编辑器检查一下存下来的5exb_monomer.pdb中5SQ残基前后是否有TER，有的话删掉，确保5SQ是完整融入了一条链。）
同时，Residues模式选择整个5SQ残基， save 5SQ.pdb, sele 保存修改原子名后的5SQ残基
在新窗口打开5SQ.pdb，点击builder，画甲氨基和乙酰基，对其封端加氢，
接下来依次点击圈出来的这几个，修改相应的原子名，最后label all, name刷新原子名的显示。 这一步是为了避免麻烦，N端的H必须叫H，其它4个H需要修改是因为pdb2gmx -ignh的时候加上的氢的原子名有规律限制。(详见hdb文件一节) - alter sele, name=’H’
- alter sele, name=’HA1’
- alter sele, name=’HA2’
- alter sele, name=’HB1’
- alter sele, name=’HB2’
检查原子名是否有冲突。没有问题的话保存这个残基 - save 5SQ_capped.mol2, all
- save 5SQ_capped.xyz, all
xyz格式喂给xtb（计算hessian，xtbopt.xyz喂给RESP_ORCA_noopt做RESP电荷计算）。 mol2格式喂给sobtop产生rtp文件。
用文本编辑器打开5SQ_capped.mol2（或者保存为pdb看着更顺眼，我这里额外save 5SQ_capped.pdb, all）确认25个重原子的顺序和原始pdb文件的5SQ残基的原子名是一致的， 至于后加的H和封端基团暂时不管，后续在rtp文件中删除即可。
二、手写hdb文件
这里引用supernova的贴子， https://zhuanlan.zhihu.com/p/87402839
hdb文件其实就是加氢规则，用于描述在什么原子上加什么多少氢，以及定义氢原子名，这6种氢够用了
标准氨基酸的hdb文件可以参考 gromacs2023/share/gromacs/top/amber14sb_parmbsc1.ff/aminoacids.hdb
比如GLN残基的加氢规则 - GLN 5
- 1 1 H N -C CA
- 1 5 HA CA N CB C
- 2 6 HB CB CA CG
- 2 6 HG CG CB CD
- 2 3 HE2 NE2 CD CG
其中第一行为残基名称和行数，接下来每行列出一类待添加氢的信息： 加氢数目，加氢方式，加氢名称，原子i，原子j，原子k，原子l。 其中“-C”表示前一个氨基酸的原子名为C的原子，也就是羧基C原子。
那么5SQ的N端H 可以写成 1 1 H N -C CA1 不写成1 1 H N CO2 CA1是因为前CO2是后加的cap，要换成一个氨基酸的位置，我标记为灰色。并且注意是CA1而不是CA，因为这个非标准氨基酸根本没标准意义上的alpha 碳。 并且gromacs不通过CA的识别去搭建拓扑，而是识别肽键的CN原子名。 对于某个原子上加了多个氢，比如type6，亚甲基氢，因为要加两个H，命名规则是在指定的H原子名后面写1和2， 比如2 6 HA CA3 N3 C ，加上的氢名为HA1和HA2。
5SQ.hdb完整内容如下（13表示13行） - 5SQ 13
- 1 1 H N -C CA1
- 1 5 H06 CA1 N CB1 C1
- 2 6 HA CA3 N3 C
- 2 6 HB CB1 CA1 CG1
- 1 1 H21 C2H N2H CG1
- 1 1 H20 N1H CG1 C1H
- 1 1 H19 C1H N1H N2H
- 1 1 H01 CB2 CA2 CG2
- 1 1 H14 CD1 CG2 CE1
- 1 1 H16 CE1 CD1 CZ1
- 1 2 H18 OH CZ1 CE1
- 1 1 H17 CE2 CZ1 CD2
- 1 1 H15 CD2 CG2 CE2
三、获得rtp文件
包括计算hessian矩阵（我这里图省事用的xtb）
xtb *xyz --ohess --gbsa h2o -P 6
用RESP_ORCA_noopt脚本计算RESP电荷，但是要修改Multiwfn_noGUI SP部分的内容 - ### Calculate RESP charges
- echo Running Multiwfn_noGUI...
- Multiwfn_noGUI SP.molden << EOF
- 7
- 18
- 6
- 1
- chgcons.txt
- 2
- y
- 0
- 0
- q
- EOF
- # ### Finalize
- # chgname=${1//$suffix/chg}
- # mv SP.chg $chgname
- # # rm -f SP.* SP_* Nval.txt
- bash RESP_ORCA_noopt.sh xtbopt.xyz 0 1 gas
其中chgcons.txt文件用于限制封端基团都是电荷为0，那么剩下的非标准残基的电荷也自然为0，其内容为两行原子序号，后面的0表示限制这些原子组的总电荷为0。 - 2,3,25,26,27,28 0
- 1,6,7,8,9,23 0
原子序号在pymol里面看就行，label all, id, 或者写个python脚本打印选中的原子，比如名为id.py文件 - from pymol import cmd
- def printid(selection):
- """Add hydrogen to selected sp2 carbon atoms with two bonds, assuming 120° bond angles."""
- atoms = cmd.get_model(selection).atom
- print(f'There are {len(atoms)} atoms in the selection')
- ids = [atom.id for atom in atoms]
- names = [atom.name for atom in atoms]
- print(*ids, sep=',')
- print(', '.join(f'"{name}"' for name in names))
在pymol命令行(确保id.py在pymol的当前目录，控制台pwd/ls可以查看) 然后printid sele，就可以打印输出当前选择的原子的id和原子名。
原子id复制粘贴chgcons.txt文件 （这里原子名的用处是在接下来的rtp文件中删除封端基团原子参数）
下一步根据RESP电荷计算结束后生成的原子电荷文件SP.chg，用sobtop创建rtp文件
cd sobtop_1.0dev5目录， - sobtop < EOF
- /root/jojo/covalent/complete/5SQ_capped.mol2
- 7
- 10
- /root/jojo/covalent/resp/SP.chg //读取RESP电荷文件
- 0
- 3 // Generate GROMACS .rtp file
- 4 // Assign AMBER atom types as much as possibl
- 0
- 3 // 先用内置的成键参数，缺的用hessian 计算
- /root/jojo/covalent/complete/xtb/hessian
- /root/jojo/covalent/complete/gmx/01_protein_top/origin.rtp
- 0
- EOF
这里全部指认AMBER原子类型。
得到的origin.rtp文件需要以下5处修改，我后面用一个python脚本完成
1.前两部分都注释掉，从[ bondedtypes ]部分开始保留。 2.把残基名[ MOL ]改成[ 5SQ ] 3.在[ bonds ]部分最后加上 表示前一个氨基酸的羧基碳原子C连接此残基的N 4.在[ impropers]部分最后加上 - -C CA1 N H 4 180.0 4.60240 2 ; C-CX-N-H from ffbonded.itp
- CA3 +N C O 4 180.0 43.93200 2 ; X-X-C-O from ffbonded.itp
非标准残基其实没有alphaC，但是我这里用的alphaC的参数，目的是保持肽键平面。 tips: 如果不加180.0 4.60240 2 180.0 43.93200 2这两行参数，那么pdb2gmx会从ffbonded.itp文件里面找到C CT N H的参数（恰好跟C CX N H是一样的），因为CA1被sobtop指认为CT原子类型了。 CA3 +N C O倒无所谓，总是用ffbonded.itp文件里X X C O的参数（X指代任意原子）
5.删除任何封端基团参与的参数 因为甲氨基和乙酰基不会出现在pdb文件里，跟他们有关的参数必须全部删除。
以上操作我合并在一个rtp.py脚本里了
rtp.py
(3.12 KB, 下载次数 Times of downloads: 7)
执行之 生成的5SQ.rtp 文件cp到gromacs2023/share/gromacs/top/amber14sb_parmbsc1.ff/5SQ.rtp 并且在/root/jojo/software/gromacs2023/share/gromacs/top/residuetypes.dat文件第一行插入 表示5SQ是个有参数的残基
三、创建topol.top 终于到正题了： - gmx pdb2gmx -f 5exb_monomer.pdb -o 5exb_monomer.gro -ignh -water spc -nobackup
GROMACS reminds you: "I'm a Wishbone and I'm Breaking" (Pixies) 此时成功创建拓扑文件topol.top，但是还有点小毛病，[ dihedrals ] 部分多写了一个9或者4，这里用sed命令删除之，注意数字前后空格 - sed -i 's/ 9 9 / 9 9 /g' topol.top
- sed -i 's/ 9 2 / 2 /g' topol.top
- sed -i 's/ 4 4 / 4 /g' topol.top
然后就可以正常使用了。
最后基于这份topol.top文件可以进入正常的gromacs建模环节，比如加水、加离子，能量最小化，NPT，一步到底没问题 检查一下跑过NPT的两边肽键，完美:-)
| 
发表于 Post on 2025-12-8 20:43:02
|
只看该作者 Only view this author
|只看大图 View big picture
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
楼主 Author
