|
QMMM参数-讲解在分子力学（MM）力场中，电子的影响由基于实验数据或基于高级量子化学计算结果分配的经验参数表示。这些对于给定共价结构的基态是有效的，对于系统中原子之间的整体连接性保持不变的基态过程，MM近似通常足够准确。然而，对于连接确实发生变化的过程，如化学反应，或涉及多个电子态的过程，如光化学转换，电子不能再被忽视，至少需要对发生反应的系统部分进行量子力学描述。 模拟溶液或酶中化学反应的一种方法是结合量子力学（QM）和分子力学（MM）。系统的反应部分被量子力学处理，其余部分使用力场建模。GROMACS的当前版本提供了一个与流行的量子化学软件包CP2K的接口。 cp2k联用 QM和MM子系统之间的GROMACS交互是使用Laino等人所述的GEEP方法处理的。这种评估QM和MM子系统之间相互作用的方法是“静电嵌入”方案的变体。QM区域的电子与MM原子之间以及QM核与MM原子间的静电相互作用明确地包含在QM子系统的哈密顿量中：
QM区域中的电子数量是n和原子核M，M是MM区域中带电的原子数量。 右侧的第一项是孤立QM系统的原始电子哈密顿量。 第二项是QM电子和MM原子之间的总静电相互作用。 第三项是QM核与MM原子的总静电相互作用。 使用CP2K/GEEP组合的一个重要优点是，在具有周期性边界条件（PBC）的系统中，它允许评估QM-QM和QM-MM相互作用的力。 为了避免对静电相互作用和LJ的双重考虑，QM原子上的经典MM电荷被清零，QM-QM原子之间的LJ相互作用也被排除在外。 应该指出的是，QM-MM原子之间的LJ相互作用仍然由GROMACS保持并计算。 QM和MM原子之间的键合相互作用在MM水平上由适当的力场项描述。所有由2个QM原子组成的键、含有2或3个QM分子的角和沉降、含有3或4个QM元素的二羟基均不包括在力场评估中。 QM和MM子系统之间的化学键断裂需要在QM计算中加以限制。这是在CP2K内通过添加氢原子来完成QM区的化合价来完成的。该原子上的力仅存在于QM区域，分布在键的两个原子上。 cap原子通常被称为连接原子。在接口内，所有描述的拓扑修改都是在gmx-grompp预处理期间自动执行的。 CP2K 8.1版（或更高版本）应作为libcp2k链接到GROMACS。 有关具体的安装说明，请按照《使用CP2K QM/MM构建支持指南》(https://manual.gromacs.org/current/install-guide/index.html#installing-with-cp2k)进行操作。 使用CP2K接口的QM/MM模拟可以通过设置mdp文件选项来控制，在某些情况下，还可以通过命令行选项为gmx-grompp提供额外的输入文件。 与CP2K的QM/MM模拟相关的所有选项都以前缀开头-qmiqmm-cp2k 设置将触发QM/MM模拟，使用整个系统作为QM部分，并为所有其他选项设置默认参数。 选择系统的QM部分，其名称对应于GROMACS索引文件中的原子组，对应于mdp文件中的qmm-cp2k-qmgroup选项。 典型的QM部分应该由从化学角度来看有趣的原子组成，即反应发生的系统的一部分。 为了使QM部分的计算可行，它应该在空间中尽可能小和紧凑。 DFT模拟通常按QM部分中原子数的三阶进行缩放。这意味着将QM部分中的原子数量增加2倍将使模拟速度减慢8倍。 此外，用户应使用qmm-cp2k-qmcharge选项提供QM子系统的总电荷，并使用qmm-crp2k-qmdersity选项提供自旋状态（多重性）。 QM方法是通过mdp文件中的qmm-cp2k-qmmathod选择的。目前支持以下QM方法： qmm-cp2k-qmethod=PBE-使用PBE泛函和DZVP-MOLOPT基组的DFT。 qmm-cp2k-qmethod=BLYP-使用BLYP泛函和DZVP-MOLOPT基组的DFT。 一旦测试并包含在界面中，该列表将使用新方法进行更新。 此外，还可以使用qmm-CP2K-qmmethod=input的自定义外部CP2K输入文件，并使用gmx-grompp提供带有选项的文件。 外部文件将被合并到模拟的tpr文件中，并受以下限制：-qmi CP2K输入中的RUN_TYPE选项应等于。ENERGY_FORCE 应提供CHARGE选项。 应提供MULTIPILICTY选项。 COORD_FILE_NAME选项应指向pdb文件。 两者都应该存在。CHARGE_EXTENDED TRUECOORD_FILE_FORMAT PDB CP2K输入文件中不允许增量包含（指令）@INCLUDE 在gmx-mdrun模拟过程中，将使用、和生成其他文件。它们包含CP2K输入、CP2K输出和pdb文件，其中包含扩展β场中MM原子的点电荷。 默认情况下，所有与CP2K相关的文件名将通过添加后缀从tpr模拟文件名中推断出来。 为了手动更改它，应使用qmm-cp2k-qmfilenames选项..inp.out.pdb_cp2k cp2k联用局限性 QM/MM接口以两种方式限制了模拟。 首先，在QM区域的模拟过程中，不可能进行拓扑修改。 其次，接口完全忽略了拓扑中的“B”状态参数，使得双拓扑设置成为不可能，例如自由能扰动模拟（自由能实现）。 此外，应该指出的是，QM/MM对系统维里力的贡献没有被考虑在内。 由于QM-MM相互作用，对压力耦合算法的影响大小随着总合力的增加而增加，并可能在NPT系综的模拟中产生伪影。 MiMiC联用 本节描述了与新型QM/MM接口的耦合。计算化学中的多尺度建模（MiMiC）接口将GROMACS与CPMD QM代码相结合。在QM/MM方法中，通常在QM理论水平上处理系统的一小部分（例如可以发生化学反应的酶的活性位点）（因为我们在描述某些过程（如化学反应）时不能忽视电子自由度），而系统的其余部分（蛋白质、溶剂等的剩余部分）由经典力场（MM）描述。MiMiC实现了描述的U.Roethlisberger教授小组开发的QM/MM耦合方案。这种加性方案使用量子哈密顿量内经典系统的静电嵌入。QM/MM总能量计算为子系统贡献的总和：
QM贡献由CPMD计算，而MM部分由GROMACS处理，交叉项由MiMiC接口处理。 交叉项，即同时涉及QM区原子和MM区原子的项，由键合和非键合相互作用组成。 键的相互作用取自用于描述MM部分的力场。每当有化学键穿过QM/MM边界时，必须格外小心，正确处理这种情况。否则，参与切割键的QM原子会留下不饱和电子轨道，导致非物理系统行为。因此，悬空键必须用另一个QM原子封端。CPMD中有两种不同的键上限选项： 氢封端-最简单的方法是用氢原子封端键，限制其相对位置 连接原子伪势——这种策略使用一种特殊的伪势来限制键。这种伪势将代表真实的原子，因此不需要键约束。 与标准力场一样，对EQM/MM可分为范德华贡献和静电贡献。第一个贡献再次来自MM力场。非键合相互作用的第二部分由MiMiC在静电嵌入方法中处理。这为系统的哈密顿量增加了额外的项：
Nmm是MM原子的数量 Nqm是QM原子的数量 rca是MM原子的共价半径。 上述第一项对应于QM和MM区的原子中电子密度ρ(r)之间的阻尼库仑相互作用 由于CPMD使用平面波基组来扩展电子波函数，因此需要阻尼。与局域基组不同，平面波是离域的，这可能会导致所谓的电子溢出问题：带正电的MM原子可能会由于缺乏通常会阻止它的量子力学效应（例如泡利排斥）而人为地使电子云过极化（在完全量子系统中）。 由于计算上述第一项中的积分可能非常昂贵，MiMiC还实现了分层静电嵌入方案，以减轻计算所需的巨大计算工作量 在这个方案中，MM原子根据与QM区域的距离分为两个壳层：短程壳层和长程壳层。对于短程壳层中的MM原子，使用上述方程计算QM/MM相互作用。与此相反，使用QM静电势的多极展开计算了涉及来自长程壳的MM原子的相互作用。 与大多数QM/MM接口不同，MiMiC在合作伙伴代码之间使用松散耦合。 这意味着MiMiC不是将两个代码编译成一个二进制文件，而是为CPMD和GROMACS构建单独的可执行文件。然后，用户将为这两个代码准备输入并同时运行它们。 每个代码都使用单独的MPI进程池运行，并通过MPI客户端-服务器机制传递必要的数据（例如坐标、能量和力）。在MiMiC框架中，CPMD充当服务器，GROMACS成为客户端。 GROMACS version 2019+. Newer major releases may support multiple versions of MiMiC. CPMD version 4.1+. 在使用MiMiC QM/MM支持进行构建后(https://manual.gromacs.org/current/install-guide/index.html#installing-with-mimic)，要运行MiMiC QH/MM模拟，需要： 获取并编译支持MiMiC的CPMD。 用GROMACS进行正常的经典平衡。 在GROMACS中创建一个表示QM原子的索引组。请记住，该组还应包括与QM区域中的原子结合的连接原子，因为它们必须在量子水平上进行处理。 根据以下建议为CPMD和GROMACS准备输入。 将CPMD和GROMACS作为单个批处理作业中的两个独立实例运行。 为了为MiMiC模拟设置mdp文件，需要添加两个选项： integrator=mimic，以在GROMACS中启用mimic工作流。 QMMM grps=<name_of_qm_index_group>表示CPMD将处理的所有原子。 由于CPMD将执行MD集成，因此只有与力计算和输出相关的mdp选项处于活动状态。 设置mdp文件后，可以像往常一样运行grompp。grompp将把所有QM原子的电荷设置为零，以避免库仑相互作用的重复计数。 此外，它将更新非键合排除列表，以排除QM原子之间的LJ相互作用（因为它们将由CPMD描述）。 最后，它将去除QM原子之间的键（如果存在）。我们建议使用gmx-grompp-pp<prepared_topology_file>输出预处理的拓扑文件，因为它将有助于以自动化的方式为CPMD准备输入。 CPMD输入文件 本节将仅描述CPMD中与MiMiC相关的输入-有关DFT相关选项的配置-请参阅CPMD手册(https://www.cpmd.org/)。在为GROMACS准备输入并获得预处理的拓扑文件后，只需运行MiMiC发行版中提供的Python预处理器脚本，即可获得CPMD输入文件的MiMiC相关部分。脚本的用法很简单： prepare-qmmm.py <index_file> <gro_file> <preprocessed_topology_file> <qm_group_name> 请注意，对于MiMiC，力场包含每种原子类型的元素数数据至关重要！如果它没有提供它，预处理器将失败并出现错误： It looks like the forcefield that you are using has no information about the element number. The element number is needed to run QM/MM simulations. 给定所有相关信息，脚本将打印与MiMiC相关的CPMD输入部分。以下是示例输出，其中包含CPMD输入这一部分中可以找到的关键字的简短描述： &MIMIC PATHS 1 <some_absolute_path> BOX 35.77988547402689 35.77988547402689 35.77988547402689 OVERLAPS 3 2 13 1 1 2 14 1 2 2 15 1 3 &END &ATOMS O 1 17.23430225802002 17.76342557295923 18.576007806615877 H 2 18.557110545368047 19.086233860307257 18.727185896598506 17.57445296048094 16.705178943080806 17.06422690678956 &END Suggested QM box size [12.661165036045407, 13.71941166592383, 13.00131573850633] PATHS表示模拟中涉及的MM客户端代码的数量及其各自文件夹的绝对路径。 请记住，此路径必须指向要运行GROMACS的文件夹，否则将导致CPMD中的死锁！下一行包含MM代码的数量（本例中为1），下一行包含各自工作目录的路径 BOX以X Y Z格式表示玻尔中整个模拟框的大小 OVERLAPS-设置GROMACS中将由CPMD处理的原子的数量和ID。格式如下： <code_id><atom_id_in_code><host_code_id><atom_id_inthat_code> CPMD主机代码id始终为id 1。因此，在QM/MM模拟中，GROMACS将具有代码ID 2。 （可选）长程耦合LONG-RANGE COUPLING-使位于距离QM盒一定距离处的原子能够更快地进行多极耦合 （可选）截止距离CUTOFF DISTANCE-下一行包含显式库仑耦合的截止值（如果存在长程耦合，默认为20玻尔） （可选）多极顺序MULTIPOLE ORDER-下一行将包含多极展开被截断的顺序（默认2，最大20）。 CPMD输入文件的&ATOMS部分包含系统中的所有QM原子，并具有默认的CPMD格式。请参阅CPMD手册，根据您的需要进行调整（需要为每种原子种类设置正确的伪电势）。 最后，预处理器建议包含电子密度的QM盒的大小。建议值不是最终值，可能需要用户进一步调整。 MiMiC联用 为了运行模拟，需要在一个作业中运行GROMACS和CPMD。在绝大多数排队系统中，这很容易做到。 例如，在SLURM队列系统中，可以在一个作业中使用两个作业步骤。以下是运行242节点slurm作业的示例作业脚本，为GROMACS分配2个节点，为CPMD分配240个节点（这两个代码都在同一文件夹中启动）： #!/bin/bash -x #SBATCH --nodes=242 #SBATCH --output=mpi-out.%j #SBATCH --error=mpi-err.%j #SBATCH --time=00:25:00 #SBATCH --partition=batch # *** start of job script *** srun -N2 --ntasks-per-node=6 --cpus-per-task=4 -r0 gmx_mpi_d mdrun -deffnm mimic -ntomp 4 & srun -N240 --ntasks-per-node=6 --cpus-per-task=4 -r2 cpmd.x benchmark.inp <path_to_pp_folder> > benchmark-240-4.out & wait |
