金属中心参数构建＆Python金属位点建模工具

金属中心参数构建（MCPB）

1. 引言

金属中心参数构建程序提供了一种快速构建并验证金属蛋白分子力学模型参数的方法。它使用键连加上静电模型的方法扩展已经存在的对势力场。它由佛罗里达大学Kenneth Merz Jr实验室的Martin Peters开发的。MCPB在参考文献344中有完整的描述。

为什么使用分子力学模型，或者更确切的说是用键连加上静电模型，去构建金属蛋白体系更好呢？目前还不能使用基于QM或者QM/MM的方法去实现结构/函数和动力学的问题，因为无法获得参数或者系统大小。含有Zn、Cu、Ni、Fe和Pt的系统的力场已经可以通过键连和静电模型来构建。

把金属引入蛋白质的力场中看起来是一项艰巨的任务，因为有太多的量子力学哈密顿量、基组和电荷模型需要去选择，以构建参数。它需要手动构建，而且对特定的金属蛋白没有广泛的验证。MCPB开发出来后消除了后述弊端，为了测试各种方法、基组和电荷模型，在构建金属蛋白力场的过程中构建了一个框架。

MCPB程序使用MTK++应用程序接口。需要更多MTK++和MCPB的信息，参见MTK++的手册：$AMBERHOME/AmberTools/src/mtkpp/-doc/MTKpp.pdf。关于金属蛋白更为详尽的描述和MCPB相关的理论可见MTK++手册的10到12章节。

2. 运行MCPB

MCPB有两个命令行参数。一个是控制文件，它是必须的，使用-i来调用。另外一个是日志文件，它是可选的，使用-l来调用。可以使用-f来获得MCPB的所有参数设置列表。

MCPB:Semi-automated tool for metalloprotein parametrization

usage: MCPB[flags] [options]

options:

-i script file

-l log file

flags:

-h help

-f functionlist

使用MCPB构建金属蛋白的完整详情可以在MTK++手册的15.10节找到。这个例子描述了2-Zn体系活性位点的参数化过程。（PDB ID: 1AMP）。参数化被分为多部，由于MCPB操作依赖于外部程序的输出，比如Gaussian和RESP。大多数步骤可以用MCPB操作，但是有一些需要使用者的输入和指令。

Python金属位点建模工具（pyMSMT）

1. 引言

1) MCPB.py：金属中心参数生成器（MCPB）的Python版本。 MCPB.py支持各种金属离子（80多种金属离子，电荷/氧化态从+1到+8），不同的AMBER力场（ff94，ff99，ff99SB，ff03，ff03.r1，ff10，ff14ipq ，ff14SB，ff14SB.redq，ff14SBonlysc，ff15ipq，ff15ipq-vac，fb15，GAFF和当前版本中的GAFF2），不同的参数化方法（Seminario，Z-矩阵和经验化），以及不同的模型（金属离子的成键和非成键模型）。它可以帮助金属蛋白和有机金属化合物的参数化，且工作流程效率更高。上文提及的MCPB版本中的许多建模过程会自动移植到MCPB.py中（MCPB.py使用的步骤少于10个，脚本比MCPB少）。该程序的操作说明书由Li和Merz发表[345]。主要方案和参数都是基于Merz等人发表的论文[94-97,334]。

2) IPMach.py：Python版本的离子参数化工具。 IPMach.py可以在很大程度上帮助离子的12-6 LJ模型和12-6-4 LJ模型参数化。它可以根据给定的水合自由能或离子-氧原子距离自动将12-6 LJ模型参数化，它也可以根据给定的水合自由能或离子-氧原子距离自动将12-6-4 LJ模型参数化。

3) PdbSearcher.py：Pdbsearcher的Python版本。 由于金属离子的PDB命名方案的兼容性更好，PdbSearcher.py更好地支持PDB文件中金属中心的自动识别。

4) OptC4.py：一个使用AMBER拓扑和坐标文件优化12-6-4势能函数中的C4项的程序。它可以自动优化与金属位点相关的C4项，以更好地重现实验测得的结构。 它以金属位点键长、键角和二面角的无量纲误差之和为标准（其中键、角和二面角具有不同的权重）。 对于每个优化周期，将通过OpenMM程序[346,347] 使得结构最小化，然后计算无量纲误差之和。 它需要6.3版本的OpenMM和当前版本中安装的SciPy软件包。

5) CartHess2FC.py：基于Seminario方法使用笛卡尔型的 Hessian矩阵计算力常数的程序。 它可以使用Gaussian程序中的fchk文件或者GAMESS程序中的日志文件计算体系中所有的键和角度的力常数，因为这些文件中包含了笛卡尔型的Hessian矩阵。

6) espgen.py：antechamber程序中espgen的Python版本。它可以从Gaussian程序输出文件或GAMESS程序日志文件中提取静电势信息。它支持Gaussian03，Gaussian09和GAMESS程序的输出文件。

7) ProScrs.py：“蛋白质剪刀”程序，用于将蛋白质片段切割和封端成簇。

8) car_to_files.py：基于car文件生成mol2和PDB文件的程序。 该功能是为在AMBER中使用INTERFACE力场的用户设计的，可以在https://bionanostructures.com/interface-md/上查看相关信息。

9) amb2chm_psf_crd.py：将AMBER中的prmtop和inpcrd文件转化为CHARMM中的PSF和CRD文件的程序，该功能是为了在CHARMM中使用AMBER力场的用户设计的。

10) amb2chm_par.py：将AMBER中的dat/frcmod文件转化为CHARMM中的PAR文件的程序。它可以一步整合多个AMBER的dat和/或frcmod文件为一个CHARMM的par文件。该功能是为了在CHARMM中使用AMBER力场的用户设计的。

11) mol2rtf.py：将mol2文件转化为CHARMM中的RTF文件的程序。该功能是为了在CHARMM中使用AMBER力场的用户设计的。

2. 使用方法

下面概述三个应用程序的使用方法和选项：

2.1 MCPB.py

Usage: MCPB.py -i input_file -s/--step step_number

[--logf Gaussian/GAMESS-US output logfile]

[--fchk Gaussian fchk file]

Options:

-h, --help Showthis help message and exit

-i INPUTFILE Inputfile name

-s STEP, --step=STEP Stepnumber

--logf Gaussian/GAMESS-USoutput logfile

--fchk Gaussianfchk file

以下是对input_file中变量的介绍：

（注：input_file中不应该有空行，值或参数应该紧跟变量，并由空格分隔。）

必选变量：

ion_ids 复合物中心金属离子的PDB原子ID。如果只有一个金属离子在金属位点，你需要将PDB原子ID放在变量后面。如果有多个金属离子在金属位点，你需要将PDB原子ID们以空格分开放在变量名后面。每个PDB原子ID应该是一个整数值。

ion_info 只有非键模型需要该变量，不需要重新设置残基电荷（步骤编号4n2）。每种金属离子都需要合计四个数据点：1）PDB文件中金属离子的残基名称；2）PDB文件中金属离子的原子名称；3）金属离子的元素符号；4）金属离子的电荷（或氧化态，其必须是整数）。举例说明：ZN ZN Zn 2（前两个是PDB文件中Zn2 +离子的残基和原子名称，第三个是它的元素符号，最后一个是它的电荷）。

ion_mol2files mol2文件中金属中心包含的离子名称。这可以是一个或多个名称，取决于金属中心包含多少种离子。用户可以使用antechamber将单个离子PDB文件转化为mol2文件，然后手动修改mol2文件中的原子类型和金属离子的原子电荷。

original_pdb 这是原始PDB文件的文件名，它应该只有一个链。 PDB文件中应该包含氢原子和金属离子。建议用户使用像pdb4amber这样的应用程序先清理PDB文件。也建议他们在MCPB.py中执行建模之前使用H ++等网络服务器添加氢原子。

可选择的变量：

add_bonded_pairs 指定想要通过MCPB.py添加到模型构建中的键合原子对。 默认情况下，MCPB.py只检测Metal-N / O / S / F / Cl / Br / I键，如果在金属位置有其他种类的金属连接键（例如Metal-C键），两个原子之间应该有破折号。例如，如果金属位置中有两个Metal-C键，并且金属的原子序号为1001，而两个碳原子的原子序号分别为1320和1380，则可以在输入文件中使用添加以下行： “add_bonded_pairs 1001-1320 1001-1380”，或者使用分开的行比如“add_bonded_pairs 1001-1320”和“add_bonded_pairs 1001-1380”。 [这个变量的默认值是空列表]

add_redcrd 指定是否为小模型的Gaussian计算添加了附加冗余坐标。该选项是为Z矩阵方法设计的。如果您使用Seminario方法，则可以忽略此选项。默认情况下，Gaussian程序使用冗余内坐标执行几何优化，默认的内坐标可能不包含金属-配体配位键和角度（表示包含金属的键角）。如果这些键和角度不包含在内的话，则在使用Z矩阵方法时，用户无法生成相关的力常数。 0表示不为金属-配体配位键和角度添加冗余坐标。1表示为金属-配体配位键和角度添加冗余坐标以优化小模型。这样，后面的力常数计算将在读取以前生成的chk文件时使用与优化程序相同的冗余内部坐标。选择1时可能会出现几何优化不收敛的情况，那么可以选择2。 2意味着只为小模型计算力常数去进行几何优化。建议用户使用在选项1时，出现几何优化不收敛情况时，使用此选项。这样，Gaussian程序将以默认方式执行几何优化，但Z矩阵的力常数将基于更新的冗余坐标。 [默认值为0]

additional_resids 指定要添加到MCPB.py构建的模型中的残基的ID。例如，它可能是协调金属连接的残基的第二层中的一个残基。这会增加QM计算的计算成本。 [此变量的默认值为空列表]

anglefc_avg 用于指示是否根据Seminario方法中选择子矩阵的不同方式导出的角度力常数的平均值。根据Seminario方法，在参数推导过程中，选择两个原子的子矩阵有A-B和B-A两种方法。以不同方式选择子矩阵获得的角度力常数可能没有很大差异。有两个选项可用：0或1。0表示不使用平均值，使用默认方式选择子矩阵。1表示取不同形式的平均值来选择子矩阵。 [默认值为0]

bondfc_avg 用于指示是否根据Seminario方法中选择子矩阵的不同方式导出的键力常数的平均值。根据Seminario方法，在参数推导过程中，选择两个原子的子矩阵有A-B和B-A两种方法。以不同方式选择子矩阵获得的键力常数可能没有很大差异。有两个选项可用：0或1。 0表示不使用平均值，使用默认方式选择子矩阵。 1表示用不同的方式平均选择子矩阵。 [默认值为0]

chgfix_resids 指定在拟合期间电荷将要被固定的残基ID。固定的电荷值是从建模过程中使用的mol2文件中引用的。[此变量的默认值为空列表]

cut_off cutoff值用于表示金属离子与周围多长距离内的原子之间存在键，单位是埃。 [默认值为2.8]

force_field 用户指定力场的名称。当前版本支持ff94，ff99，ff99SB，ff03，ff03.r1，ff10，ff14ipq，ff14SB，ff14SB.redq，ff14SBonlysc，ff15ipq，ff15ipq-vac和fb15。 [默认值为ff14SB]

frcmod_files 用于指示金属络合物中非标准残基（如parmchk为配体分子生成的frcmod文件）的参数修改文件的变量。它可以是由空间分隔的一个名字或几个名字。 [这个变量的默认值是空列表]

gaff 用于表示在建模过程中使用GAFF力场的变量。 0表示不使用GAFF，1表示使用GAFF，2表示使用GAFF2。[默认值是1]

group_name 用户指定的组名。组名是各种建模文件的前缀，例如不同类型的PDB，fingerprint和Gaussian输入文件。[默认为MOL]

ion_paraset 用户设置指定的用于构建非键模型的离子参数集。（该选项对键合模型中的金属离子VDW参数没有影响，其中作者已经为不同离子选择了某些VDW参数集）。这个变量有四个选项：HFE，CM，IOD和12_6_4（提示：数字之间的“_“是下划线）。 如果您使用12-6 Lennard-Jones非键合模型，推荐设置为：HFE设置为+1和-1离子，CM设置为+2离子，IOD设置为+3和+4离子。 它们也是这些金属离子的默认设置。

large_opt 用于指示是否在Gaussian输入文件中进行几何优化。有三种选择：0,1或2。0表示不优化，1表示仅优化氢位置，2表示完全几何优化。 [默认值为0.]

lgmodel_chg 确定大模型的总电荷。 [电荷的默认值将由程序自动分配，但不保证是正确的。建议仔细检查运行的Gaussian / GAMESS-US程序。如果不正确，可以将此选项和正确的电荷添加到MCPB.py输入文件中，并重新生成建模文件]

lgmodel_spin 确定大模型的自旋状态。 [自旋的默认值将根据电子数量由程序自动分配为1或2。这并不保证是正确的。建议仔细检查运行的Gaussian / GAMESS-US程序。如果不正确，可以将此选项和正确的自旋状态添加到MCPB.py输入文件中，并重新生成建模文件]

naa_mol2files 用于指示金属络合物中非氨基酸mol2文件中是否有任何非标准残基。非标准残基的例子包括羟基和配体分子。对于这些残基，用户可以通过antechamber首先进行AM1-BCC或HF / 6-31G * RESP电荷拟合，然后分配AMBER原子类型（推荐用于水或羟基）或GAFF原子类型（推荐用于配体）。 [这个变量的默认值是空列表]

scale_factor 根据QM方法指定力常数推导的频率比例因子。这个比例因子会将QM计算中确定所有的键合力和角度力常数进行缩放。提示：力常数比例因子通常等于频率比例因子的平方。例如，如果使用HF / 6-31G *理论级别进行计算，并且其频率比例因子为0.9，则需要使用的力常数比例因子为0.9 ^ 2 = 0.81。 [默认值为1.0（不执行缩放）]

smmodel_chg 确定小模型的总电荷。 [电荷的默认值将由程序自动分配，但不保证是正确的。建议仔细检查运行的Gaussian / GAMESS-US程序。如果不正确，可以将此选项和正确的电荷添加到MCPB.py输入文件中，并重新生成建模文件]

smmodel_spin 确定小模型的自旋状态。[自旋的默认值将根据电子数量由程序自动分配为1或2。这并不保证是正确的。建议仔细检查运行的Gaussian / GAMESS-US程序。如果不正确，可以将此选项和正确的自旋状态添加到MCPB.py输入文件中，并重新生成建模文件]

software_version 用户执行QM计算所用的软件版本。有三个选项可用，g03（代表Gaussian03），g09（代表Gaussian09）和gms（代表GAMESS-US）。 [默认为g03]

sqm_opt 用于指示在使用Gaussian执行计算之前，使用AmberTools中的SQM对侧链和/或大模型进行模拟的变量。请注意：如果选择1,2或3，建模过程的第一步将花费额外的时间（侧链模型大概要花费数分钟和大模型大概要花费数小时）。 [默认值是0.]

0- 表示没有使用SQM

1- 表示仅对侧链模型进行优化

2- 表示仅对大模型进行优化

3- 表示对侧链模型和大模型进行优化

water_model 用于分子建模的用户指定水模型。选项是TIP3P，SPCE和TIP4PEW。 [默认为TIP3P]

xstru 指定原始PDB文件中的结构是否用于生成frcmod文件中的平衡键距离和角度值。 0表示不使用，但使用QM优化结构。 1表示使用。 [默认值是0.]

以下是step_number变量的说明：

以下是step_number的选项：

对于第1步，有三个选项：1a（默认，与指定1相同），1m和1n。

对于第2步，有四个选项：2b，2e，2s（默认，与指定2相同）和2z。

对于第3步，有四个选项：3a，3b（默认，与指定3相同），3c和3d。

对于第4步，有三个选项：4b（默认，与指定4相同），4n1和4n2。

以下是在建模过程中使用的步骤的详细说明：

Step 1 为不同模型（例如侧链，标准和大模型）生成建模文件（例如PDB，fingerprint和Gaussian输入文件）。有三个选项可用，其解释如下所示。默认是1a。

•1a - 用于自动重命名标准fingerprint文件中中心金属离子和周围键合原子的原子类型。

•1m - 用于仅将中心金属离子的原子类型自动重命名为标准fingerprint文件中的AMBER原子离子原子类型样式。

•1n - 用于生成标准fingerprint文件，无需重命名原子类型。用户可以在标准fingerprint文件中手动重命名金属离子及其结合原子的原子类型。

请注意：在使用Step1和Step2之间，应使用Gaussian输入文件对侧链模型（计算力常数）和大模型（进行RESP电荷计算）进行Gaussian计算（如果需要）。在计算之前，用户可以根据自己的喜好更改Gaussian输入文件中的参数（例如计算方法，基本设置等）。完成此过程后，用户可以转到步骤2。

Step 2 用于生成建模的frcmod文件。在此步骤中，预先生成一个frcmod文件（文件名为pre.frcmod名称）。该文件包含除与金属离子相关的键长和键角参数外的所有参数。稍后，将生成最终的frcmod文件，其中将包含所有参数。有三种方法可供选择：经验型，Seminario和Z矩阵。每种方法都会生成与金属离子相关的键长和键角参数。如果没有QM优化结构，你也可以生成一个frcmod文件，其与金属相关的键长和键角参数为零（如步骤2b），然后稍后手动修改以供进一步使用。默认是2s（Seminario方法）。

•2b - “空白”方法。生成一个frcmod文件，其中与金属相关的键长和键角的平衡值和力常数为零。如果在MCPB.py输入文件中使用选项“xstru 1”，它将根据原始PDB结构和力常数为零生成平衡值。用户可以通过修改生成的frcmod文件，手动分配键长和键角的参数，以供进一步使用。

•2e - 经验法，[348]可以高效地生成与金属离子相关的键长和键角的参数，而无需进行Gaussian计算。当前版本只支持Zn2+离子。

•2s -Seminario方法[349]根据从量子计算中获得的笛卡尔型Hessian矩阵的子矩阵生成力场参数。 这种方法需要一个Gaussianfchk文件（可以通过使用Gaussian程序中的formchk命令从chk文件生成）。提醒：需要几何优化和力常数计算程序来生成最终的chk文件，之后用fchk文件通过Seminario方法进行力常数计算。

•2z - Z矩阵方法通过使用从量子计算中获得的笛卡尔Hessian矩阵生成力场参数。 该方法需要在几何优化和力常数计算之后使用力常数Gaussian输出文件（通常命名为log文件）。

Step 3 用于执行RESP电荷拟合并生成金属离子络合物内残基的mol2文件。这一步有几种适合的方案。下面显示了四个选项。默认值为3b，因为在Peters等人[328]的工作中，Seminario / ChgModB被确定为的最佳组合。提醒：chgfix_resids变量在此过程中有效，如果指定了该变量，则也会使用电荷限制，与以下选项之一共同作用。

•3a -允许连接残基中的原子的所有电荷没有任何限制地改变。

•3b -根据用户选择的力场，限制配体残基中重骨架原子的电荷。

•3c -根据用户选择的力场，限制连接残基中骨架原子（重原子和氢原子）的电荷。

•3d -根据用户选择的力场，限制连接残基中骨架原子（重原子和氢原子）和C β原子的电荷

Step 4 生成leap输入文件。 默认值是4b。

•4b - 生成键合模型的leap输入文件。

•4n1 - 生成非键合模型的leap输入文件并重新调整连接残基的电荷。

•4n2 - 为非键合模型生成leap输入文件，而不必重新调整连接残基的电荷。

下面是对参数化过程的一些建议：

1) 对于成键模型来说，以下步骤通常是必不可少的（4步）：1a/1n→2e/2s/2z→3a/3b/3c/3d→4b

2) 对于需要重新拟合电荷的非键模型，用户可以使用以下流程（3步）：

1m→3a/3b/3c/3d→4n1

3) 对于常规的非键模型建模（不需要拟合任何电荷），用户通常只需一步就可以进行建模（1步）：4n2

以下是对logf和fchk变量的解释：

这些变量是可选的，如果提供了，它们只会在步骤2和或3中有效。默认的log文件名是group_name + ’_sidechain_fc.log’（对于步骤2）和group_name + ’_large_mk.log’（对于步骤3）。默认的fchk文件名为group_name + ’_sidechain_opt.fchk’，它只在步骤2中有效，当使用Gaussian程序和Seminario方法去获得力常数参数。

如果你使用Gaussian软件和Seminario方法去生成力常数参数，它将使用侧链模型的fchk文件去储存笛卡尔型的Hessian矩阵。如果你使用Gaussian软件和Z矩阵方法去生成力常数参数，它将使用侧链模型的log文件去储存笛卡尔型的Hessian矩阵。当你使用GAMESS-US软件和Seminario方法去生成力常数，侧链模型的log文件去储存笛卡尔型的Hessian矩阵。当前版本不支持GAMESS-US和Z矩阵方法去生成力常数。Gaussian和GAMESS-US软件使用大模型的log文件去储存ESP电荷。