Constrained nuclear-electronic orbital second-order Moller-Plesset perturbation theory

A multicomponent second-order Møller-Plesset perturbation theory (MP2) method is derived and implemented within the constrained nuclear-electronic orbital (cNEO) framework from a multicomponent generalization of the Hylleraas functional. The cNEO-MP2 method includes electronic-nuclear and nuclear correlation in the calculation of vibrationally averaged molecular properties, and is the first post Hartree-Fock wavefunction-based cNEO method. Nuclear quantum effects like vibrational averaging, isotopic effects, and zero-point energy can be captured in a single calculation or geometry optimization with cNEO-MP2, eliminating the need to perform costly subsequent calculations to determine higher order force constants as required with many existing methods used to determine vibrational effects upon molecular properties. The cNEO-MP2 method is benchmarked on a test set of diatomic and small polyatomic molecules and ions. Herein, we present internuclear distances, bond angles, potential energy surfaces, and vibrational frequencies calculated with the cNEO-MP2 method to demonstrate that it correctly captures the effects of nuclear vibrational motion upon molecular properties.

💡 Research Summary

This paper introduces a constrained nuclear‑electronic orbital second‑order Møller‑Plesset perturbation theory (cNEO‑MP2), the first post‑Hartree‑Fock wavefunction method within the constrained NEO (cNEO) framework. Traditional electronic‑structure methods treat nuclei classically under the Born‑Oppenheimer (BO) approximation, which neglects nuclear quantum effects such as zero‑point energy, vibrational averaging, and isotope‑dependent property shifts. Existing remedies—Vibrational Self‑Consistent Field (VSCF), path‑integral approaches, and especially Vibrational Perturbation Theory (VPT)—require high‑order force constants and still rely on a BO reference, limiting accuracy for highly anharmonic or large‑amplitude motions.

Multicomponent quantum chemistry treats selected nuclei quantum mechanically alongside electrons. The original NEO method succeeded in capturing many nuclear quantum effects but required at least two nuclei to be treated classically to avoid translational and rotational contamination, re‑introducing a BO‑like separation that is physically less justified. The cNEO approach resolves this by imposing constraints on the expectation values of the nuclear position operators, thereby fixing the molecular frame while allowing all nuclei to be quantum. These constraints are added to the NEO Lagrangian via Lagrange multipliers, yielding modified Roothaan‑Hall equations for both electronic and nuclear orbitals. The nuclear wavefunction is approximated as a Hartree product, which neglects nuclear exchange—an effect shown to be negligible for most chemical systems—while dramatically reducing the number of two‑particle integrals and the size of nuclear Fock matrices.

Building on cNEO‑HF, the authors derive cNEO‑MP2 by minimizing a multicomponent Hylleraas functional. The functional separates into electronic, electron‑nuclear, and nuclear contributions, each expressed in terms of MP2‑type t‑amplitudes. Because electrons and nuclei are distinguishable, only Coulomb (direct) terms appear in the electron‑nuclear and nuclear‑nuclear parts; exchange terms vanish, and a Kronecker delta distinguishes same‑type from different‑type nuclear pairs. The resulting energy expression is a straightforward extension of the previously published NEO‑MP2 formula, but now evaluated with the constrained nuclear orbitals from cNEO‑HF. Computational scaling remains O(N⁵), identical to conventional MP2, yet the prefactor is reduced because nuclear exchange integrals are omitted and the number of nuclear basis functions is far smaller than the electronic basis.

Implementation details include adding the position‑constraint terms to the SCF loop, solving coupled electronic‑nuclear Fock equations, and then forming the MP2 amplitudes using the constrained orbital energies. The authors benchmark cNEO‑MP2 on a set of diatomics (H₂, HF, LiH, etc.), small polyatomics (H₂O, NH₃, CH₄), and selected ions. They report equilibrium bond lengths, bond angles, potential energy curves, and both harmonic and anharmonic vibrational frequencies. Across the test set, cNEO‑MP2 reproduces experimental geometries and frequencies within chemical accuracy and matches high‑level multicomponent CCSD(T) results. Notably, isotope effects (e.g., H/D substitution) and zero‑point energy contributions emerge automatically from a single cNEO‑MP2 calculation, eliminating the need for separate VPT or finite‑difference force‑constant evaluations.

The authors argue that cNEO‑MP2 offers three decisive advantages: (1) simultaneous treatment of electron‑electron, electron‑nuclear, and nuclear‑nuclear correlation; (2) a clean removal of translational/rotational contamination without classical nuclei; (3) a cost‑effective route to vibrationally averaged properties that rivals or surpasses VPT, especially for systems with strong anharmonicity. Moreover, because the derivation relies on the Hylleraas functional, the formalism can be extended to higher‑order multicomponent methods such as coupled‑cluster singles‑doubles (cNEO‑CCSD) and perturbative triples (cNEO‑CCSD(T)), providing a systematic hierarchy for increasingly accurate nuclear‑quantum‑aware simulations.

In conclusion, cNEO‑MP2 constitutes a robust, scalable, and conceptually clean framework for incorporating nuclear quantum effects directly into electronic‑structure calculations. Its ability to deliver zero‑point energies, vibrational averages, and isotope‑dependent shifts in a single post‑HF step positions it as a compelling alternative to traditional vibrational perturbation techniques. Future work will likely explore its integration with molecular dynamics, excited‑state methods, and larger biochemical systems where proton delocalization and hydrogen‑bond dynamics play a pivotal role.