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

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的概率改变Ne0.28的概率扰动原子坐标R0.32的概率在PIMC中扰动量子珠的内部自由度R_Δ^(k)随机决定本次尝试移动trial move的类型。执行尝试移动与接受判断改变NeMetropolis准则提议一个新的电子数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{-β Σ_{k1}^P [ (1/2) ω_P^2 (R^(k1)-R^(k))^T m (R^(k1)-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_solIS的构型蓝色部分呈现出H2FS的构型黄色。这种“既在此处又在彼处”的量子特性正是隧穿效应的体现。质子不再需要像经典粒子那样必须翻越势垒的顶点而是可以通过其波函数的非零概率幅直接“穿透”势垒从而有效降低了反应的活化能垒。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。我们测试了P16和P32两种情况图S6(a)发现自由能剖面差异很小证明P16对于室温下的质子已足够。温度效应验证作为合理性检查我们计算了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等更复杂电催化体系中质子/氢原子量子行为的标准工具之一。它让我们意识到在追求“高大上”的催化剂设计之前回归到反应最基本的粒子——质子——的量子本性或许能发现那些被经典面纱所掩盖的关键细节。