The evolution of machine learning potentials for molecules, reactions and materials

Recent years have witnessed the fast development of machine learning potentials (MLPs) and their widespread applications in chemistry, physics, and material science. By fitting discrete ab initio data faithfully to continuous and symmetry-preserving mathematical forms, MLPs have enabled accurate and efficient atomistic simulations in a large scale from first principles. In this review, we provide an overview of the evolution of MLPs in the past two decades and focus on the state-of-the-art MLPs proposed in the last a few years for molecules, reactions, and materials. We discuss some representative applications of MLPs and the trend of developing universal potentials across a variety of systems. Finally, we outline a list of open challenges and opportunities in the development and applications of MLPs.

Molecular vibrations of acetylacetone enol. Infrared polarization spectroscopy and theoretical predictions

The IR polarization spectrum of acetylacetone enol (AAe, (3Z)-4-hydroxy-3-penten-2-one) was recorded in the region 2000 - 450 cm-1 using stretched polyethylene as an anisotropic solvent. The measured orientation factors were consistent with Cs molecular symmetry of AAe and provided an experimental distinction between in-plane and out-of-plane polarized spectral features. The results suggest the assignment of at least one previously unrecognized fundamental transition. Excellent agreement with presently and previously observed transition wavenumbers in the investigated region was obtained by theoretical predictions in the harmonic approximation using B3LYP-D3 Density Functional Theory (DFT); calculated diagrams are provided for all 39 normal modes of AAe. The prediction of in-plane transition moment directions was found to be very sensitive to calculational details. The wavenumbers computed with the anharmonic models VPT2 (2nd-order Vibrational Perturbation Theory) and DVPT2 (Deperturbed VPT2) were consistent with those of the harmonic analysis, but the predictions obtained using GVPT2 (Generalized VPT2) deviated significantly.

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.

Fundamental Invariant Neural Network (FI-NN) Potential Energy Surface for the OH + CH3OH Reaction with Analytical Forces

The hydrogen abstraction reaction of OH + CH3OH plays a great role in combustion and atmospheric and interstellar chemistry and has been extensively studied theoretically and experimentally. Theoretically, the numerical gradients with respect to the Cartesian coordinates of atoms in molecular simulations on our recent potential energy surface (PES) for the title reaction trained using the permutationally invariant polynomial neural network (PIP-NN) approach hinder the extensive calculation because of the unaffordable computation cost. To address this issue, we in this work report a new full-dimensional accurate analytical PES for the title reaction using the fundamental invariant neural network (FI-NN) approach based on 140,192 points of the quality UCCSD(T)-F12a/AVTZ. Besides, the spin-orbit (SO) corrections of OH in the entrance channel were determined at the level of complete active space self-consistent field with the AVTZ basis set. As a compromise between computational cost and efficiency, the Δ-machine learning approach was employed to construct the SO-corrected PES. Based on this new FI-NN PES with analytical forces, thermal rate coefficients and various dynamic properties, including the integral cross sections, the differential cross sections, and the product energy partitioning, were determined by running a total of 5.5 million trajectories. The use of analytical gradients of the FI-NN PES accelerated simulations and about 99% of computation cost was saved, compared to that for the PIP-NN PES with numerical gradients. Such a significant acceleration is achieved mainly by replacing PIPs with FIs.

Quantum Monte Carlo Simulations of the Vibrational Wavefunction of the Aromatic Cyclo[10]carbon Using a Full Dimensional Permutationally Invariant Potential Energy Surface

New experimental measurements [Sun et al., Nature 2023, 623, 972] of the cyclic C10 reveal a cumulenic pentagon-like D5h structure at ∼5 K. However, the long-standing presumption that a large zero-point vibrational energy combined with an extremely flat D5h ↔ D10h ↔ D5h isomerization pathway washes out the pentagonal D5h structure and yields a symmetric D10h decagon remains at odds with the experiment. We resolve this issue with our fitting approach based on a bond-order charge-density matrix expressed in permutationally invariant polynomials. We train the model on τHCTH/cc-pVQZ data morphed to reproduce a relativistic all-electron CCSDT(Q)/CBS D5h-D10h potential energy barrier (benchmarked previously by others). Large scale diffusion Monte Carlo simulations in full dimensionality show that the vibrational ground state of C10 has compositional character of more than 96% D5h, fully reflecting the experimental imaging data. Quantum mechanical variational calculations in 1-D further suggest persistence of the D5h symmetry structure at higher temperatures.

No Headache for PIPs: A PIP Potential for Aspirin Runs Much Faster and with Similar Precision Than Other Machine-Learned Potentials

Assessments of machine-learning (ML) potentials are an important aspect of the rapid development of this field. We recently reported an assessment of the linear-regression permutationally invariant polynomial (PIP) method for ethanol, using the widely used (revised) rMD17 data set. We demonstrated that the PIP approach outperformed numerous other methods, e.g., ANI, PhysNet, sGDML, and p-KRR, with respect to precision and notably with respect to speed [Houston et al., J. Chem. Phys. 2022, 156, 044120]. Here, we extend this assessment to the 21-atom aspirin molecule, using the rMD17 data set, with a focus on the speed of evaluation. Both energies and forces are used for training, and the precision of several PIPs is examined for both. Normal mode frequencies, the methyl torsional potential, and 1d vibrational energies for an OH stretch are presented. We show that the PIP approach achieves the level of precision obtained from other ML methods, e.g., atom-centered neural network methods, linear regression ACE, and kernel methods, as reported by Kovács et al. in J. Chem. Theory Comput. 2021, 17, 7696-7711. More significantly, we show that the PIP PESs run much faster than all other ML methods, whose timings were evaluated in that paper. We also show that the PIP PES extrapolates well enough to describe several internal motions of aspirin, including an OH stretch.

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.

Investigating Intramolecular H Atom Transfer Dynamics in β-Diketones with Ultrafast Infrared Spectroscopies and Theoretical Modeling

The vibrational signatures and ultrafast dynamics of the intramolecular H-bond in a series of β-diketones are investigated with 2D IR spectroscopy and computational modeling. The chosen β-diketones exhibit a range of H atom donor-acceptor distances and asymmetry along the H atom transfer coordinate that tunes the intramolecular H-bond strength. The species with the strongest H-bonds are calculated to have very soft H atom potentials, resulting in highly red-shifted OH stretch fundamental frequencies and dislocation of the H atom upon vibrational excitation. These soft potentials lead to significant coupling to the other normal mode coordinates and give rise to the very broad vibrational signatures observed experimentally. The 2D IR spectra in both the OH and OD stretch regions of the light and deuterated isotopologues reveal broadened and long-lived ground-state bleach signatures of the vibrationally hot molecules. Polarization-sensitive transient absorption measurements in the OH and OD stretch regions reveal notable isotopic differences in orientational dynamics. Orientational relaxation was measured to occur on ∼600 fs and ∼2 ps time scales for the light and deuterated isotopologues, respectively. The orientational dynamics are interpreted in terms of activated H/D atom transfer events driven by collective intramolecular structural rearrangements.

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.

Δ-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.

Graph-Theoretic Molecular Fragmentation for Potential Surfaces Leads Naturally to a Tensor Network Form and Allows Accurate and Efficient Quantum Nuclear Dynamics

Molecular fragmentation methods have revolutionized quantum chemistry. Here, we use a graph-theoretically generated molecular fragmentation method, to obtain accurate and efficient representations for multidimensional potential energy surfaces and the quantum time-evolution operator, which plays a critical role in quantum chemical dynamics. In doing so, we find that the graph-theoretic fragmentation approach naturally reduces the potential portion of the time-evolution operator into a tensor network that contains a stream of coupled lower-dimensional propagation steps to potentially achieve quantum dynamics with reduced complexity. Furthermore, the fragmentation approach used here has previously been shown to allow accurate and efficient computation of post-Hartree-Fock electronic potential energy surfaces, which in many cases has been shown to be at density functional theory cost. Thus, by combining the advantages of molecular fragmentation with the tensor network formalism, the approach yields an on-the-fly quantum dynamics scheme where both the electronic potential calculation and nuclear propagation portion are enormously simplified through a single stroke. The method is demonstrated by computing approximations to the propagator and to potential surfaces for a set of coupled nuclear dimensions within a protonated water wire problem exhibiting the Grotthuss mechanism of proton transport. In all cases, our approach has been shown to reduce the complexity of representing the quantum propagator, and by extension action of the propagator on an initial wavepacket, by several orders, with minimal loss in accuracy.

Quantum Gaussian process model of potential energy surface for a polyatomic molecule

With gates of a quantum computer designed to encode multi-dimensional vectors, projections of quantum computer states onto specific qubit states can produce kernels of reproducing kernel Hilbert spaces. We show that quantum kernels obtained with a fixed ansatz implementable on current quantum computers can be used for accurate regression models of global potential energy surfaces (PESs) for polyatomic molecules. To obtain accurate regression models, we apply Bayesian optimization to maximize marginal likelihood by varying the parameters of the quantum gates. This yields Gaussian process models with quantum kernels. We illustrate the effect of qubit entanglement in the quantum kernels and explore the generalization performance of quantum Gaussian processes by extrapolating global six-dimensional PESs in the energy domain.

Vibrational Dynamics of the Intramolecular H-Bond in Acetylacetone Investigated with Transient and 2D IR Spectroscopy

Acetylacetone (AcAc) has proven to be a fruitful but highly challenging model system for the experimental and computational interrogation of strong intramolecular hydrogen bonds. Key questions remain, however, regarding the identity of the minimum-energy structure of AcAc and the dynamics of intramolecular proton transfer. Here, we investigate the OH/OD stretch and bend regions of the enol tautomer of AcAc and its deuterated isotopologue with transient absorption and 2D IR spectroscopy. The OH bend region reveals a single dominant diagonal transition near 1625 cm^-1 with intense cross peaks to lower-frequency modes, demonstrating highly mixed fingerprint transitions that contain OH bend character. The anharmonic coupling of the OH bend results in a highly elongated OH bend excited-state absorption transition that indicates a large manifold of OH bend overtone/combination bands in the OH stretch region that leads to strong bend-stretch Fermi resonance interactions. The OH and OD stretch regions consist of broad ground-state bleach signals, but there is no clear evidence of ω₂₁ excited-state absorptions due to rapid population relaxation arising from strong intramolecular coupling to bending, fingerprint, and low-frequency H-bond modes. Orientational relaxation dynamics persist for timescales longer than the vibrational lifetimes, with polarization anisotropy components decaying within approximately 2 and 10 periods of the O-O oscillation for the OH and OD stretch, respectively. The significant isotopic dependence of the orientational dynamics is discussed in the context of intramolecular mode coupling, diffusional processes, and contributions from proton/deuteron transfer dynamics.

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.

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.

Permutationally Invariant Polynomial Expansions with Unrestricted Complexity

A general strategy is presented for constructing and validating permutationally invariant polynomial (PIP) expansions for chemical systems of any stoichiometry. Demonstrations are made for three categories of gas-phase dynamics and kinetics: collisional energy-transfer trajectories for predicting pressure-dependent kinetics, three-body collisions for describing transient van der Waals adducts relevant to atmospheric chemistry, and nonthermal reactivity via quasiclassical trajectories. In total, 30 systems are considered with up to 15 atoms and 39 degrees of freedom. Permutational invariance is enforced in PIP expansions with as many as 13 million terms and 13 permutationally distinct atom types by taking advantage of petascale computational resources. The quality of the PIP expansions is demonstrated through the systematic convergence of in-sample and out-of-sample errors with respect to both the number of training data and the order of the expansion, and these errors are shown to predict errors in the dynamics for both reactive and nonreactive applications. The parallelized code distributed as part of this work enables the automation of PIP generation for complex systems with multiple channels and flexible user-defined symmetry constraints and for automatically removing unphysical unconnected terms from the basis set expansions, all of which are required for simulating complex reactive systems.

Machine Learning for Chemical Reactions

Machine learning (ML) techniques applied to chemical reactions have a long history. The present contribution discusses applications ranging from small molecule reaction dynamics to computational platforms for reaction planning. ML-based techniques can be particularly relevant for problems involving both computation and experiments. For one, Bayesian inference is a powerful approach to develop models consistent with knowledge from experiments. Second, ML-based methods can also be used to handle problems that are formally intractable using conventional approaches, such as exhaustive characterization of state-to-state information in reactive collisions. Finally, the explicit simulation of reactive networks as they occur in combustion has become possible using machine-learned neural network potentials. This review provides an overview of the questions that can and have been addressed using machine learning techniques, and an outlook discusses challenges in this diverse and stimulating field. It is concluded that ML applied to chemistry problems as practiced and conceived today has the potential to transform the way with which the field approaches problems involving chemical reactions, in both research and academic teaching.

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.

Transfer Learning to CCSD(T): Accurate Anharmonic Frequencies from Machine Learning Models

The calculation of the anharmonic modes of small- to medium-sized molecules for assigning experimentally measured frequencies to the corresponding type of molecular motions is computationally challenging at sufficiently high levels of quantum chemical theory. Here, a practical and affordable way to calculate coupled-cluster quality anharmonic frequencies using second-order vibrational perturbation theory (VPT2) from machine-learned models is presented. The approach, referenced as "NN + VPT2", uses a high-dimensional neural network (PhysNet) to learn potential energy surfaces (PESs) at different levels of theory from which harmonic and VPT2 frequencies can be efficiently determined. The NN + VPT2 approach is applied to eight small- to medium-sized molecules (H₂CO, trans-HONO, HCOOH, CH₃OH, CH₃CHO, CH₃NO₂, CH₃COOH, and CH₃CONH₂) and frequencies are reported from NN-learned models at the MP2/aug-cc-pVTZ, CCSD(T)/aug-cc-pVTZ, and CCSD(T)-F12/aug-cc-pVTZ-F12 levels of theory. For the largest molecules and at the highest levels of theory, transfer learning (TL) is used to determine the necessary full-dimensional, near-equilibrium PESs. Overall, NN + VPT2 yields anharmonic frequencies to within 20 cm^-1 of experimentally determined frequencies for close to 90% of the modes for the highest quality PES available and to within 10 cm^-1 for more than 60% of the modes. For the MP2 PESs only ∼60% of the NN + VPT2 frequencies were within 20 cm^-1 of the experiment, with outliers up to ∼150 cm^-1, compared to the experiment. It is also demonstrated that the approach allows to provide correct assignments for strongly interacting modes such as the OH bending and the OH torsional modes in formic acid monomer and the CO-stretch and OH-bend mode in acetic acid.

Application of Reaction Path Search Calculations to Potential Energy Surface Fits

There has been significant progress in recent years in the use of machine learning techniques to model high-dimensional reactive potential energy surfaces using large-scale data obtained from ab initio electronic structure calculations. In these methods, the strategy used to gather data becomes a key issue as the molecular size increases. In this work, we examine the applicability of the reaction path search algorithm implemented in the Global Reaction Route Mapping (GRRM) code as a data-gathering approach. The electronic energies and gradients sampled by using the GRRM calculation are directly used in potential energy surface fitting to a permutationally invariant polynomial function. This simple approach was applied to the HNS and HCNO reaction systems, and we found that the fitted potential energy surfaces reasonably reproduce the features of the electronic structure calculations used in the GRRM calculations. This suggests that the GRRM sampling scheme can be used to construct an initial potential energy surface.