捐赠本论坛 (Donate this site)
开启辅助访问 切换到宽版

计算化学公社

 找回密码 Forget password
 注册 Register
在11月21-24日期间，本论坛仅限等级≥level 4的成员发贴，带来的不便请谅解！其他成员若有学术问题请等到能发贴的时候再发帖。此期间不要给管理员sobereva发私信问学术问题
计算化学公社»论坛首页 › 理论与计算化学 (Theoretical and Computational Chemistry) › 分子模拟 (Molecular Modeling) › 使用MCPB.py生成二茂铁中Fe的力场参数
返回列表 Return 发新帖
Views: 6226|回复 Reply: 5
打印 Print 上一主题 Last thread 下一主题 Next thread

[综合交流] 使用MCPB.py生成二茂铁中Fe的力场参数

[复制链接 Copy URL]
大村驴

27

帖子

0

威望

732

eV
积分
759

Level 4 (黑子)
跳转到指定楼层 Go to specific reply
楼主
发表于 Post on 2022-2-17 01:38:32 | 只看该作者 Only view this author |只看大图 View big picture 回帖奖励 |倒序浏览 Reverse view |阅读模式 Reading model
本帖最后由 大村驴 于 2022-2-17 01:46 编辑

这几天研究了一下二茂铁的力场构建，踩了一些坑，好在最后总算是做出来了。这里就把过程总结一下，当作后来人的教程。
先在论坛里搜了一下，关于如何用MCPB.py构建金属结合蛋白中的金属力场的问题，这个帖子已经总结过了（推测是我师弟）。但是关于二茂铁这种特殊成键方式的还没有，所以就把我做的这个例子也发一下（知乎上也同步发了，是本人）。

本文是根据 MCPB.py 的作者的两篇教程，在自己的二茂铁体系上操作得到的。推荐看看原文。
ambermd.org/tutorials/advanced/tutorial20/mcpbpy.htm
ambermd.org/tutorials/advanced/tutorial20/mcpbpy_heme.htm
注：本文改编自以上两个教程，如有错误之处，还请批评指教

使用 Amber18 版。仅用到 AmberTools 功能。

前言
引用

使用 MCPB.py 需引用原文

Pengfei Li and Kenneth M. Merz, Jr., "MCPB.py: A Python Based Metal Center Parameter Builder."J. Chem. Inf. Model.,2016,56, 599-604 One typo in the paper: should be step 2b (not 2c) is the default for step 2.

使用方法

  1. MCPB.py [flags] [options]
  2. 参数：
  3. -i 输入文件
  4. -s/--step 第几步
  5. 选项：
  6. -h/--help 帮助信息
复制代码

应用

MCPB.py 可以对金属蛋白和有机金属化合物系统进行参数化。MCPB.py 参考论文的 Supporting 中给出了两个例子，一个是针对金属蛋白的，另一个是针对有机金属化合物的。

用 MCPB 生成二茂铁的Fe参数

注：本文的所有操作都是在 Python 2.7 下完成的。用 conda 新建了一个虚拟环境，名为 py2，并指定 Python 2.7 为默认 python，然后在这里安装的 Amber，因此后面在执行 .py 文件的时候，默认都是在 py2 虚拟环境下运行的（文件首行指定的 python 环境也改成了 py2 中的 python 路径）。

关于使用什么 python 版本的问题，我和师弟得出的结论恰好相反。我在生成金属那里只能用 py2 环境，用 py3 就会出错。师弟则只能用 py3，用 py2 就会错。但这不重要，能用就行。

二茂铁简介

二茂铁是一种典型的夹心配合物，由一个二价铁和两个环戊二烯负离子（Cp）组成。二茂铁中心铁原子的氧化态为 +2，每个茂环带有一个单位负电荷，因此每个环含有 6 个 π 电子（五个 C 的 p 电子和一个额外负电荷），符合休克尔规则中4n+2电子数的要求（n为非负整数），故每个环都有芳香性。两个茂环以 π 电子与 Fe 成键，每个 Cp 拿出六个单电子与二价 Fe 构成六配位，形成 18 电子等排结构（每个环的 6 个电子 * 2，再加上二价铁离子的 6 个 d 电子），符合18电子规则，因此二茂铁非常稳定。

创建一个新目录存放所有文件。MCPB.py 在建模过程中需要一些中间文件。对于这里的残基中含有金属离子的系统，可以将 Fe 提取到独立的残基中进行建模。



首先对所建立的 Fc 结构在 TPSSh/def2-SVP 的条件下进行了几何优化，如上图所示。

1. 准备 PDB 和 mol2 文件
准备蛋白质的 PDB 文件，该文件应采用 PDB 3.0 版 的格式，并与残基的 AMBER 命名方案一致（例如 HID、HIE、HIP 用于 HIS）。 我们还需要为非标准残基（如金属离子、水和配体）准备 mol2 文件。

这一步相当重要，大部分错误都是在这里出现的，请仔细阅读以下内容！

对于 PDB 文件，应确保在任何残基中没有原子重名，这种情况在配体残基中可能出现。如果有，请在进行下一步之前手动更正。

将优化后的 FC 结构保存为 pdb，修正其原子名如下：

  1. REMARK   Accelrys Discovery Studio PDB file
  2. REMARK   Created:  2022-02-14T18:56:43Z
  3. HETATM    1 FE    FC A   1       6.103   6.652   9.882  1.00  0.00          FE  
  4. HETATM    2  C1   FC A   1       4.581   5.401   9.337  1.00  0.00           C  
  5. HETATM    3  C2   FC A   1       4.627   6.533   8.473  1.00  0.00           C  
  6. HETATM    4  C3   FC A   1       4.460   7.703   9.269  1.00  0.00           C  
  7. HETATM    5  C4   FC A   1       4.311   7.294  10.626  1.00  0.00           C  
  8. HETATM    6  C5   FC A   1       4.385   5.871  10.668  1.00  0.00           C  
  9. HETATM    7  C6   FC A   1       7.821   7.433   9.095  1.00  0.00           C  
  10. HETATM    8  C7   FC A   1       7.625   7.904  10.426  1.00  0.00           C  
  11. HETATM    9  C8   FC A   1       7.579   6.772  11.291  1.00  0.00           C  
  12. HETATM   10  C9   FC A   1       7.747   5.602  10.495  1.00  0.00           C  
  13. HETATM   11  C10  FC A   1       7.896   6.010   9.138  1.00  0.00           C  
  14. HETATM   12  H1   FC A   1       4.699   4.361   9.037  1.00  0.00           H  
  15. HETATM   13  H2   FC A   1       4.787   6.509   7.395  1.00  0.00           H  
  16. HETATM   14  H3   FC A   1       4.469   8.730   8.908  1.00  0.00           H  
  17. HETATM   15  H4   FC A   1       4.186   7.954  11.484  1.00  0.00           H  
  18. HETATM   16  H5   FC A   1       4.328   5.254  11.563  1.00  0.00           H  
  19. HETATM   17  H6   FC A   1       7.879   8.050   8.200  1.00  0.00           H  
  20. HETATM   18  H7   FC A   1       7.507   8.944  10.726  1.00  0.00           H  
  21. HETATM   19  H8   FC A   1       7.420   6.796  12.368  1.00  0.00           H  
  22. HETATM   20  H9   FC A   1       7.737   4.575  10.857  1.00  0.00           H  
  23. HETATM   21  H10  FC A   1       8.021   5.350   8.281  1.00  0.00           H  
  24. TER      22       FC A   1
  25. CONECT    2   12    6    3
  26. CONECT    3    2   13    4
  27. CONECT    4    3    5   14
  28. CONECT    5    4    6   15
  29. CONECT    6    2    5   16
  30. CONECT    7    8   11   17
  31. CONECT    8    7    9   18
  32. CONECT    9    8   10   19
  33. CONECT   10    9   20   11
  34. CONECT   11    7   10   21
  35. CONECT   12    2
  36. CONECT   13    3
  37. CONECT   14    4
  38. CONECT   15    5
  39. CONECT   16    6
  40. CONECT   17    7
  41. CONECT   18    8
  42. CONECT   19    9
  43. CONECT   20   10
  44. CONECT   21   11
  45. END
复制代码
MCPB.py 现在支持 80 多种金属离子（参见参考论文中的图 3）。对于所有这些金属离子和卤化物离子（如果是离子格式但不是中性格式），建议将大写的元素名称用于残基名称和原子名称。请确保金属离子的残基名称和原子名称全部大写（例如，此处均等于“FE”，而不是“Fe”、“fE”或“fe”），仅在这样 MCPB.py 会将其识别为金属离子。同时，请确保每个金属离子或卤化物离子（如果它是离子形式但不是中性形式）在 PDB 文件中作为独立残基单独处理。如果您的金属位点在配体残基中嵌入了金属离子，例如 PDB 文件中的 HEME 基团，请在 PDB 文件中将金属离子分离成一个独立的残基（也具有唯一的残基编号和原子编号） .同时，为了使 MCPB.py 能更好地识别原子，如果您的原始 PDB 文件中的原子名称以数字大写（适用于蛋白质或配体），请将这些原子名称更改为以大写元素符号开头的原子名称（例如，将 PDB 文件中的“2HA1”更改为“HA12”，将“1HAA”更改为“HAA1”），然后再执行以下步骤。

1）首先为非标准残基准备 PDB 和 mol2 文件：

请注意，正确获取每个 mol2 文件的总电荷很重要。 因为 MCPB.py 会将这些电荷加在一起以确定小型和大型模型的总电荷，并在 RESP 电荷拟合步骤中对大型模型应用总电荷的电荷约束。

① 处理配体

如果是从 PDB 数据库中获得的结构，需要为配体添加H。教程中使用的是 Amber 的 reduce 命令。也可以使用GaussView，Discovery Studio Visualizer 等。

二茂铁的两个配体就是两个 Cp 负离子。由于结构较为简单，直接用 Discovery Studio 提取相应的原子出来，并保存为新的 pdb 文件。

FE.pdb：
  1. REMARK   Accelrys Discovery Studio PDB file
  2. REMARK   Created:  2022-02-12T18:58:15Z
  3. HETATM    1 FE    FE A   1       6.103   6.652   9.882  1.00  0.00          FE  
  4. TER       2       FE A   1
  5. END
复制代码

CPA.pdb：
  1. REMARK   Accelrys Discovery Studio PDB file
  2. REMARK   Created:  2022-02-14T19:10:34Z
  3. HETATM    1  C1  CPA A   2       4.581   5.401   9.337  1.00  0.00           C  
  4. HETATM    2  C2  CPA A   2       4.627   6.533   8.473  1.00  0.00           C  
  5. HETATM    3  C3  CPA A   2       4.460   7.703   9.269  1.00  0.00           C  
  6. HETATM    4  C4  CPA A   2       4.311   7.294  10.626  1.00  0.00           C  
  7. HETATM    5  C5  CPA A   2       4.385   5.871  10.668  1.00  0.00           C  
  8. HETATM    6  H1  CPA A   2       4.699   4.361   9.037  1.00  0.00           H  
  9. HETATM    7  H2  CPA A   2       4.787   6.509   7.395  1.00  0.00           H  
  10. HETATM    8  H3  CPA A   2       4.469   8.730   8.908  1.00  0.00           H  
  11. HETATM    9  H4  CPA A   2       4.186   7.954  11.484  1.00  0.00           H  
  12. HETATM   10  H5  CPA A   2       4.328   5.254  11.563  1.00  0.00           H  
  13. TER      11      CPA A   2
  14. CONECT    1    6    5    2
  15. CONECT    2    1    7    3
  16. CONECT    3    2    4    8
  17. CONECT    4    3    5    9
  18. CONECT    5    1    4   10
  19. CONECT    6    1
  20. CONECT    7    2
  21. CONECT    8    3
  22. CONECT    9    4
  23. CONECT   10    5
  24. END
复制代码

CPB.pdb：
  1. REMARK   Accelrys Discovery Studio PDB file
  2. REMARK   Created:  2022-02-14T19:11:27Z
  3. HETATM    1  C6  CPB A   2       7.821   7.433   9.095  1.00  0.00           C  
  4. HETATM    2  C7  CPB A   2       7.625   7.904  10.426  1.00  0.00           C  
  5. HETATM    3  C8  CPB A   2       7.579   6.772  11.291  1.00  0.00           C  
  6. HETATM    4  C9  CPB A   2       7.747   5.602  10.495  1.00  0.00           C  
  7. HETATM    5  C10 CPB A   2       7.896   6.010   9.138  1.00  0.00           C  
  8. HETATM    6  H6  CPB A   2       7.879   8.050   8.200  1.00  0.00           H  
  9. HETATM    7  H7  CPB A   2       7.507   8.944  10.726  1.00  0.00           H  
  10. HETATM    8  H8  CPB A   2       7.420   6.796  12.368  1.00  0.00           H  
  11. HETATM    9  H9  CPB A   2       7.737   4.575  10.857  1.00  0.00           H  
  12. HETATM   10  H10 CPB A   2       8.021   5.350   8.281  1.00  0.00           H  
  13. TER      11      CPB A   2
  14. CONECT    1    2    5    6
  15. CONECT    2    1    7    3
  16. CONECT    3    2    8    4
  17. CONECT    4    3    5    9
  18. CONECT    5    1    4   10
  19. CONECT    6    1
  20. CONECT    7    2
  21. CONECT    8    3
  22. CONECT    9    4
  23. CONECT   10    5
  24. END
复制代码

然后我们可以使用 antechamber 生成两个 Cp 的 mol2 文件，这里我们为配体生成 -1 的 AM1-BCC 电荷。使用 GAFF 原子类型处理原子。
  1. antechamber -fi pdb -fo mol2 -i CPA.pdb -o CPA.mol2 -c bcc -pf y -nc -1
  2. antechamber -fi pdb -fo mol2 -i CPB.pdb -o CPB.mol2 -c bcc -pf y -nc -1
复制代码

生成的两个mol2文件如下所示：

CPA.mol2：
  1. @<TRIPOS>MOLECULE
  2. CPA
  3.    10    10     1     0     0
  4. SMALL
  5. bcc


  8. @<TRIPOS>ATOM
  9.       1 C1           4.5810     5.4010     9.3370 c2         2 CPA      -0.268000
  10.       2 C2           4.6270     6.5330     8.4730 c2         2 CPA      -0.268000
  11.       3 C3           4.4600     7.7030     9.2690 c2         2 CPA      -0.268000
  12.       4 C4           4.3110     7.2940    10.6260 c2         2 CPA      -0.268000
  13.       5 C5           4.3850     5.8710    10.6680 c2         2 CPA      -0.268000
  14.       6 H1           4.6990     4.3610     9.0370 ha         2 CPA       0.068000
  15.       7 H2           4.7870     6.5090     7.3950 ha         2 CPA       0.068000
  16.       8 H3           4.4690     8.7300     8.9080 ha         2 CPA       0.068000
  17.       9 H4           4.1860     7.9540    11.4840 ha         2 CPA       0.068000
  18.      10 H5           4.3280     5.2540    11.5630 ha         2 CPA       0.068000
  19. @<TRIPOS>BOND
  20.      1     1     2 1   
  21.      2     1     5 1   
  22.      3     1     6 1   
  23.      4     2     3 1   
  24.      5     2     7 1   
  25.      6     3     4 1   
  26.      7     3     8 1   
  27.      8     4     5 1   
  28.      9     4     9 1   
  29.     10     5    10 1   
  30. @<TRIPOS>SUBSTRUCTURE
  31.      1 CPA         1 TEMP              0 ****  ****    0 ROOT
复制代码

CPB.mol2：
  1. @<TRIPOS>MOLECULE
  2. CPB
  3.    10    10     1     0     0
  4. SMALL
  5. bcc


  8. @<TRIPOS>ATOM
  9.       1 C6           7.8210     7.4330     9.0950 c2         2 CPB      -0.268000
  10.       2 C7           7.6250     7.9040    10.4260 c2         2 CPB      -0.268000
  11.       3 C8           7.5790     6.7720    11.2910 c2         2 CPB      -0.268000
  12.       4 C9           7.7470     5.6020    10.4950 c2         2 CPB      -0.268000
  13.       5 C10          7.8960     6.0100     9.1380 c2         2 CPB      -0.268000
  14.       6 H6           7.8790     8.0500     8.2000 ha         2 CPB       0.068000
  15.       7 H7           7.5070     8.9440    10.7260 ha         2 CPB       0.068000
  16.       8 H8           7.4200     6.7960    12.3680 ha         2 CPB       0.068000
  17.       9 H9           7.7370     4.5750    10.8570 ha         2 CPB       0.068000
  18.      10 H10          8.0210     5.3500     8.2810 ha         2 CPB       0.068000
  19. @<TRIPOS>BOND
  20.      1     1     2 1   
  21.      2     1     5 1   
  22.      3     1     6 1   
  23.      4     2     3 1   
  24.      5     2     7 1   
  25.      6     3     4 1   
  26.      7     3     8 1   
  27.      8     4     5 1   
  28.      9     4     9 1   
  29.     10     5    10 1   
  30. @<TRIPOS>SUBSTRUCTURE
  31.      1 CPB         1 TEMP              0 ****  ****    0 ROOT
复制代码

检查生成的原子类型是 c2 和 ha，在 GAFF 力场中分别代表 sp2 杂化 C 和其连接的 H。

注意：非标准残基（本例中为配体）的 mol2 文件名应与其在PDB文件中的残基名称相同（字母和大写），后缀为“.mol2”，才能用于 MCPB.py。因此，后面的处理过程中，要始终确认两个 Cp 的残基名分别为 CPA 和 CPB。

注：保险起见，残基名建议不要带数字。

然后我们可以执行以下命令来获取配体的 frcmod 文件：
  1. parmchk2 -i CPA.mol2 -o CPA.frcmod -f mol2
  2. parmchk2 -i CPB.mol2 -o CPB.frcmod -f mol2
复制代码

这里使用 parmchk2，因为它比 parmchk 具有更强的性能。（但原教程做HEME的时候用的是parmchk，因为parmchk2给的参数不全。具体用的时候保证参数完整合理就行了，这里用的是parmchk2）

提醒：请确保 frcmod 文件中每个参数部分（BOND、ANGLE、DIHEDRAL、IMPROPER、NONBOND）后面只有一个空行。

frcmod 文件也将用于 leap 建模。

② 处理金属离子

对于金属离子，我们可以将 FE 离子复制到一个单独的 PDB 文件（FE.pdb）中，然后使用脚本 metalpdb2mol2.py 生成 mol2 文件。可能需要修改此脚本的第一行以指定要使用的 python 程序的路径。该脚本将包含在 AmberTools 2020 版本中，因此 AmberTools 2020 的用户无需下载。我们可以运行以下命令：
  1. ./metalpdb2mol2.py -i FE.pdb -o FE.mol2 -c 2
复制代码

这里“-c 2”表示离子的电荷为 2.0。金属离子的氧化态可以通过使用此选项或在之后修改 mol2 文件中的电荷信息来指定。

FE.mol2：
  1. @<TRIPOS>MOLECULE
  2. FE
  3.     1     0     1     0     0
  4. SMALL
  5. User Assigned Charge


  8. @<TRIPOS>ATOM
  9.       1 FE          6.1030    6.6520    9.8820 FE        1 FE       2.000000
  10. @<TRIPOS>BOND
  11. @<TRIPOS>SUBSTRUCTURE
  12.      1 FE          1 TEMP              0 ****  ****    0 ROOT
复制代码

注意：这个脚本使用了 Amber 自带的脚本，而Amber自带的脚本是用 Python 2 写的，与 Python 3 不兼容。

③处理水

本例不包含水。如果要在建模过程中保留结晶水，可以将其复制到一个 PDB 文件中，然后使用 tleap 将氢原子添加到水分子中，然后使用 antechamber 为水分子生成 mol2 文件。

2) 准备包含所有标准残基的 PDB 文件

本例只有二茂铁，不包含标准残基。

如果要做金属结合蛋白，可以使用网络服务器 H++ 将氢原子添加到 PDB 文件中。H++ 将在建模过程中删除非标准残基。除了生成带有氢原子的 PDB 文件外，它还将为系统生成 AMBER 拓扑和坐标文件。然后我们可以使用拓扑和坐标文件来生成 PDB 文件，因为这个 PDB 文件使用 AMBER 命名方案来命名残基（例如组氨酸是 HID、HIE、HIP，而不是 HIS）。具体参考开头的两篇官方文档。

3) 将所有 PDB 文件组合成一个 PDB 文件

为了与 PDB 文件中的通常顺序保持一致，一般先放置标准残基，然后是金属离子，然后是配体，最后是水（如果有的话）。

这里不建议用 Discovery Studio 合成，因为 DS 保存的时候会自动添加 CONECT 字段，往往导致结果不合理。可以使用下面命令：
  1. cat FE.pdb CPA.pdb CPB.pdb > FC_pre.pdb
复制代码

然后使用 pdb4amber 重新编号 PDB 文件。
  1. pdb4amber -i FC_pre.pdb -o FC.pdb
复制代码

注意删除所得 FC.pdb 里面的 TER 行。

现在我们有了用于下一阶段的 PDB 文件 FC.pdb 和 mol2 文件 FE.mol2，CPA.mol2 和 CPB.mol2。

FC.pdb：
  1. HETATM    1 FE   FE  A   1       6.103   6.652   9.882  1.00  0.00          FE  
  2. HETATM    2  C1  CPA A   2       4.581   5.401   9.337  1.00  0.00           C  
  3. HETATM    3  C2  CPA A   2       4.627   6.533   8.473  1.00  0.00           C  
  4. HETATM    4  C3  CPA A   2       4.460   7.703   9.269  1.00  0.00           C  
  5. HETATM    5  C4  CPA A   2       4.311   7.294  10.626  1.00  0.00           C  
  6. HETATM    6  C5  CPA A   2       4.385   5.871  10.668  1.00  0.00           C  
  7. HETATM    7  H1  CPA A   2       4.699   4.361   9.037  1.00  0.00           H  
  8. HETATM    8  H2  CPA A   2       4.787   6.509   7.395  1.00  0.00           H  
  9. HETATM    9  H3  CPA A   2       4.469   8.730   8.908  1.00  0.00           H  
  10. HETATM   10  H4  CPA A   2       4.186   7.954  11.484  1.00  0.00           H  
  11. HETATM   11  H5  CPA A   2       4.328   5.254  11.563  1.00  0.00           H  
  12. HETATM   12  C6  CPB A   3       7.821   7.433   9.095  1.00  0.00           C  
  13. HETATM   13  C7  CPB A   3       7.625   7.904  10.426  1.00  0.00           C  
  14. HETATM   14  C8  CPB A   3       7.579   6.772  11.291  1.00  0.00           C  
  15. HETATM   15  C9  CPB A   3       7.747   5.602  10.495  1.00  0.00           C  
  16. HETATM   16  C10 CPB A   3       7.896   6.010   9.138  1.00  0.00           C  
  17. HETATM   17  H6  CPB A   3       7.879   8.050   8.200  1.00  0.00           H  
  18. HETATM   18  H7  CPB A   3       7.507   8.944  10.726  1.00  0.00           H  
  19. HETATM   19  H8  CPB A   3       7.420   6.796  12.368  1.00  0.00           H  
  20. HETATM   20  H9  CPB A   3       7.737   4.575  10.857  1.00  0.00           H  
  21. HETATM   21  H10 CPB A   3       8.021   5.350   8.281  1.00  0.00           H  
  22. TER      21      CPB A   3
  23. END
复制代码

检查一下这个 FC.pdb 和开头的区别：由原来的一个残基变为三个残基，两个 Cp 的残基名大写，原子顺序有所调整。

2. 生成PDB、Gaussian、GAMESS-US和指纹建模文件

我们需要手动创建输入文件。这是输入文件 FC.in：
  1. software_version g09
  2. original_pdb FC.pdb
  3. group_name FC
  4. cut_off 2.5
  5. ion_ids 1
  6. ion_mol2files FE.mol2
  7. naa_mol2files CPA.mol2 CPB.mol2
  8. frcmod_files CPA.frcmod CPB.frcmod
  9. additional_resids 2 3
复制代码

FC.in 的参数决定了 MCPB.py 的行为。参数解释可以查看 AMBER 2016 参考手册 中的第 287-291 页，这里列出其中关键部分：
ion_ids The PDB atom ID of the complex's central metal ion. It should be an integer value, depending on how many metal ions are included in the metal complex.
ion_mol2files The name(s) of the ion(s) in the mol2 file(s) contained in the metal center. This can be one or several name(s), depending on how many kinds of ions are included in the metal center. The user can use antechamber to transfer the single ion PDB file to a mol2 file and then manually modify the atom type and the atomic charge of the metal ion in the mol2 file.
original_pdb This is the file name of the original PDB file, which should have only one chain. The PDB file should have hydrogen atoms and metal ions in it. Users are advised to use an application like pdb4amber to clean up the PDB file first. They are also advised to add the hydrogen atoms by using a webserver such as H++ before performing the modeling in MCPB.py.
cut_off The cutoff value is used to indicate there is a bond between the metal ion and the surrounding atoms. The unit is Angstroms. [The default is 2.8.]
group_name The group name the user has specified. The group name is the prefix for different kinds of modeling files e.g. PDB, fingerprint and Gaussian input files for different models. [The default is MOL.]
naa_mol2files The variable used to indicate non-amino acid mol2 file(s) in the metal complex if there are any nonstandard residue(s) in the metal complex. Examples of nonstandard residues include hydroxyl group and ligand molecules. For these residues, the user can use antechamber to generate the mol2 file(s) by first doing an AM1-BCC or HF/6-31G* RESP charge fit and then assigning an AMBER atom type (recommended for water or hydroxyl group) or a GAFF atom type (recommended for ligand). [The default value of this variable is the null list.]
large_opt A variable used to indicate whether to do an geometry optimization in the Gaussian input file. Three options are aviable: 0, 1, or 2. 0 means no optimization, 1 means only optimizing the hydrogen positions, 2 means full geometry optimization. [The default is 0.]
force_field The user-designated name of the force field. The current version supports ff94, ff99, ff99SB, ff03, ff03.r1, ff10, ff14ipq, ff14SB, ff14SB.redq, ff14SBonlysc, ff15ipq, ff15ipq-vac, and fb15. [The default is ff14SB.]
frcmod_files The variable used to indicate the parameter modification file(s) for the nonstandard residue(s) (e.g. frcmod file generated by parmchk for a ligand molecule) in the metal complex. It can be one name or several names seperated by space. [The default value of this variable is the null list.]
software_version The version of software the user used to perform the QM calculations. Three options are available, g03 (which represents Gaussian03), g09 (which represents Gaussian09) and gms (which represents GAMESS-US). [The default is g03.]
additional_resids Specify the residues' IDs for which you want to add to the models built by MCPB.py. For example, it may be a residue in the second layer which coordinates a metal bonded residue. It will increase the computational cost for QM calculations. [The default value of this variable is the null list.]

第一行 software_version 指定用于 QM 计算的软件版本，可以用 Gaussian09，Gaussian03 或 GAMESS-US。原教程中设置 large_opt 等于 1 来优化 RESP 计算的大模型的氢位置，因为所用的氢位置可能不稳定。我们改为 0，因为我们的模型已经进行过充分的优化了。使用 ff14SB 力场（默认）来执行建模。

additional_resids 指定了要向模型中添加的残基，这里指定添加 2 号和 3 号残基。（重要）如果不写这个，会导致第 2 步的时候 Cp 识别不出来。

执行以下命令：
  1. MCPB.py -i FC.in -s 1
复制代码

在此步骤中，将生成小型（small）、标准（standard）和大型（large）模型的 PDB 和指纹（fingerprint）文件。小型和大型模型的高斯输入文件（com）也会生成。大型模型的高斯输入文件将首先对氢原子执行优化，以纠正任何位置不佳的氢原子。标准模型指纹文件（FC_standard.fingerprint）包含标准模型中原子的原子类型信息（倒数第三列是原始原子类型，最后一列是最终分配的原子类型）。金属离子的原子类型（名称以 M 开头）和连接原子（名称以 Y/Z 开头）会自动重命名，以将它们与 AMBER 力场中的其他原子类型区分开来。您可以将 step_number 改为 1n 以生成指纹文件而不重命名任何原子类型，然后手动修改它（在此之前，请检查 AMBER 力场 parm*.dat 文件以确保您的金属离子和连接原子不存在原子类型）。指纹文件的末尾还包含金属离子与周围原子的连接信息，开头字母为“LINK”。用户还可以根据需要手动修改它们（例如，如果 QM 几何优化后有任何连接更改）。

考虑到 FC 体系的成键方式较为特殊，程序没有识别出 FE 与周围 C 的连接关系。由于环配位难以用传统动力学方法来描述，这里将形成配位键的方式修改为 FE 与每个 C 配位，故修改 FC_standard.fingerprint 如下（将两个环的 C 的原子类型改为 Y1 和 Z1，并在末尾增加一些LINK字段）：

参数含义：
残基编号-残基名-原子名 原子编号 原子原类型 -> 原子新类型
LINK 原子编号1-原子名1 原子编号2-原子名2
  1. 1-FE-FE 1 FE -> M1
  2. 2-CP1-C1 2 c2 -> Y1
  3. 2-CP1-C2 3 c2 -> Y1
  4. 2-CP1-C3 4 c2 -> Y1
  5. 2-CP1-C4 5 c2 -> Y1
  6. 2-CP1-C5 6 c2 -> Y1
  7. 2-CP1-H1 7 ha -> ha
  8. 2-CP1-H2 8 ha -> ha
  9. 2-CP1-H3 9 ha -> ha
  10. 2-CP1-H4 10 ha -> ha
  11. 2-CP1-H5 11 ha -> ha
  12. 3-CP2-C6 12 c2 -> Z1
  13. 3-CP2-C7 13 c2 -> Z1
  14. 3-CP2-C8 14 c2 -> Z1
  15. 3-CP2-C9 15 c2 -> Z1
  16. 3-CP2-C10 16 c2 -> Z1
  17. 3-CP2-H6 17 ha -> ha
  18. 3-CP2-H7 18 ha -> ha
  19. 3-CP2-H8 19 ha -> ha
  20. 3-CP2-H9 20 ha -> ha
  21. 3-CP2-H10 21 ha -> ha
  22. LINK 1-FE 2-C1
  23. LINK 1-FE 3-C2
  24. LINK 1-FE 4-C3
  25. LINK 1-FE 5-C4
  26. LINK 1-FE 6-C5
  27. LINK 1-FE 12-C6
  28. LINK 1-FE 13-C7
  29. LINK 1-FE 14-C8
  30. LINK 1-FE 15-C9
  31. LINK 1-FE 16-C10
复制代码

FC_large.fingerprint：
  1. 1-FE-FE
  2. 2-CPA-C1
  3. 2-CPA-C2
  4. 2-CPA-C3
  5. 2-CPA-C4
  6. 2-CPA-C5
  7. 2-CPA-H1
  8. 2-CPA-H2
  9. 2-CPA-H3
  10. 2-CPA-H4
  11. 2-CPA-H5
  12. 3-CPB-C6
  13. 3-CPB-C7
  14. 3-CPB-C8
  15. 3-CPB-C9
  16. 3-CPB-C10
  17. 3-CPB-H6
  18. 3-CPB-H7
  19. 3-CPB-H8
  20. 3-CPB-H9
  21. 3-CPB-H10
复制代码

3. 执行 Gaussian/GAMESS-US 计算

在这里，我们将使用在上一步中获得的高斯输入文件执行高斯计算。您可以根据需要更改输入文件中的参数（例如，CPU 数量、内存使用情况、方法、基组、隐式求解模型等）。确保高斯输入文件的多重度正确（对小型和大型模型）。错误的多重度可能导致错误结果、收敛失败或其他问题。

原教程用 B3LYP/6-31G* 来执行小型和大型模型的计算，这个级别其实不适合算过渡金属，另外 Pople 基组的准确性不如def系列，所以这里相应做了修改。

考虑到我们所用的初始模型已经进行了充分的优化，这里就不再做几何优化了。删去 FC_small_opt.com 文件的 opt 关键词，改用 TPSSh/def2-SVP 级别做单点计算（为 FC_small_fc.com 准备 chk 文件）；将 FC_small_fc.com 的级别也改为 TPSSh/def2-SVP。将 FC_large_mk.com 的级别改为 TPSSh/def2-TZVP。

可以使用不同层次的计算级别对小型和大型模型进行计算。但小模型的几何优化和力常数计算（振动分析）应使用相同的计算级别，以确保找到局部最小值。为了确保参数化合理，建议在优化后检查小模型的结构。这是因为一些与金属的配位键在几何优化后可能会断裂，在这种情况下，最终的参数化结果可能会出现问题，故应尝试不同的计算级别或其他策略（例如执行多步优化）以获得有意义的结果。这里在大模型的高斯计算中，我们不进行优化。可以通过修改 MCPB.py 输入文件中的“large_opt”变量来更改此设置。建议并行进行计算，因为这些作业可能需要一段时间。
  1. # 将 small 的高斯输入文件的计算级别改为 TPSSh/def2SVP, large 的改为 TPSSh/def2TZVP
  2. g16 < FC_small_opt.com > FC_small_opt.log
  3. g16 < FC_small_fc.com > FC_small_fc.log
  4. formchk FC_small_opt.chk FC_small_opt.fchk
  5. # Perform the Merz-Kollman RESP charge calculation for the large model:
  6. g16 < FC_large_mk.com > FC_large_mk.log
复制代码

小模型QM优化后，如果金属与配位原子的配位键断裂，可能需要改变方法或基组重新计算。如果这个问题经过多次试验仍不能解决，仍然可以使用 MCPB.py 建模，但使用另一种方法（如 PES 扫描）获取金属涉及键和角度参数，并手动将其添加到下一步生成的 PREFIX_mcpbpy.frcmod 文件中 并仔细检查最终的 tleap 输入文件以确保一切顺利。 如果您参数化FE金属位点并遇到上述问题，您可以使用经验方法来确定涉及键和角度参数的金属（只需将 --step 改为 2e 而不是 2）。

4. 执行最终建模

这里我们使用 Seminario 方法来生成力场参数。 其他选项可用：Z 矩阵（使用 step_number 2z）和经验（使用 step_number 2e）方法。 Empirical 方法不需要任何高斯计算来获得力常数（但仍然需要高斯计算来获得 RESP 电荷），但它仅支持当前版本的锌离子建模。
  1. MCPB.py -i FC.in -s 2
复制代码

注意观察程序的输出。如果没有键长键角的输出，证明成键没有被正确识别（没有原子与金属成键）。

我们可以得到一个 FC_mcpbpy.frcmod 文件，该文件将用于 leap 建模。

在这里，我们使用 ChgModB 执行 RESP 电荷拟合并为金属位点残基生成 mol2 文件。 其他选项也可用：ChgModA、ChgModC 和 ChgModD（分别为 3a、3c 和 3d）。
  1. MCPB.py -i FC.in -s 3a
复制代码

因为配体不是氨基酸，这里把参数改成 3a，即所有原子电荷都可变。具体看原文。

在这一步中，我们可以获得金属部位的 3 个 mol2 文件：FE1.mol2，CA1.mol2，CB1.mol2（CPA和CPB被重命名了）。这些文件中的电荷由 MK RESP 电荷拟合算法重新拟合。 我们将在 leap 建模中使用这些 mol2 文件。

生成 tleap 输入文件：
  1. MCPB.py -i FC.in -s 4
复制代码

至此程序处理完成。输出的金属位点残基文件为 FC_mcpbpy.pdb，同时还输出了 leap 输入文件 http://FC_tleap.in。之后我们可以使用 tleap 生成拓扑和坐标文件：
  1. tleap -s -f FC_tleap.in > FC_tleap.out
复制代码

现在我们生成了溶液中的拓扑和坐标文件（FC_solv.prmtop 和 FC_solv.inpcrd），以及真空中的拓扑和坐标文件（FC_dry.prmtop 和 FC_dry.inpcrd）。

可以通过 acpype 程序来将拓扑转换为 gromacs 拓扑的格式：
  1. acpype -p FC_dry.prmtop -x FC_dry.inpcrd
复制代码

5. 检查模型（对成功建模很重要）

在执行最小化和 MD 模拟之前，我们可以使用 VMD 和 cpptraj 对建模进行快速检查。

1）首先我们可以使用VMD来检查拓扑文件中是否存在与金属离子的配位键。

我们可以使用以下命令：
  1. vmd -parm7 FC_solv.prmtop -rst7 FC_solv.inpcrd
复制代码

得到结果如下图所示：


发现 FE 和 C9，C10 少了两个键没形成，而是定义在 [ pairs ] 字段。重复了几次也是这个问题，没有找到合适的理由，故手动将其全部修改为 [ bonds ]。此外，生成的键角参数有大问题，所有的 C-FE-C 键都需要重新测量键角并修改：
  1. [ bonds ]
  2. ;   ai     aj funct   r             k
  3.      1      2   1    2.0443e-01    7.3806e+04 ;     FE - C1   
  4.      1      3   1    2.0443e-01    7.3806e+04 ;     FE - C2   
  5.      1      4   1    2.0443e-01    7.3806e+04 ;     FE - C3   
  6.      1      5   1    2.0443e-01    7.3806e+04 ;     FE - C4   
  7.      1      6   1    2.0443e-01    7.3806e+04 ;     FE - C5   
  8.      1     12   1    2.0443e-01    7.3806e+04 ;     FE - C6   
  9.      1     13   1    2.0443e-01    7.3806e+04 ;     FE - C7   
  10.      1     14   1    2.0443e-01    7.3806e+04 ;     FE - C8   
  11.      1     15   1    2.0443e-01    7.3806e+04 ;     FE - C9   
  12.      1     16   1    2.0443e-01    7.3806e+04 ;     FE - C10   
  13. ...

  15. [ pairs ]
  16. ;   ai     aj    funct
  17.      2      9      1 ;     C1 - H3   
  18. ...

  20. [ angles ]
  21. ...
  22.      1     15     14      1    6.9600e+01    6.0208e+02 ;     FE - C9     - C8   
  23.      1     15     16      1    6.9600e+01    6.0208e+02 ;     FE - C9     - C10   
  24.      1     15     20      1    1.2526e+02    2.6920e+02 ;     FE - C9     - H9   
  25.      1     16     15      1    6.9600e+01    6.0208e+02 ;     FE - C10    - C9   
  26.      1     16     12      1    6.9600e+01    6.0208e+02 ;     FE - C10    - C6   
  27.      1     16     21      1    1.2526e+02    2.6920e+02 ;     FE - C10    - H10   
  28. ; 下面的 C-FE-C 键角是我重新测量过的, 原始输出全是100°, 这显然是不对的
  29.      2      1      3      1    40.81637    7.5262e+02 ;     C1 - FE     - C2   
  30.      2      1      4      1    68.67354    7.5262e+02 ;     C1 - FE     - C3   
  31.      2      1      5      1    68.65093    7.5262e+02 ;     C1 - FE     - C4   
  32.      2      1      6      1    40.79887    7.5262e+02 ;     C1 - FE     - C5   
  33.      2      1     12      1    139.17721    7.5262e+02 ;     C1 - FE     - C6   
  34.      2      1     13      1    179.96157    7.5262e+02 ;     C1 - FE     - C7   
  35.      2      1     14      1    139.23110    7.5262e+02 ;     C1 - FE     - C8   
  36.      2      1     15      1    111.36477    7.5262e+02 ;     C1 - FE     - C9   
  37.      2      1     16      1    111.33097    7.5262e+02 ;     C1 - FE     - C10   
  38.      2      3      4      1    1.2181e+02    5.7990e+02 ;     C1 - C2     - C3   
  39.      2      3      8      1    1.2043e+02    4.1756e+02 ;     C1 - C2     - H2   
  40.      2      6      5      1    1.2181e+02    5.7990e+02 ;     C1 - C5     - C4   
  41.      2      6     11      1    1.2043e+02    4.1756e+02 ;     C1 - C5     - H5   
  42.      3      1      4      1    40.81637    7.5262e+02 ;     C2 - FE     - C3   
  43.      3      1      5      1    68.67354    7.5262e+02 ;     C2 - FE     - C4   
  44.      3      1      6      1    68.67354    7.5262e+02 ;     C2 - FE     - C5   
  45.      3      1     12      1    111.33097    7.5262e+02 ;     C2 - FE     - C6   
  46.      3      1     13      1    139.17721    7.5262e+02 ;     C2 - FE     - C7   
  47.      3      1     14      1    179.96157    7.5262e+02 ;     C2 - FE     - C8   
  48.      3      1     15      1    139.23110    7.5262e+02 ;     C2 - FE     - C9   
  49.      3      1     16      1    111.36477    7.5262e+02 ;     C2 - FE     - C10   
  50.      3      2      6      1    1.2181e+02    5.7990e+02 ;     C2 - C1     - C5   
  51.      3      2      7      1    1.2043e+02    4.1756e+02 ;     C2 - C1     - H1   
  52.      3      4      5      1    1.2181e+02    5.7990e+02 ;     C2 - C3     - C4   
  53.      3      4      9      1    1.2043e+02    4.1756e+02 ;     C2 - C3     - H3   
  54.      4      1      5      1    40.81637    7.5262e+02 ;     C3 - FE     - C4   
  55.      4      1      6      1    68.67354    7.5262e+02 ;     C3 - FE     - C5   
  56.      4      1     12      1    111.36477    7.5262e+02 ;     C3 - FE     - C6   
  57.      4      1     13      1    111.33097    7.5262e+02 ;     C3 - FE     - C7   
  58.      4      1     14      1    139.17721    7.5262e+02 ;     C3 - FE     - C8   
  59.      4      1     15      1    179.96157    7.5262e+02 ;     C3 - FE     - C9   
  60.      4      1     16      1    139.23110    7.5262e+02 ;     C3 - FE     - C10   
  61.      4      3      8      1    1.2043e+02    4.1756e+02 ;     C3 - C2     - H2   
  62.      4      5      6      1    1.2181e+02    5.7990e+02 ;     C3 - C4     - C5   
  63.      4      5     10      1    1.2043e+02    4.1756e+02 ;     C3 - C4     - H4   
  64.      5      1      6      1    40.81637    7.5262e+02 ;     C4 - FE     - C5   
  65.      5      1     12      1    139.23110    7.5262e+02 ;     C4 - FE     - C6   
  66.      5      1     13      1    111.36477    7.5262e+02 ;     C4 - FE     - C7   
  67.      5      1     14      1    111.33097    7.5262e+02 ;     C4 - FE     - C8   
  68.      5      1     15      1    139.17721    7.5262e+02 ;     C4 - FE     - C9   
  69.      5      1     16      1    179.96157    7.5262e+02 ;     C4 - FE     - C10   
  70.      5      4      9      1    1.2043e+02    4.1756e+02 ;     C4 - C3     - H3   
  71.      5      6     11      1    1.2043e+02    4.1756e+02 ;     C4 - C5     - H5   
  72.      6      1     12      1    179.96157    7.5262e+02 ;     C5 - FE     - C6   
  73.      6      1     13      1    139.23110    7.5262e+02 ;     C5 - FE     - C7   
  74.      6      1     14      1    111.36477    7.5262e+02 ;     C5 - FE     - C8   
  75.      6      1     15      1    111.33097    7.5262e+02 ;     C5 - FE     - C9   
  76.      6      1     16      1    139.17721    7.5262e+02 ;     C5 - FE     - C10   
  77.      6      2      7      1    1.2043e+02    4.1756e+02 ;     C5 - C1     - H1   
  78.      6      5     10      1    1.2043e+02    4.1756e+02 ;     C5 - C4     - H4   
  79.     12      1     13      1    40.81637    7.5262e+02 ;     C6 - FE     - C7   
  80.     12      1     14      1    68.67354    7.5262e+02 ;     C6 - FE     - C8   
  81.     12      1     15      1    68.67354    7.5262e+02 ;     C6 - FE     - C9   
  82.     12      1     16      1    40.81637    7.5262e+02 ;     C6 - FE     - C10   
  83.     12     13     14      1    1.2181e+02    5.7990e+02 ;     C6 - C7     - C8   
  84.     12     13     18      1    1.2043e+02    4.1756e+02 ;     C6 - C7     - H7   
  85.     12     16     15      1    1.2181e+02    5.7990e+02 ;     C6 - C10    - C9   
  86.     12     16     21      1    1.2043e+02    4.1756e+02 ;     C6 - C10    - H10   
  87.     13      1     14      1    40.81637    7.5262e+02 ;     C7 - FE     - C8   
  88.     13      1     15      1    68.67354    7.5262e+02 ;     C7 - FE     - C9   
  89.     13      1     16      1    68.67354    7.5262e+02 ;     C7 - FE     - C10   
  90.     13     12     16      1    1.2181e+02    5.7990e+02 ;     C7 - C6     - C10   
  91.     13     12     17      1    1.2043e+02    4.1756e+02 ;     C7 - C6     - H6   
  92.     13     14     15      1    1.2181e+02    5.7990e+02 ;     C7 - C8     - C9   
  93.     13     14     19      1    1.2043e+02    4.1756e+02 ;     C7 - C8     - H8   
  94.     14      1     15      1    40.81637    7.5262e+02 ;     C8 - FE     - C9   
  95.     14      1     16      1    68.67354    7.5262e+02 ;     C8 - FE     - C10   
  96.     14     13     18      1    1.2043e+02    4.1756e+02 ;     C8 - C7     - H7   
  97.     14     15     16      1    1.2181e+02    5.7990e+02 ;     C8 - C9     - C10   
  98.     14     15     20      1    1.2043e+02    4.1756e+02 ;     C8 - C9     - H9   
  99.     15      1     16      1    40.81637    7.5262e+02 ;     C9 - FE     - C10   
  100.     15     14     19      1    1.2043e+02    4.1756e+02 ;     C9 - C8     - H8   
  101.     15     16     21      1    1.2043e+02    4.1756e+02 ;     C9 - C10    - H10   
  102.     16     12     17      1    1.2043e+02    4.1756e+02 ;    C10 - C6     - H6   
  103.     16     15     20      1    1.2043e+02    4.1756e+02 ;    C10 - C9     - H9   
  104. ...
复制代码

二面角项虽然也缺失，但包含 FE 的二面角力常数全是 0，推测这个并不重要，不加了。

注：如果键角项没有补充，可能会造成 FE 原子在模拟过程中倒向一侧的 Cp。

2）其次，我们可以使用 cpptraj 检查原子编号问题。

使用以下命令：
  1. cpptraj -p FC_solv.prmtop
复制代码

如果没有报错，我们可以继续进行最小化和MD模拟。

如果报告如下错误：

  • Error: Atom XXX was assigned a lower molecule # than previous atom. This can
  • Error: happen if 1) bond information is incorrect or missing, or 2) if the
  • Error: atom numbering in molecules is not sequential. If topology did not
  • Error: originally contain bond info, 1) can potentially be fixed by
  • Error: increasing the bondsearch cutoff offset (currently 0.200). 2) can be
  • Error: fixed by either using the 'fixatomorder' command, or using
  • Error: the 'setMolecules' command in parmed.py.
  • Error: Could not determine molecule information for FC_solv.prmtop.
  • Error: SetSolventInfo: No molecule information.
  • Error: Could not determine solvent information for FC_solv.prmtop.

我们可以使用 cpptraj 来修复原子顺序：

fixatord.in 文件内容：
  1. fixatomorder outprefix fixed
  2. trajout restart fixed.FC_solv.inpcrd
  3. run
  4. quit
复制代码

命令：
  1. cpptraj -p FC_solv.prmtop -c FC_solv.inpcrd -i fixatord.in > fixatord.out
复制代码

然后我们可以试试下面的命令看看问题是否解决了：
  1. cpptraj -p fixed.FC_solv.prmtop
复制代码

如果没有如上所示的错误，我们可以使用 fixed.FC_solv.prmtop 和 fixed.FC_solv.inpcrd 进行建模。

3）最后，我们可以使用 ParmEd 来检查金属位点参数：

输入文件：mcpbpy_parmed.in
  1. printBonds :FE1
  2. printAngles :FE1
  3. printDihedrals :FE1
  4. printDetails :FE1
复制代码

其中我们有“:FE1”作为亚铁离子的掩码。更多信息可以在 AMBER 手册的 Atom 和 Residue Selections 一章的 Amber Masks 部分中找到。 这里我们使用 ParmEd 及以上输入文件来检查参数：
  1. parmed -i mcpbpy_parmed.in -p FC_solv.prmtop
复制代码

通常金属离子相关参数遵循以下标准：

  • 金属离子与其配位原子的键力常数小于 200 kcal/(mol*Angstrom^2)，平衡键距离小于 2.8 埃；
  • 与金属离子相关的角力常数通常小于100 kcal/(mol*Rad^2)，而平衡角值大于100度；
  • 对于涉及金属的二面角，所有或大部分二面角势垒为零；
  • 金属离子的RESP电荷小于其氧化态，通常甚至小于+1；
  • 一种金属离子的 LJ 半径通常大于 1.0 埃。

我们可以看到 ParmEd 的输出结果与我们上面列出的条件一致，因此我们可以继续使用拓扑和坐标文件进行相关建模。

如果在最小化或分子动力学模拟后金属位点结构扭曲或在这些过程中遇到错误，这可能意味着正常成键模型无法描述这种情况下的金属位点。例如，对于有两个水分子连接到同一个金属离子的情况，由于一种水中的氧和另一种水中的氢之间的 1-4 VDW 相互作用为零，而 1-4 静电它们之间的相互作用是非零且有吸引力的。应对策略是删除 LEaP 建模中的相关键，以非键的方式处理键，同时在最小化和分子动力学模拟中对这些键使用谐振约束。以这种方式，两个水分子之间的完全非键相互作用将被视为解决问题。用户可以以这种方式处理所有的配位键，或者以成键方式处理其中一些，而以这种方式处理其他配位键（混合成键/约束非键模型）。可以在此处找到相关教程：受约束的非粘合模型教程）。

结果测试

用该拓扑跑一个MD，结果如下图所示


比较正常。


评分 Rate

参与人数
Participants 3		eV +13 收起 理由
Reason
wildon + 4 好物！
量化小菜鸟 + 4
puzhongji + 5 你太可爱

查看全部评分 View all ratings

回复 Reply

举报 Report

心向暖阳

68

帖子

0

威望

1637

eV
积分
1705

Level 5 (御坂)
2#
发表于 Post on 2022-2-17 09:31:43 | 只看该作者 Only view this author
想我当初折腾MPCB.py的时候，不记得花了多久才搞清楚怎么用
麻了 麻了
回复 Reply

举报 Report

sobereva

5万

帖子

99

威望

5万

eV
积分
112351

管理员

公社社长
3#
发表于 Post on 2022-2-17 12:52:01 | 只看该作者 Only view this author
用Sobtop产生过渡金属配合物的拓扑文件比MCPB.py方便得多得多
http://bbs.keinsci.com/thread-27803-1-1.html
北京科音自然科学研究中心（http://www.keinsci.com）致力于计算化学的发展和传播，长期开办高质量的各种计算化学类培训：初级量子化学培训班、中级量子化学培训班、高级量子化学培训班、量子化学波函数分析与Multiwfn程序培训班、分子动力学与GROMACS培训班、CP2K第一性原理计算培训班，内容介绍以及往届资料购买请点击相应链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的高速路！培训各种常见问题见《北京科音办的培训班FAQ》。
欢迎加入“北京科音”微信公众号获取北京科音培训的最新消息、避免错过网上有价值的计算化学文章！
欢迎加入人气非常高、专业性特别强的综合性理论与计算化学交流QQ群“思想家公社QQ群”：1号：18616395，2号：466017436，3号：764390338，搜索群号能搜到哪个说明目前哪个能加，合计9000人。北京科音培训班的学员在群中可申请VIP头衔，提问将得到群主Sobereva的最优先解答。
思想家公社的门口Blog：http://sobereva.com（发布大量原创计算化学相关博文）
Multiwfn主页：http://sobereva.com/multiwfn（十分强大的量子化学波函数分析程序）
ResearchGate：https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!
回复 Reply

举报 Report

大村驴

27

帖子

0

威望

732

eV
积分
759

Level 4 (黑子)
4#
 楼主 Author| 发表于 Post on 2022-2-17 13:09:17 | 只看该作者 Only view this author
sobereva 发表于 2022-2-17 12:52
用Sobtop产生过渡金属配合物的拓扑文件比MCPB.py方便得多得多
http://bbs.keinsci.com/thread-27803-1-1.h ...

谢谢老师提醒！我学习一下Sobtop
回复 Reply

举报 Report

jiangyy

6

帖子

0

威望

117

eV
积分
123

Level 2 能力者
5#
发表于 Post on 2022-6-16 15:41:24 | 只看该作者 Only view this author
你好，我在做酮戊二酸依赖的Fe家族的酶时，用MCPB生成的frcmod文件产生的键长键角的力常数都为0，是怎么回事呢？谢谢
BOND
M1-Y4    0.0    1.7140      Created by Seminario method using MCPB.py
M1-Y5    0.0    1.4436      Created by Seminario method using MCPB.py
M1-Y6    0.0    2.9193      Created by Seminario method using MCPB.py
Y1-M1    0.0    2.6213      Created by Seminario method using MCPB.py
Y2-M1    0.0    2.7930      Created by Seminario method using MCPB.py
Y3-M1    0.0    2.8003      Created by Seminario method using MCPB.py
回复 Reply

举报 Report

chw985550192

72

帖子

0

威望

400

eV
积分
472

Level 3 能力者
6#
发表于 Post on 2023-5-4 10:52:33 | 只看该作者 Only view this author
本帖最后由 chw985550192 于 2023-5-4 11:00 编辑

你好，你的这个例子是形成二茂铁的prmtop和inpcrd文件。如果我的体系中既有二茂铁，又有蛋白，那蛋白的prmtop和inpcrd文件可以用如FF14SB力场形成，那岂不是有两个prmtop文件和两个inpcrd文件，后边跑MD，要怎么操作呢，谢谢。
回复 Reply

举报 Report

返回列表 Return 发新帖

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2024-11-23 10:38 , Processed in 0.216281 second(s), 25 queries , Gzip On.

快速回复 返回顶部 返回列表 Return to list