Permutationally invariant polynomial potential energy surfaces for tropolone and H and D atom tunneling dynamics

Targeted Transferable Machine-Learned Potential for Linear Alkanes Trained on C14H30 and Tested for C4H10 to C30H62

Given the great importance of linear alkanes in fundamental and applied research, an accurate machine-learned potential (MLP) would be a major advance in computational modeling of these hydrocarbons. Recently, we reported a novel, many-body permutationally invariant model that was trained specifically for the 44-atom hydrocarbon C14H30 on roughly 250,000 B3LYP energies (Qu, C.; Houston, P. L.; Allison, T.; Schneider, B. I.; Bowman, J. M. J. Chem. Theory Comput. 2024, 20, 9339-9353). Here, we demonstrate the accuracy of the transferability of this potential for linear alkanes ranging from butane C4H10 up to C30H62. Unlike other approaches for transferability that aim for universal applicability, the present approach is targeted for linear alkanes. The mean absolute error (MAE) for energy ranges from 0.26 kcal/mol for butane and rises to 0.73 kcal/mol for C30H62 over the energy range up to 80 kcal/mol for butane and 600 kcal/mol for C30H62. These values are unprecedented for transferable potentials and indicate the high performance of a targeted transferable potential. The conformational barriers are shown to be in excellent agreement with high-level ab initio calculations for pentane, the largest alkane for which such calculations have been reported. Vibrational power spectra of C30H62 from molecular dynamics calculations are presented and briefly discussed. Finally, the evaluation time for the potential is shown to vary linearly with the number of atoms.

MOLPIPx: An end-to-end differentiable package for permutationally invariant polynomials in Python and Rust

In this work, we present MOLPIPx, a versatile library designed to seamlessly integrate permutationally invariant polynomials with modern machine learning frameworks, enabling the efficient development of linear models, neural networks, and Gaussian process models. These methodologies are widely employed for parameterizing potential energy surfaces across diverse molecular systems. MOLPIPx leverages two powerful automatic differentiation engines-JAX and EnzymeAD-Rust-to facilitate the efficient computation of energy gradients and higher-order derivatives, which are essential for tasks such as force field development and dynamic simulations. MOLPIPx is available at https://github.com/ChemAI-Lab/molpipx.

New Algorithms to Generate Permutationally Invariant Polynomials and Fundamental Invariants for Potential Energy Surface Fitting

Symmetric functions, such as Permutationally Invariant Polynomials (PIPs) and Fundamental Invariants (FIs), are effective and concise descriptors for incorporating permutation symmetry into neural network (NN) potential energy surface (PES) fitting. The traditional algorithm for generating such symmetric polynomials has a factorial time complexity of N!, where N is the number of identical atoms, posing a significant challenge to applying symmetric polynomials as descriptors of NN PESs for larger systems, particularly with more than 10 atoms. Herein, we report a new algorithm which has only linear time complexity for identical atoms. It can tremendously accelerate generation process of symmetric polynomials for molecular systems. The proposed algorithm is based on graph connectivity analysis following the action of the generation set of molecular permutational group. For instance, in the case of calculating the invariant polynomials for a 15-atom molecule, such as tropolone, our algorithm is approximately 2 million times faster than the previous method. The efficiency of the new algorithm can be further enhanced with increasing molecular size and number of identical atoms, making the FI-NN approach feasible for systems with over 10 atoms and high symmetry demands.

DFT-Based Permutationally Invariant Polynomial Potentials Capture the Twists and Turns of C14H30

Hydrocarbons are ubiquitous as fuels, solvents, lubricants, and as the principal components of plastics and fibers, yet our ability to predict their dynamical properties is limited to force-field mechanics. Here, we report two machine-learned potential energy surfaces (PESs) for the linear 44-atom hydrocarbon C14H30 using an extensive data set of roughly 250,000 density functional theory (DFT) (B3LYP) energies for a large variety of configurations, obtained using MM3 direct-dynamics calculations at 500, 1000, and 2500 K. The surfaces, based on Permutationally Invariant Polynomials (PIPs) and using both a many-body expansion approach and a fragmented-basis approach, produce precise fits for energies and forces and also produce excellent out-of-sample agreement with direct DFT calculations for torsional and dihedral angle potentials. Going beyond precision, the PESs are used in molecular dynamics calculations that demonstrate the robustness of the PESs for a large range of conformations. The many-body PIPs PES, although more compute intensive than the fragmented-basis one, is directly transferable for other linear hydrocarbons.

Fourier transform microwave spectroscopy of the 13C- and 18O-substituted tropolone. Proton tunneling effect for the isotopic species with the asymmetric potential wells

Fourier-transform microwave spectroscopy has been applied for the 13C/18O-substituted tropolone to observe tunneling-rotation transitions as well as pure rotational transitions. The tunneling-rotation transitions were observed between the 13C-4 and -6 forms as well as between 13C-3 and -7, between 13C-1 and -2, and between 18O-8 and -9 (we denote these tunneling pairs as 13C-46, etc., below) although they have an asymmetric tunneling potential due to the difference in the zero point energy (ZPE). From the observed tunneling splittings ΔEij (0.9800-1.6824 cm-1), the differences in ZPE Δij for the 13C-46, -37, -12, and 18O-89 species are derived to be -0.1104, 0.5652, -1.3682, and 1.3897 cm-1 to agree well with the DFT calculation. The state mixing ratio of the tunneling states decreases drastically from (44%:56%) to (8.7%:91.3%) for 13C-46 and 18O-89 with an increase in the asymmetry Δij of the tunneling potential function. The observed tunneling-rotation interaction constants Fij decrease from 16.001 to 9.224 cm-1 as the differences in ZPE Δij increase, while the diagonal tunneling-rotation interaction constants Fu increase from 1.767 to 13.70 cm-1, explained well by the mixing ratio of the tunneling states.

Machine Learning of Reactive Potentials

In the past two decades, machine learning potentials (MLPs) have driven significant developments in chemical, biological, and material sciences. The construction and training of MLPs enable fast and accurate simulations and analysis of thermodynamic and kinetic properties. This review focuses on the application of MLPs to reaction systems with consideration of bond breaking and formation. We review the development of MLP models, primarily with neural network and kernel-based algorithms, and recent applications of reactive MLPs (RMLPs) to systems at different scales. We show how RMLPs are constructed, how they speed up the calculation of reactive dynamics, and how they facilitate the study of reaction trajectories, reaction rates, free energy calculations, and many other calculations. Different data sampling strategies applied in building RMLPs are also discussed with a focus on how to collect structures for rare events and how to further improve their performance with active learning.

Formic Acid-Ammonia Heterodimer: A New Δ-Machine Learning CCSD(T)-Level Potential Energy Surface Allows Investigation of the Double Proton Transfer

The formic acid-ammonia dimer is an important example of a hydrogen-bonded complex in which a double proton transfer can occur. Its microwave spectrum has recently been reported and rotational constants and quadrupole coupling constants were determined. Calculated estimates of the double-well barrier and the internal barriers to rotation were also reported. Here, we report a full-dimensional potential energy surface (PES) for this complex, using two closely related Δ-machine learning methods to bring it to the CCSD(T) level of accuracy. The PES dissociates smoothly and accurately. Using a 2d quantum model the ground vibrational-state tunneling splitting is estimated to be less than 10-4 cm-1. The dipole moment along the intrinsic reaction coordinate is calculated along with a Mullikan charge analysis and supports the mildly ionic character of the minimum and strongly ionic character at the double-well barrier.

Evaluation of Rate Coefficients in the Gas Phase Using Machine-Learned Potentials

We assess the capability of machine-learned potentials to compute rate coefficients by training a neural network (NN) model and applying it to describe the chemical landscape on the C5H5 potential energy surface, which is relevant to molecular weight growth in combustion and interstellar media. We coupled the resulting NN with an automated kinetics workflow code, KinBot, to perform all necessary calculations to compute the rate coefficients. The NN is benchmarked exhaustively by evaluating its performance at the various stages of the kinetics calculations: from the electronic energy through the computation of zero point energy, barrier heights, entropic contributions, the portion of the PES explored, and finally the overall rate coefficients as formulated by transition state theory.

Accurate fundamental invariant-neural network representation of ab initio potential energy surfaces

Highly accurate potential energy surfaces are critically important for chemical reaction dynamics. The large number of degrees of freedom and the intricate symmetry adaption pose a big challenge to accurately representing potential energy surfaces (PESs) for polyatomic reactions. Recently, our group has made substantial progress in this direction by developing the fundamental invariant-neural network (FI-NN) approach. Here, we review these advances, demonstrating that the FI-NN approach can represent highly accurate, global, full-dimensional PESs for reactive systems with even more than 10 atoms. These multi-channel reactions typically involve many intermediates, transition states, and products. The complexity and ruggedness of this potential energy landscape present even greater challenges for full-dimensional PES representation. These PESs exhibit a high level of complexity, molecular size, and accuracy of fit. Dynamics simulations based on these PESs have unveiled intriguing and novel reaction mechanisms, providing deep insights into the intricate dynamics involved in combustion, atmospheric, and organic chemistry.

Proton-Tunneling Dynamics along Low-Barrier Hydrogen Bonds: A Full-Dimensional Instanton Study of 6-Hydroxy-2-formylfulvene

Understanding the dynamics of proton transfer along low-barrier hydrogen bonds remains an outstanding challenge of great fundamental and practical interest, reflecting the central role of quantum effects in reactions of chemical and biological importance. Here, we combine ab initio calculations with the semiclassical ring-polymer instanton method to investigate tunneling processes on the ground electronic state of 6-hydroxy-2-formylfulvene (HFF), a prototypical neutral molecule supporting low-barrier hydrogen-bonding. The results emerging from a full-dimensional ab initio instanton analysis reveal that the tunneling path does not pass through the instantaneous transition-state geometry. Instead, the tunneling process involves a multidimensional reaction coordinate with concerted reorganization of the heavy-atom skeletal framework to substantially reduce the donor-acceptor distance and drive the ensuing intramolecular proton-transfer event. The predicted tunneling-induced splittings for HFF isotopologues are in good agreement with experimental findings, leading to percentage deviations of only 20-40%. Our full-dimensional results allow us to characterize vibrational contributions along the tunneling path, highlighting the intrinsically multidimensional nature of the attendant hydron-migration dynamics.

PESPIP: Software to fit complex molecular and many-body potential energy surfaces with permutationally invariant polynomials

We wish to describe a potential energy surface by using a basis of permutationally invariant polynomials whose coefficients will be determined by numerical regression so as to smoothly fit a dataset of electronic energies as well as, perhaps, gradients. The polynomials will be powers of transformed internuclear distances, usually either Morse variables, exp(-r_i,j/λ), where λ is a constant range hyperparameter, or reciprocals of the distances, 1/r_i,j. The question we address is how to create the most efficient basis, including (a) which polynomials to keep or discard, (b) how many polynomials will be needed, (c) how to make sure the polynomials correctly reproduce the zero interaction at a large distance, (d) how to ensure special symmetries, and (e) how to calculate gradients efficiently. This article discusses how these questions can be answered by using a set of programs to choose and manipulate the polynomials as well as to write efficient Fortran programs for the calculation of energies and gradients. A user-friendly interface for access to monomial symmetrization approach results is also described. The software for these programs is now publicly available.

Wavepacket dynamical study of H-atom tunneling in catecholate monoanion: the role of intermode couplings and energy flow

We present a study of H-atom tunneling in catecholate monoanion through wavepacket dynamical simulations. In our earlier study of this symmetrical double-well system [Phys. Chem. Chem. Phys., 2022, 24, 10887], a limited number of transition state modes were identified as being important for the tunneling process. These include the imaginary frequency mode Q₁, the CO scissor mode Q₁₀, and the OHO bending mode Q₂₉. In this work, starting from non-stationary initial states prepared with excitations in these modes, we have carried out wavepacket dynamics in two and three dimensional spaces. We analyse the dynamical effects of the intermode couplings, in particular the role of energy flow between the studied modes on H-atom tunneling. We find that while Q₁₀ strongly modulates the donor-acceptor distance, it does not exchange energy with Q₁. However, excitation in Q₂₉ or Q₁ does lead to rapid energy exchange between these modes, which modifies the tunneling rate at early times.

Δ-Machine Learned Potential Energy Surfaces and Force Fields

There has been great progress in developing machine-learned potential energy surfaces (PESs) for molecules and clusters with more than 10 atoms. Unfortunately, this number of atoms generally limits the level of electronic structure theory to less than the "gold standard" CCSD(T) level. Indeed, for the well-known MD17 dataset for molecules with 9-20 atoms, all of the energies and forces were obtained with DFT calculations (PBE). This Perspective is focused on a Δ-machine learning method that we recently proposed and applied to bring DFT-based PESs to close to CCSD(T) accuracy. This is demonstrated for hydronium, N-methylacetamide, acetyl acetone, and ethanol. For 15-atom tropolone, it appears that special approaches (e.g., molecular tailoring, local CCSD(T)) are needed to obtain the CCSD(T) energies. A new aspect of this approach is the extension of Δ-machine learning to force fields. The approach is based on many-body corrections to polarizable force field potentials. This is examined in detail using the TTM2.1 water potential. The corrections make use of our recent CCSD(T) datasets for 2-b, 3-b, and 4-b interactions for water. These datasets were used to develop a new fully ab initio potential for water, termed q-AQUA.

Semiclassical and VSCF/VCI Calculations of the Vibrational Energies of trans- and gauche-Ethanol Using a CCSD(T) Potential Energy Surface

A recent full-dimensional Δ-Machine learning potential energy surface (PES) for ethanol is employed in semiclassical and vibrational self-consistent field (VSCF) and virtual-state configuration interaction (VCI) calculations, using MULTIMODE, to determine the anharmonic vibrational frequencies of vibration for both the trans and gauche conformers of ethanol. Both semiclassical and VSCF/VCI energies agree well with the experimental data. We find significant mixing between the VSCF basis states due to Fermi resonances between bending and stretching modes. The same effects are also accurately described by the full-dimensional semiclassical calculations. These are the first high-level anharmonic calculations using a PES, in particular a "gold-standard" CCSD(T) one.

Quantum Calculations on a New CCSD(T) Machine-Learned Potential Energy Surface Reveal the Leaky Nature of Gas-Phase Trans and Gauche Ethanol Conformers

Ethanol is a molecule of fundamental interest in combustion, astrochemistry, and condensed phase as a solvent. It is characterized by two methyl rotors and trans (anti) and gauche conformers, which are known to be very close in energy. Here we show that based on rigorous quantum calculations of the vibrational zero-point state, using a new ab initio potential energy surface (PES), the ground state resembles the trans conformer, but substantial delocalization to the gauche conformer is present. This explains experimental issues about identification and isolation of the two conformers. This “leak” effect is partially quenched when deuterating the OH group, which further demonstrates the need for a quantum mechanical approach. Diffusion Monte Carlo and full-dimensional semiclassical dynamics calculations are employed. The new PES is obtained by means of a Δ-machine learning approach starting from a pre-existing low level density functional theory surface. This surface is brought to the CCSD(T) level of theory using a relatively small number of ab initio CCSD(T) energies. Agreement between the corrected PES and direct ab initio results for standard tests is excellent. One- and two-dimensional discrete variable representation calculations focusing on the trans–gauche torsional motion are also reported, in reasonable agreement with experiment.

The MD17 datasets from the perspective of datasets for gas-phase “small” molecule potentials

There has been great progress in developing methods for machine-learned potential energy surfaces. There have also been important assessments of these methods by comparing so-called learning curves on datasets of electronic energies and forces, notably the MD17 database. The dataset for each molecule in this database generally consists of tens of thousands of energies and forces obtained from DFT direct dynamics at 500 K. We contrast the datasets from this database for three “small” molecules, ethanol, malonaldehyde, and glycine, with datasets we have generated with specific targets for the potential energy surfaces (PESs) in mind: a rigorous calculation of the zero-point energy and wavefunction, the tunneling splitting in malonaldehyde, and, in the case of glycine, a description of all eight low-lying conformers. We found that the MD17 datasets are too limited for these targets. We also examine recent datasets for several PESs that describe small-molecule but complex chemical reactions. Finally, we introduce a new database, “QM-22,” which contains datasets of molecules ranging from 4 to 15 atoms that extend to high energies and a large span of configurations.

Quantitative molecular simulations

All-atom simulations can provide molecular-level insights into the dynamics of gas-phase, condensed-phase and surface processes. One important requirement is a sufficiently realistic and detailed description of the underlying intermolecular interactions. The present perspective provides an overview of the present status of quantitative atomistic simulations from colleagues' and our own efforts for gas- and solution-phase processes and for the dynamics on surfaces. Particular attention is paid to direct comparison with experiment. An outlook discusses present challenges and future extensions to bring such dynamics simulations even closer to reality.

Multidimensional H-atom tunneling in the catecholate monoanion

We present the catecholate monoanion as a new model system for the study of multidimensional tunneling. It has a symmetrical O-H double-well structure, and the H atom motion between the two wells is coupled to both low and high frequency modes with different strengths. With a view to studying mode-specific tunneling in the catecholate monoanion, we have developed a full (33) dimensional potential energy surface in transition state (TS) normal modes using a Distributed Gaussian Empirical Valence Bond (DGEVB) based approach. We have computed eigenstates in different subspaces using both unrelaxed and relaxed potentials based on the DGEVB model. With unrelaxed potentials, we present results up to 7D subspaces that include the imaginary frequency mode and six modes coupled to it. With relaxed potentials, we focus on the two most important coupling modes. The structures of the ground and vibrationally excited eigenstates are discussed for both approaches and mode-specific tunneling splitting and their trends are presented.

Enabling Rapid and Accurate Construction of CCSD(T)-Level Potential Energy Surface of Large Molecules Using Molecular Tailoring Approach

The construction of a potential energy surface (PES) of even a medium-sized molecule employing correlated theory, such as CCSD(T), is arduous due to the high computational cost involved. The present study reports the possibility of efficiently constructing such a PES of molecules containing up to 15 atoms and 550 basis functions by employing the fragment-based molecular tailoring approach (MTA) on off-the-shelf hardware. The MTA energies at the CCSD(T)/aug-cc-pVTZ level for several geometries of three test molecules, viz., acetylacetone, N-methylacetamide, and tropolone, are reported. These energies are in excellent agreement with their full calculation counterparts with a time advantage factor of 3-5. The energy barrier from the ground to transition state is also accurately captured. Further, we demonstrate the accuracy and efficiency of MTA for estimating the energy gradients at the CCSD(T) level. As a further application of our MTA methodology, the energies of acetylacetone at ∼430 geometries are computed at the CCSD(T)/aug-cc-pVTZ level and used for generating a Δ-machine learning (Δ-ML) PES. This leads to the H-transfer barrier of 3.02 kcal/mol, well in agreement with the benchmarked barrier of 3.19 kcal/mol. The fidelity of this Δ-ML PES is examined by geometry optimization and normal mode frequency calculations of global minima and saddle point geometries. We trust that the present work is a major development for the rapid and accurate construction of PES at the CCSD(T) level for molecules containing up to 20 atoms and 600 basis functions using off-the-shelf hardware.

Transfer learned potential energy surfaces: accurate anharmonic vibrational dynamics and dissociation energies for the formic acid monomer and dimer

The vibrational dynamics of the formic acid monomer (FAM) and dimer (FAD) is investigated from machine-learned potential energy surfaces at the MP2 (PES_MP2) and transfer-learned (PES_TL) to the CCSD(T) levels of theory. The normal mode (MAEs of 17.6 and 25.1 cm⁻¹) and second order vibrational perturbation theory (VPT2, MAEs of 6.7 and 17.1 cm⁻¹) frequencies from PES_TL for all modes below 2000 cm⁻¹ for FAM and FAD agree favourably with experiment. For the OH stretch mode the experimental frequencies are overestimated by more than 150 cm⁻¹ for both FAM and FAD from normal mode calculations. Conversely, VPT2 calculations on PES_TL for FAM reproduce the experimental OH frequency to within 22 cm⁻¹. For FAD the VPT2 calculations find the high-frequency OH stretch at 3011 cm⁻¹, compared with an experimentally reported, broad (∼100 cm⁻¹) absorption band with center frequency estimated at ∼3050 cm⁻¹. In agreement with earlier reports, MD simulations at higher temperature shift the position of the OH-stretch in FAM to the red, consistent with improved sampling of the anharmonic regions of the PES. However, for FAD the OH-stretch shifts to the blue and for temperatures higher than 1000 K the dimer partly or fully dissociates using PES_TL. Including zero-point energy corrections from diffusion Monte Carlo simulations for FAM and FAD and corrections due to basis set superposition and completeness errors yields a dissociation energy of D₀ = −14.23 ± 0.08 kcal mol⁻¹ compared with an experimentally determined value of −14.22 ± 0.12 kcal mol⁻¹.

Neural network based PESs are constructed for formic acid monomer and dimer at the MP2 and transfer learned to the CCSD(T) level of theory. The PESs are used to study the vibrational dynamics and dissociation energy of the molecules.

A CCSD(T)-Based 4-Body Potential for Water

High-level, ab initio calculations find that the 4-body (4-b) interaction is needed to account for near-100% of the total interaction energy for water clusters as large as the 21-mer. Motivated by this, we report a permutationally invariant polynomial potential energy surface (PES) for the 4-body interaction. This machine-learned PES is a fit to 2119 symmetry-unique, CCSD(T)-F12a/haTZ 4-b interaction energies. Configurations for these come from tetramer direct-dynamics calculations, fragments from an MD water simulation at 300 K, and tetramer fragments in a variety of water clusters. The PIP basis is purified to ensure that the PES goes rigorously to zero in monomer+trimer and dimer+dimer dissociations. The 4-b energies of isomers of the hexamer calculated with the new PES are shown to be in better agreement with benchmark CCSD(T) results than those from the MB-pol potential. Tests on larger clusters further validate the high-fidelity of the PES. The PES is shown to be fast to evaluate, taking 2.4 s for 10⁵ evaluations on a single core of 2.4 GHz Intel Xeon processor, and significantly faster using a parallel version of the PES.

Parsimonious Potential Energy Surface Expansions Using Dictionary Learning with Multipass Greedy Selection

Potential energy surfaces fit with basis set expansions have been shown to provide accurate representations of electronic energies and have enabled a variety of high-accuracy dynamics, kinetics, and spectroscopy applications. The number of terms in these expansions scales poorly with system size, a drawback that challenges their use for systems with more than ∼10 atoms. A solution is presented here using dictionary learning. Subsets of the full set of conventional basis functions are optimized using a newly developed multipass greedy regression method inspired by forward and backward selection methods from the statistics, signal processing, and machine learning literatures. The optimized representations have accuracies comparable to the full set but are 1 or more orders of magnitude smaller, and notably, the number of terms in the optimized multipass greedy expansions scales approximately linearly with the number of atoms.

Breaking the Coupled Cluster Barrier for Machine-Learned Potentials of Large Molecules: The Case of 15-Atom Acetylacetone

Machine-learned potential energy surfaces (PESs) for molecules with more than 10 atoms are typically forced to use lower-level electronic structure methods such as density functional theory (DFT) and second-order Møller-Plesset perturbation theory (MP2). While these are efficient and realistic, they fall short of the accuracy of the "gold standard" coupled-cluster method, especially with respect to reaction and isomerization barriers. We report a major step forward in applying a Δ-machine learning method to the challenging case of acetylacetone, whose MP2 barrier height for H-atom transfer is low by roughly 1.1 kcal/mol relative to the benchmark CCSD(T) barrier of 3.2 kcal/mol. From a database of 2151 local CCSD(T) energies and training with as few as 430 energies, we obtain a new PES with a barrier of 3.5 kcal/mol in agreement with the LCCSD(T) barrier of 3.5 kcal/mol and close to the benchmark value. Tunneling splittings due to H-atom transfer are calculated using this new PES, providing improved estimates over previous ones obtained using an MP2-based PES.

Machine learning meets chemical physics

Over recent years, the use of statistical learning techniques applied to chemical problems has gained substantial momentum. This is particularly apparent in the realm of physical chemistry, where the balance between empiricism and physics-based theory has traditionally been rather in favor of the latter. In this guest Editorial for the special topic issue on "Machine Learning Meets Chemical Physics," a brief rationale is provided, followed by an overview of the topics covered. We conclude by making some general remarks.

Full-dimensional potential energy surface for acetylacetone and tunneling splittings

New, full-dimensional potential energy surface for acetylacetone allows for description of H-tunneling dynamics and characterization of stationary points.

Nuclear-Electronic Orbital Multistate Density Functional Theory

Hydrogen tunneling is essential for a wide range of chemical and biological processes. The description of hydrogen tunneling with multicomponent quantum chemistry approaches, where the transferring hydrogen nucleus is treated on the same level as the electrons, is challenging due to the importance of both static and dynamical electron-proton correlation. Herein the nuclear-electronic orbital multistate density functional theory (NEO-MSDFT) method is presented as a strategy to include both types of correlation. In this approach, two localized nuclear-electronic wave functions obtained with the NEO-DFT method are combined with a nonorthogonal configurational interaction approach to produce bilobal, delocalized ground and excited vibronic states. By including a correction function, the NEO-MSDFT approach can produce quantitatively accurate hydrogen tunneling splittings for fixed geometries of systems such as malonaldehyde and acetoacetaldehyde. This approach is computationally efficient and can be combined with methods such as vibronic coupling theory to describe tunneling dynamics and to compute vibronic couplings in many types of systems.