Hydrogen tunneling underlies many chemical and biological processes, from enzymatic reactions to proton transfer in materials. Measuring kinetic isotope effects (KIEs), sensitive indicators of tunneling, requires accurate theoretical models. Existing approaches, such as wave packets or path integrals, are accurate but computationally expensive for large systems. The NEO-MSDFT method, developed earlier, offers an efficient alternative, but until now there was no way to systematically explore the configuration space where tunneling actually occurs. This work fills that gap by introducing V-SHAKE, inspired by a similar approach for electronic conical intersections. It fits into the line of research developing NEO methods for dynamics and kinetics, but stands out by directly sampling the diabatic seam for vibronic states for the first time, rather than just optimizing individual points.
V-SHAKE is based on nuclear-electronic orbital (NEO) theory, where selected nuclei, typically protons or deuterons, are quantized on equal footing with electrons. To describe tunneling, NEO-MSDFT (multistate density functional theory within NEO) is used: two diabatic vibronic states, reactant and product, are defined by localization of the proton density at the donor or acceptor, respectively. The diabatic seam is the subspace of configurations where the energies of these states are equal. V-SHAKE imposes this condition as a rigid constraint in the equations of motion for classical nuclei, using Lagrange multipliers (similar to the RATTLE algorithm). Trajectories were initialized from the minimum energy crossing point (MECP) found earlier and thermostatted at 300 K and 1000 K to ensure broad coverage of the seam. Calculations were performed in the Q-Chem program with the B3LYP functional for electrons and epc17-2 for electron-proton correlation, basis aug-cc-pVDZ for electrons and PB4-D for the proton (for deuteron, exponents were scaled by √2). For each molecule and isotope, at least 1000 (CBO) and 2000 (HEO) configurations were accumulated with a time step of 1 fs.
Analysis of the distribution of vibronic couplings along the diabatic seam revealed significant variation: for CBO and HEO, the coupling varied widely, with configuration energies relative to MECP reaching 40–60 kcal/mol. For HEO, the coupling generally decreased with increasing energy, consistent with an increase in donor-acceptor distance. For CBO, in contrast, a noticeable number of structures with energy close to MECP but with significantly larger coupling were found, indicating greater geometric freedom in the seam. In both cases, the Condon approximation (constancy of vibronic coupling in the crossing region) is violated, but for HEO this does not affect kinetics because structures with small coupling do not contribute to the rate. For CBO, non-Condon effects may be important. The key coordinate modulating the coupling turned out to be the donor-acceptor distance, confirmed both by analysis of the coupling gradient at MECP and by the sliced inverse regression (SIR) method. The reaction coordinate vector g⊥ at MECP mainly involves compression/extension of covalent bonds at the donor and acceptor, while the coupling gradient involves compression of the donor-acceptor distance. These vectors were found to be nearly orthogonal (cosine of the angle close to zero).
This work provides the first direct protocol for comprehensive sampling of diabatic seams in systems with hydrogen tunneling. This allows testing fundamental assumptions of rate theories, such as the Condon approximation and the separation of the reaction coordinate and the donor-acceptor distance. Unlike previous works that either optimized individual MECP points or relied on a priori assumptions, V-SHAKE gives a complete picture of the seam and identifies dominant coordinates. The validity of thermal averaging over the donor-acceptor distance in NEO-GRT and similar theories is confirmed. Furthermore, the SIR method offers a general way to identify important motions without prior knowledge of the system, which could be extended to more complex cases.
The V-SHAKE method can be generalized to sample seams between excited diabatic vibronic states, which will require further methodological developments. It may be combined with NEO approaches for the condensed phase (e.g., NEO-QM/MM) to study tunneling in complex chemical and biological systems. Also promising is the use of more accurate NEO methods, such as NEO-CASSCF or NEO-MRCI, to improve the quantitative reliability of vibronic couplings.
The results are important for calculating rate constants and kinetic isotope effects in chemistry, biology, and materials science. The method could be useful in studying enzymatic reactions, proton-coupled electron transfer, and other processes where hydrogen tunneling plays a key role.
Immediate next steps include applying V-SHAKE to a broader range of molecular systems to test the universality of the conclusions. It is also necessary to integrate the method with more accurate electronic structure methods and extend it to multinuclear tunneling.
In chemical physics, the question of accurate quantitative description of vibronic couplings in the strong tunneling regime remains open, especially for systems with low proton density. This work offers a tool for systematic investigation of this problem and for identifying the limits of applicability of the Condon approximation. Also unsolved is the inclusion of dynamic environmental effects in tunneling processes, which requires further development of sampling methods.
This section is addressed to the author of the paper and their colleagues. Neighbouring works were found by meaning-based search across our archive, not by word overlap. We do not review: we show the field around the work.
What the work does
The paper introduces the V-SHAKE method for systematic sampling of the diabatic seam within NEO-MSDFT, applied to two molecules (CBO and HEO). It is shown that the vibronic coupling varies significantly along the seam, mainly due to changes in the donor–acceptor distance, and that the reaction coordinate is nearly orthogonal to the gradient of the vibronic coupling at the minimum-energy crossing point.
Where its strength is
The V-SHAKE method allows direct exploration of the distribution of vibronic couplings on the diabatic seam, which paves the way for refining theories of tunneling reaction rates without costly dynamical calculations.
Where it could go next
Worth trying
Why this may matter
Such an approach could provide a practical tool for estimating tunneling contributions to kinetics in large systems, where accurate dynamical methods are too expensive, thereby improving predictions of kinetic isotope effects in biology and materials science.
Nearby works in our archive
The archive is sparse around this work: the closest review is only 0.50, against a usual 0.64. It looks like few have drilled nearby.
The vibronic coupling between reactant and product states can change by an order of magnitude when the donor-acceptor distance changes by only 0.1 Å, which is comparable to the thickness of a single atomic layer.