大正则路径积分框架：揭示电催化中质子核量子效应的关键作用-编程实验室

1. 项目概述与核心挑战

在电催化领域，尤其是析氢反应（HER）这类涉及质子转移的关键过程中，我们一直面临一个根本性的理论难题：如何精确地描述和模拟在电极/溶液界面处，质子耦合电子转移（PCET）步骤的真实物理图景。传统的计算模拟方法，无论是基于计算氢电极（CHE）模型的静态计算，还是基于第一性原理的约束分子动力学（CMD），都或多或少地做出了妥协。它们要么牺牲了构型采样的充分性，要么采用了固定的电荷或近似恒电位模型，更重要的是，它们几乎无一例外地将质子视为经典的质点来处理。

这带来了一个显著的矛盾：氢是自然界最轻的元素，其核量子效应（NQEs）——简单说就是质子的“波粒二象性”和隧穿行为——在低温下已被证实是主导机制。然而，大量针对室温电催化过程的理论研究却默认忽略了这一效应，部分原因是其计算成本高昂，部分则源于一种潜在的假设：室温下的热涨落足以掩盖量子效应。但近年来越来越多的证据表明，在液态水、生物大分子乃至氧化物材料中，质子的量子行为在室温下依然扮演着重要角色。那么，在电化学界面这个复杂环境中，NQEs是否同样不可忽视？它会对我们理解HER的决速步骤、乃至整个反应路径的选择产生多大影响？

我们的工作正是为了直接回答这个问题。我们开发了一个统一的计算框架，旨在在精确的恒电位（大正则系综）条件下，显式地处理质子的核量子效应。这个框架的核心由三部分组成：1)大正则混合蒙特卡洛（GC-HMC）算法，用于在恒定电化学势下对系统进行采样；2)路径积分蒙特卡洛（PIMC）方法，用于精确描述质子的量子行为；3) 一个为电化学条件量身定制的机器学习势能（MLP）模型——DP-Ne，它将系统的总电子数作为一个额外的自由度，从而能够高效且准确地描述在扩展构型空间[R, Ne]中的势能面。

通过将这个框架应用于Pt(111)表面的HER反应（聚焦于Volmer和Heyrovsky这两个PCET步骤），我们首次在原子尺度上定量揭示了室温下质子隧穿对反应活化能的显著降低作用，并提供了一个全新的物理图像：量子化的质子其波函数可以“弥散”在经典势垒的两侧，这种“不确定性”使得质子能够以隧穿的方式更高效地穿越能垒。这不仅调和了此前理论计算与实验在活化能数值上的长期差异，也深刻影响了我们对HER反应路径（Volmer-Heyrovsky vs. Volmer-Tafel）竞争关系的理解。

2. 方法论核心：构建一个兼顾量子效应与恒电位条件的采样框架

2.1 为何需要大正则（GC）系综？

在电化学模拟中，电极与溶液界面构成一个开放系统，它与外部电子库（电源）处于平衡状态。这意味着系统的总电子数Ne不是一个固定值，而是一个会在其平均值附近波动的统计变量。系统的控制参数是外电子库的电化学势μe（对应于我们施加的电极电位U），而非总电荷。

传统的“固定电荷”或“固定功函数”方法，实际上是对这一开放系统物理图像的近似。固定电荷法忽略了电位涨落；而固定功函数法（如某些GC-DFT实现）则强制每个采样构型的功函数严格等于-μe，这并不符合大正则系综的基本统计规律。在大正则系综中，瞬时功函数W和总电子数Ne是一对共轭的热力学变量，两者都会围绕其平均值涨落，正如等温等压（NPT）系综中体积V和压强P的关系一样。

因此，一个严格的处理必须允许Ne在采样中变化，并确保其系综平均功函数<W>_GC等于-μe。我们的GC-HMC算法正是基于这一原理构建的。

2.2 GC-HMC算法：如何实现恒电位采样？

我们的GC-HMC算法本质上是将系统的总电子数Ne作为一个额外的自由度引入蒙特卡洛（MC）采样。算法的核心流程如下：

初始化：从一个初始构型[R, Ne]开始，其中R是所有原子坐标的集合。
随机选择扰动类型：在每个MC步，根据预设的概率（例如，0.4的概率改变Ne，0.28的概率扰动原子坐标R，0.32的概率在PIMC中扰动量子珠的内部自由度R_Δ^(k)），随机决定本次尝试移动（trial move）的类型。
执行尝试移动与接受判断：
- 改变Ne（Metropolis准则）：提议一个新的电子数Ne'。接受概率为：A(Ne'|Ne) = min[1, exp(-β[E(R, Ne') - E(R, Ne) - μe(Ne' - Ne)])]这里的关键是，能量E是Ne的函数，这要求我们的势能模型必须能处理可变电子数。
- 扰动原子坐标R（混合蒙特卡洛，HMC）：在固定Ne下，使用HMC方法高效地采样原子构型。HMC结合了分子动力学的动力学演化和MC的Metropolis接受/拒绝步骤，能有效克服能垒，在复杂势能面上高效采样。
- 扰动量子珠内部自由度R_Δ^(k)（PIMC）：当考虑NQEs时，使用staging算法更新路径积分中量子珠的相对位置。
估计物理量：对采样轨迹中每个构型计算感兴趣的物理量（如力、能量、功函数等）。
循环与平均：重复步骤2-4，直至达到足够的采样步数，最后对物理量取系综平均。

通过这种方式，我们实现了在恒定μe下，对系统构型空间和电子数空间的联合采样，得到了真正符合大正则分布的微观状态集合。

2.3 路径积分（PI）：如何将质子“量子化”？

经典力学将原子核视为遵循牛顿运动定律的质点。而路径积分表述则将一个量子粒子等价于一个由P个经典“珠子”（beads）通过谐振子耦合连接而成的环状聚合物（ring-polymer）。每个珠子感受的势能是经典势能E(R^(k))的1/P，相邻珠子之间通过频率为ω_P = sqrt(P)/(βħ)的谐振势耦合。

系统的量子配分函数可以映射为这个经典环状聚合物系统的配分函数：Q_qtm(β, Ne) = lim_{P→∞} ∫ dR^(1)...dR^(P) exp{-β Σ_{k=1}^P [ (1/2) ω_P^2 (R^(k+1)-R^(k))^T m (R^(k+1)-R^(k)) + (1/P) E(R^(k), Ne) ] }

其中，R^(k)是第k个珠子的坐标。当P足够大时（室温下对氢原子通常取16或32），这个经典系统的统计平均就能精确再现量子粒子的行为。质子的量子隧穿效应，在这个图像下表现为环状聚合物珠子的空间弥散——珠子可以同时分布在反应坐标上初始态（IS）和末态（FS）对应的区域，即使其质心被约束在过渡态（TS）附近。

在我们的GC-PIHMC框架中，我们将需要量子化处理的质子（如转移中的H+或形成H2的H原子）用环状聚合物表示，其质心坐标R被用作定义反应坐标（RC）的变量。然后，我们在GC-HMC的采样循环中，加入了对这些量子珠内部自由度R_Δ^(k)的PIMC采样。

2.4 势能面引擎：DP-Ne机器学习势能

上述采样框架的可行性，严重依赖于能否快速、准确地计算任意构型[R, Ne]下的能量E和原子受力F。直接使用第一性原理计算（如DFT）进行数百万步的采样是完全不现实的。因此，我们开发了DP-Ne机器学习势能。

DP-Ne的创新之处在于，它将系统的总电子数Ne作为一个明确的输入特征，与原子坐标R一同输入到深度势能（Deep Potential）神经网络中。网络不仅输出总能量E和原子力F，还通过自动微分输出∂E/∂Ne，这个量正好等于该构型瞬时功函数W的负值。这一点至关重要，因为它使我们能够直接验证大正则采样的自洽性：<∂E/∂Ne>_GC = -<W>_GC = μe。

DP-Ne的训练通过DP-GEN（深度势能生成器）流程完成，这是一个“训练-探索-标记”的迭代过程。我们使用ABACUS软件进行第一性原理计算，为不同[R, Ne]构型生成标签数据。最终得到的DP-Ne模型在测试集上表现出色，能量误差低至~1 meV/atom，力误差~70 meV/Å，在保证DFT精度的同时，将计算速度提升了数个数量级，使得大规模GC-PIHMC采样成为可能。

实操心得：DP-Ne训练的关键
数据集的代表性：探索采样必须在目标反应的构型空间和电子数范围内充分进行。对于PCET反应，需要覆盖从反应物到产物整个路径上不同RC值、以及在不同施加电位（对应不同平均Ne）下的构型。
Ne的输入处理：Ne是一个全局标量。在DP框架中，我们将其作为一个额外的环境特征，与每个原子的局部描述符（descriptor）结合后，再输入到拟合网络（fitting net）。确保网络能学习到E对Ne的平滑依赖关系是关键。
收敛判据：在DP-GEN迭代中，我们设定当新探索轨迹中超过85%的构型其模型偏差低于0.10 eV/Å时，认为训练收敛。这个阈值需要根据所研究体系的能量尺度进行调整。

3. 模型构建与计算细节

3.1 界面原子模型

我们选择Pt(111)表面作为模型电极，构建了一个(5×5)的超胞，包含4层Pt原子，共100个Pt原子。在Pt顶位（atop site）吸附了1个单层（ML）的氢原子（H*），以模拟HER反应中常见的表面氢覆盖度。在电极上方，我们放置了两层共36个显式水分子，以模拟溶剂化效应并提供质子源（H3O+）。在水分子上方设置了15 Å的真空层，以消除周期性镜像之间的相互作用。

实现恒电位的关键：补偿电荷板由于我们在GC采样中会改变系统的总电子数Ne，为了在周期性边界条件（PBC）下保持整个超胞的电中性，必须引入补偿电荷。我们采用了一种广泛使用的方案：在真空层中、水层上方约2 Å处放置一个均匀的补偿电荷板。这个电荷板所带电荷与系统净电荷大小相等、符号相反，从而模拟了电极/溶液界面处的双电层（EDL）效应。同时，我们在补偿电荷板上方约7.8 Å处加入了偶极修正，以消除由于平板模型和PBC引入的虚假偶极场。我们在ABACUS代码中实现了这一功能，并验证了其与主流DFT软件Quantum ESPRESSO计算结果的一致性。

3.2 反应坐标（RC）的定义

为了使用热力学积分（TI）方法计算自由能面，我们需要为每个反应定义合适的反应坐标（RC）。

Volmer步骤 (H+_sol + e- + * → H)*：我们定义RC为吸附位点Pt原子与转移质子之间的距离，减去质子给体水分子中O原子与该质子之间的距离。q_Volmer = |r_PtH| - |r_OH|当q_Volmer为较大的负值时，质子更靠近O（水合质子态）；当其为较大的正值时，质子更靠近Pt（吸附氢原子态）。q_Volmer = 0附近对应过渡态。
Heyrovsky步骤 (H+_sol + e- + H→ H2)*：我们定义RC为两个结合形成H2的氢原子之间的距离。q_Heyrovsky = |r_HH*|当距离较大时，对应分离的H*和H+_sol；当距离接近H2分子的平衡键长（~0.74 Å）时，对应产物H2。

在量子情况下，对于Volmer步骤，转移的质子被量子化为环状聚合物；对于Heyrovsky步骤，两个形成H2的氢原子都被量子化。此时，上述距离定义中的原子位置取的是对应环状聚合物所有珠子的质心坐标。

3.3 采样与自由能计算流程

对于每个反应，我们在感兴趣的电位下（例如U = -0.9 V vs. SHE，对应μe = -3.5 eV），选取一系列固定的RC值。在每个RC值下，我们运行独立的GC-PIHMC采样轨迹。

约束采样：使用我们之前开发的约束HMC方法，将系统的RC固定在目标值q = s。此时，系统的其他自由度（原子坐标、电子数、量子珠内部坐标）被充分采样。
计算平均力：在约束采样中，我们记录沿RC方向的平均力<dF/dq>_cond。根据TI原理，自由能F(q)对q的导数正是这个约束平均力。
数值积分：在所有离散的RC点计算得到平均力后，通过数值积分（如辛普森法则）得到整个反应路径的自由能剖面ΔF(q)。F(q2) - F(q1) = ∫_{q1}^{q2} < (dF/dq)_estm >_cond dq
误差评估：为了获得可靠的统计结果，我们对每个RC点进行了多次（Volmer 15次，Heyrovsky 10次）独立的采样轨迹（每次30万步），并计算平均力的标准误差作为误差棒。对于Heyrovsky反应在H2形成的关键区域（q从0.85到0.95 Å），我们甚至将单次采样步数增加到了100万步以确保收敛。

4. 结果分析：核量子效应如何重塑HER反应图像？

4.1 恒电位采样的自洽性验证

首先，我们验证了GC采样框架的正确性。图3(a)展示了在固定RC和电位下，系统瞬时功函数W随MC步数的波动。可以看到，W在一个平均值附近涨落，其系综平均<W>_GC精确等于我们设定的控制参数-μe（-3.5 eV）。涨落范围约为±0.5 eV，这与之前基于电位计的恒电位MD方法的结果一致。同时，系统的总电子数Ne也呈现正态分布（图3(b)），其平均值<Ne>_GC对应于该电位下系统的平均电荷态。这些结果完美符合大正则系综的理论预期。

随着反应从初始态（IS）向末态（FS）进行，系统的平均电子数<Ne>_GC逐渐增加（图3(c)）。这直观地反映了还原反应的本质：在恒定还原电位驱动下，电极需要从外电路持续获取电子来推动反应进行。电位越负（还原性越强），系统所需的平均电子数就越多，这与物理直觉完全吻合。

4.2 NQEs对活化能的定量影响

图3(d)展示了在U = -0.9 V vs. SHE下，Volmer和Heyrovsky步骤的经典与量子自由能剖面对比。结果令人印象深刻：

Volmer步骤：经典计算的活化能约为0.50 eV，而考虑NQEs后，活化能降至约0.37 eV，降低了0.13 eV。
Heyrovsky步骤：经典活化能约为0.46 eV，量子活化能约为0.37 eV，降低了0.09 eV。

根据阿伦尼乌斯公式r ∝ exp(-Ea/k_B T)，在室温300 K下，0.1 eV的活化能降低意味着反应速率将提升50至100倍。这是一个不可忽视的效应！它强烈暗示，过去许多基于经典质子近似计算得到的HER反应速率，可能存在严重的低估。

深度解读：0.1 eV的差异是否显著？有些读者可能会质疑：DFT计算本身就有~0.1-0.2 eV的不确定性，这个0.1 eV的NQEs效应是否有物理意义？我们的回答是：有，且至关重要。原因在于：
内部一致性：我们比较的是在同一套计算框架、相同DFT泛函和设置下，仅切换“经典质子”与“量子质子”模型得到的结果。DFT的系统误差在两者中是相同的，因此差值反映的是真实的物理效应。
理论与实验的矛盾：长期以来，针对Pt表面HER的PCET步骤，诸多经典理论计算预测的活化能普遍高于实验测量值。我们的工作指出，忽略NQEs可能是导致这一偏差的关键因素之一。考虑量子效应后，理论预测更接近实验，这增强了计算模型的预测能力。

4.3 质子隧穿的直观图像：环状聚合物的“弥散”

为什么量子化会导致活化能降低？图3(e)和(f)提供了直观的答案。我们分析了在过渡态（TS）附近RC值时，量子珠的RC值分布。

在经典图像中，质子在TS处被约束在一个确定的位置（图3(e)中的红色虚线）。然而，在PIMC采样中，一个量子质子被表示为16个珠子的环状聚合物。当我们统计每个珠子构型对应的RC值（基于珠子质心定义）时，发现其分布是弥散的（图3(e)中的直方图）。这意味着，即使在质心被约束在TS附近时，仍有一部分珠子分布在对应于IS（反应物）和FS（产物）的RC区域。

图3(f)更生动地展示了这一点。我们固定所有经典原子和量子质子质心的位置，仅对量子珠的内部自由度进行PIMC采样，然后随机选取一个构型进行可视化。可以看到，对于Heyrovsky反应的TS，部分量子珠呈现出H* + H+_sol（IS）的构型（蓝色），部分呈现出H2（FS）的构型（黄色）。这种“既在此处，又在彼处”的量子特性，正是隧穿效应的体现。质子不再需要像经典粒子那样必须翻越势垒的顶点，而是可以通过其波函数的非零概率幅直接“穿透”势垒，从而有效降低了反应的活化能垒。

4.4 对HER反应机理的启示：路径选择的转变

HER生成H2有两条主要路径：Volmer步骤后接Tafel步骤（化学脱附），或接Heyrovsky步骤（电化学脱附）。我们的计算揭示了NQEs对这两条路径竞争关系的深刻影响（图4）。

Tafel步骤不受NQEs影响：我们对Tafel步骤（H* + H* → H2）也进行了计算，发现其活化能无论在经典还是量子情况下都约为0.54 eV，且不随电位变化。这是因为Tafel步骤是表面化学重组过程，不涉及质子的长程转移，量子隧穿效应不显著。
Heyrovsky步骤因隧穿而增强：如上所述，Heyrovsky步骤的活化能因NQEs降低了约0.1 eV。
机理转变：在经典图像下，当电位U ≥ -0.9 V时，Volmer-Tafel路径（活化能更高但不受电位影响）占主导。然而，考虑NQEs后，Heyrovsky步骤的活化能曲线整体下移，使得Volmer-Heyrovsky路径开始占优的电位区间向更正的方向移动了约0.5 V。这意味着在更小的过电位下，电化学脱附路径就可能压倒化学脱附路径。

这一发现与一些实验-理论联合分析的结论相符，也为解决Kronberg等人工作中的困惑提供了线索：他们仅模拟了Volmer和Tafel步骤，发现无法复现实验观测到的Tafel斜率。我们的结果表明，很可能是因为在考虑NQEs后，Volmer-Heyrovsky路径在更宽的电位范围内成为了主导路径，从而改变了整个反应的动力学特征。

5. 技术实现要点与避坑指南

5.1 约束GC-PIHMC算法的实现细节

坐标变换与雅可比行列式：在约束采样中，将反应坐标q从广义坐标中分离出来时，会引入一个雅可比行列式J(q, q_trans)。在计算平均力估计量(dF/dq)_estm时（见主文方法部分公式17-19），必须包含-k_B T ∂/∂q ln J项。对于简单的距离或距离差型RC，这个项有解析表达式。忽略它将导致自由能计算出现系统误差。
采样概率分配：在GC-PIHMC中，需要合理设置改变Ne、移动原子坐标质心R、以及移动量子珠内部坐标R_Δ^(k)这三类尝试移动的概率。我们的经验是，给Ne的变化分配较高的概率（如40%），因为电子数的变化通常能更有效地探索构型空间并加速平衡。具体比例需要通过测试不同体系的采样效率来调整。
刚性墙约束：在界面模型中，为了防止水分子逃逸到真空区，或吸附氢原子迁移到水层中干扰反应，我们在相应位置设置了“刚性墙”势。这是一种简化的处理，但能有效保持界面结构的稳定性，对于获得合理的自由能剖面至关重要。

5.2 DP-Ne模型训练与使用的注意事项

初始数据生成：启动DP-GEN迭代需要初始数据集。一个有效的策略是：先在不考虑NQEs和恒电位（即固定电荷）的条件下，用常规DP模型对反应路径进行初步采样，获得一批构型。然后，对这些构型在几个不同的Ne值（对应不同的电位）下进行DFT单点计算，构建初始的[R, Ne]数据集。
Ne的归一化：输入神经网络的Ne需要进行适当的归一化处理，例如减去体系的中性态电子数，并除以一个特征尺度，使其与其他描述符处于相近的数量级。
验证∂E/∂Ne的准确性：DP-Ne的关键输出之一是∂E/∂Ne。必须用有限的差分法（ΔE/ΔNe）在测试集上验证其准确性。如图S3所示，两者应高度一致。
注意PIMC中的Ne一致性：一个容易被忽略但至关重要的问题是，在路径积分中，所有珠子必须共享同一个Ne值。这是因为在量子配分函数的推导中，不同电子数的电子基态是正交的（公式S2）。如果像某些GC-DFT-PI尝试那样，允许每个珠子有不同的Ne以保持各自功函数恒定，将违反量子统计力学的基本原理，导致错误的结果。

5.3 自由能计算的收敛性判断

平均力的收敛：自由能剖面的准确性直接取决于每个RC点平均力计算的收敛性。除了观察平均力随MC步数的波动是否平稳（图S7），还必须进行多次独立重复计算，用标准误差来评估统计不确定性。对于自由能变化剧烈的区域（如Heyrovsky反应中H-H键形成时），需要更长的采样步数。
珠子数P的收敛性测试：路径积分的精度依赖于珠子数P。我们测试了P=16和P=32两种情况（图S6(a)），发现自由能剖面差异很小，证明P=16对于室温下的质子已足够。
温度效应验证：作为合理性检查，我们计算了250 K下的自由能剖面（图S6(b)(c)）。结果显示，随着温度降低，NQEs导致的活化能降低效应更加显著，这与理论预期一致，也侧面印证了方法在物理上的正确性。

6. 框架的局限性与未来展望

尽管我们的GC-PIHMC框架结合DP-Ne在电化学NQEs模拟上迈出了重要一步，但它仍存在一些局限，这也是未来可以拓展的方向：

热力学而非动力学信息：目前我们通过热力学积分得到的是自由能剖面，即反应的平衡态热力学性质。要获得真实的动力学速率常数，需要结合过渡态理论或进行量子动力学模拟。框架本身可以扩展，例如结合环状聚合物过渡态理论（RPMD）来估算量子速率。
绝热近似：当前工作采用了绝热近似，即假设电子始终处于基态。对于某些非绝热效应显著的PCET反应（如涉及电荷转移与质子运动强烈耦合的情况），需要考虑电子-质子非绝热耦合。Hammes-Schiffer组发展的vibronic耦合理论可以作为一个潜在的整合方向。
恒定pH条件：目前我们只控制了电位（μe）。一个更完整的电化学模拟还需要控制pH值，即允许溶液中的质子数N_H+也作为变量进行GC采样。这可以在现有框架中增加一个独立的“质子数变化”MC模块来实现，如图1(a)所示，具有良好的可扩展性。
溶剂模型：我们使用了显式水分子模型，这虽然准确但计算量大。对于更大的体系或更长的模拟时间，发展结合显式-隐式溶剂模型或更高效的机器学习溶剂势将是必要的。

7. 总结与个人体会

回顾整个工作，从方法论构建到具体应用，我最大的体会是：处理复杂问题需要“分而治之”的思维，但最终的突破往往来自于“合而为一”的整合。核量子效应、恒电位条件、充分构型采样、第一性原理精度，每一个都是电化学模拟中的硬骨头。我们并没有发明一个全新的理论去解决其中一个，而是巧妙地将成熟的路径积分、大正则蒙特卡洛、约束采样和机器学习势能这几个强大的工具编织在一起，形成了一个既严谨又高效的计算框架。

在实际操作中，有几个“坑”是值得后来者特别注意的：

切勿混淆控制变量：在大正则系综中，μe是控制参数，Ne和W是涨落的共轭变量。任何试图在采样中固定瞬时W的做法（如某些GC-DFT实现），从统计力学基本原理上就是错误的。我们的GC-HMC算法严格遵循了系综理论。
量子与经典的边界要清晰：在PIMC中，哪些粒子需要量子化（通常是氢），哪些可以当作经典粒子（如重金属原子），需要根据其德布罗意波长和温度来判断。对于HER中的质子，量子化是必须的；对于Pt原子，在室温下其量子效应通常可以忽略。
MLP是使能技术，但非黑箱：DP-Ne的强大毋庸置疑，但它并非万能。训练数据的质量决定了它的上限。必须确保采样数据覆盖了反应路径、电位范围以及可能出现的异常构型（如质子转移中间态）。主动学习（DP-GEN）流程是构建高质量数据集的利器，但初始种子的选择和探索策略的设置需要结合化学直觉。

这项工作的价值不仅在于给出了HER中NQEs的定量数据，更在于它提供了一套可移植、可扩展的方法论。我相信，这套结合了机器学习加速的大正则路径积分框架，将成为未来研究ORR、CO2RR、NRR等更复杂电催化体系中质子/氢原子量子行为的标准工具之一。它让我们意识到，在追求“高大上”的催化剂设计之前，回归到反应最基本的粒子——质子——的量子本性，或许能发现那些被经典面纱所掩盖的关键细节。