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.

Introduction

In principle, quantum mechanics (QM) exactly describes the energetics of chemical systems of any size. However, despite the ever increasing computational power available, determining the interactions and forces for large systems required to follow their molecular dynamics becomes computationally prohibitive. This is due to the unfavourable scaling of electronic structure calculations with the number of electrons and the large number of basis functions required for accurately solving the Schrödinger equation although linear scaling methods provide some remedies for these disadvantages.[1] A quantum mechanical treatment for the nuclear degrees of freedom is even more computationally demanding and typical system sizes are ∼10 heavy atoms for which rigorous calculations are feasible.[2]

On the other hand, empirical molecular mechanics (MM) energy functions (“force fields”) are computationally advantageous to evaluate. Propagating the dynamics based on Newton's equations of motion allows to access long time scales for large molecular systems. In this context, “long time scale” is sub-microsecond and “large systems” means 109 atoms for which performance of up to 8 ns day−1 can be achieved.[3,4] Such molecular dynamics (MD) simulations have been used to investigate processes ranging from protein folding,[5,6] ligand binding[7-11] crowding in cellular environments[12] to characterizing spectroscopic properties of solutes and peptides and reactions in the gas phase and in solution.[13-15] However one of the challenges remains to develop suitable energy functions that retain the precision of the QM methods they are often based on and that are suitable to follow bond breaking and bond formation.

One- or multi-dimensional vibrational spectroscopy is a powerful means to characterize the structural dynamics of complex systems.[16,17] Experiments greatly benefit and often require accompanying MD simulations for their molecular-level interpretation. Such atomistic simulations rely on the bonded and non-bonded interactions to be described in a meaningful way for making direct contact with experiments. Traditionally, empirical energy functions use harmonic springs for chemical bonds and valence angles, periodic functions for dihedrals, an atom-centered point charge-based model for charges and a Lennard-Jones representation for van der Waals interactions together with additional, more purpose-tailored terms.[18] For applications in spectroscopy the chemical bonds (stretching vibration) need to be described somewhat more realistically to account for mechanical anharmonicity for which Morse-oscillators are often sufficient because experiments at ambient temperatures are not sensitive to highly vibrationally excited states. Nevertheless, there is scope to use more accurate representations, for example based on machine learning-type approaches, in particular if reference data from high-level electronic structure calculations are available. For the non-bonded interactions – which include electrostatic and van der Waals terms – more physics-based models that go beyond the standard representations have been developed.

The first-order treatment of the electrostatic interaction is based on atom-centered point charges for Coulomb interactions. Such pair interactions can be rapidly computed but lack the accuracy for describing anisotropic contributions to the charge density.[19] Including higher-order atomic multipoles improves the accuracy but at the expense of increased computational cost and implementation complexity.[20-23] Accounting for polarizability is another contribution that has been recently included in empirical force fields and shows much promise for further improvements of the computational models.[24] From an empirical force field perspective the van der Waals interactions are often represented as Lennard-Jones terms with ad hoc (Lorentz–Berthelot) combination rules. Alternative and potentially improved representations are the buffered 14–7 parametrization[25] and modified combination rules.[26,27]

Current methods to investigate reactive systems in the gas- and condensed-phase include mixed quantum mechanical/molecular mechanics (QM/MM),[28-30] reactive force fields ReaxFF[31] and reactive MD (RMD).[32] The general ansatz for reactive force fields is to describe reactant and product states with a separate energy function and to connect the two either by mixing functions or by diagonalizing an n × n matrix, where n is the number of states considered. This reactive part is embedded in an environment that is described by a more empirical energy function, akin to mixed QM/MM calculations and simulations. In recent years, a variety of sophisticated Machine Learning (ML) potentials were developed to accurately represent ab initio results from high level reference calculations (see ref. 33, 34 and references therein). Such ML-based energy functions were also extended to reactive systems in the gas phase,[35] for simulations in solution[36] and on surfaces,[37] and recently they were also presented for composite systems in the gas phase to study complex combustion processes.[38]

Computational investigations of surface reactions are of particular interest from a technological perspective. A substantial amount of heterogeneous catalytic reactions are performed in chemical industry and make up an important economic factor.[39] The focus of theoretical and experimental investigation are not only on formation and breaking of chemical bonds but also the diffusion processes as well as the prediction of scattering experiments.[40-42] The systems considered range from monocrystalline surfaces of metals, metal oxides, minerals or ionic compounds to graphene sheets, amorphous porous carbon[43] and water surfaces.[44,45] They can be modified by adsorbing ultra-thin layers on metal support, using metal alloy or nanostructuring by metal cluster or single atoms or vacancies.[46-48] Of particular relevance are surface sites with coordination-unsaturated atoms on surface edges, kinks or vacancies which have been found to be sites of increased reactivity.[49,50] However, modeling such surfaces with impurities from high-level electronic quantum methods is still a challenging task, in particular if dynamics information is sought.

A combination of methods with the accuracy of a quantum mechanical treatment with a computational performance comparable to an empirical force field would open up possibilities to investigate reactive processes and spectroscopic properties at a quantitative level. Advances in this direction were acknowledged by the 2013 Nobel Prize in Chemistry.[51,52] However, application of accurate, physics-based force fields for condensed-phase simulations is still not routine. The present work highlights specific applications from our own work and from colleagues in the field, demonstrates the opportunities of such approaches, and discusses future prospects and reasons which slow down more rapid adaptation of such methods in a broader sense. This overview of the field of quantitative MD simulations focuses on applications to questions arising in physical chemistry and biophysics. In the current context, the quantitative aspect is judged as from comparing with available experimental data. Therefore, the close relationship between experimental and computational characterization of the systems is essential.

Computational methods for quantitative MD

Following the time evolution of a chemical system by means of atomistic simulations is in principle possible through ab initio MD (AIMD) simulations. However, computational feasibility and the limited accuracy of density functional theory (DFT) methods for certain applications, such as reactions, often preclude using such an approach for quantitative studies. Production runs for full AIMD simulations are typically limited to tens or hundreds of trajectories with simulation times up to hundreds of picoseconds at the semiempirical or DFT level.[53-55] On the other hand, empirical force fields provide a computationally efficient way to determine the total energy and forces of systems in the condensed phase to carry out MD simulations. However, such (parametrized) representations are often not sufficiently accurate for quantitative comparisons or even predictions.

Potential energy surfaces for gas- and condensed-phase processes

With the advent of computationally efficient and accurate electronic structure calculations, empirical potential energy surfaces (PESs) for small molecules can now be replaced with more accurate representations. It is now possible to determine energies for ∼104 geometries at the coupled cluster singles/doubles and perturbative triples (CCSD(T)) or at the multi reference configuration interaction (MRCI) levels of theory. This shifted the main problem to representing this information such that the PES can be evaluated efficiently and with comparable accuracy as the underlying quantum chemical calculations. With this, even “spectroscopically accurate” calculations are now possible for small molecules.[56]

For small molecules, in particular atom–diatom or diatom–diatom van der Waals complexes, expansions in terms of products of Legendre polynomials P(cos θ), single Y(θ, ϕ) or coupled Yl (θ1, θ2, ϕ) spherical harmonics together with radial functions V (R) is convenient.[57,58] Such an approach was used together with explicit fits to experimental data or to reference electronic structure calculations. Alternatively, for the long-range, distance-dependent part explicit electrostatics using experimentally determined or computed atomic and molecular multipoles and polarizabilities and hyperpolarizabilities[59,60] provides very accurate PESs.

Alternatively, permutationally invariant polynomials (PIPs) can be used to describe the total PES.[61] This has the added benefit that it lends itself to more straightforward generalization to globally reactive PESs as has been done to study reactive collisions for N2 + N2 → N2 + 2N and N2 + N2 → 4N.[62] PIPs use a basis of Morse-type functions and fit products of such basis functions to the reference electronic structure data. Recently, PIPs have been used for systems as large as N-methyl acetamide[63] or tropolone.[64]

PESs can also be represented by kernel-based approaches such as a reproducing kernel Hilbert space (RKHS).[65-67] By construction, the RKHS reproduces the reference energies from electronic structure calculations exactly on the grid points. In addition, the physical long-range decay for large separations → ∞ can be explicitly included in the functional dependence of the kernel. Such RKHS-based representations have been used for entire PESs[68] or in a QM/MM-type approach to treat part of an extended system with higher accuracy.[69-71]

Finally, recent efforts for high-accuracy representations of energy functions for individual molecules have revolved around machine learning techniques including neural networks (NNs),[72,73] PIPs combined with NNs,[74] Gaussian processes,[75] or kernel-based methods.[76,77] Recent reviews provide a concise status of this field.[34,78-83]

For condensed phase simulations accurate electrostatic and van der Waals interactions are essential. The most widely adopted approximation for electrostatics is to represent the molecular charge distribution as a superposition of point charges located on the atoms.[18] However, such an approach neglects the local anisotropy of the electrostatic potential (ESP) which can be particularly relevant for halogens or charged groups. Specifically for halogen atoms the σ-hole requires particular attention.[84,85] For such problems, and for more generally representing a molecule's ESP in a realistic fashion, multipolar[21,86] or distributed charge models (DCM)[87,88] were developed over the past 20 years. A non-exhaustive list of notable efforts in this direction are the “Sum of Interactions Between Fragments Ab Initio” (SIBFA),[89] “atomic multipole optimized energetics for biomolecular applications” (AMOEBA),[90] and the “atomic multipole” (MTP)[91] approaches. In addition, energy functions accounting for polarizability have been developed which also allow to obtain more realistic, physics-based models, in particular for condensed-phase simulations.[90,92]

The charge distribution ρ(x) also depends on the geometry x. Developing charge models capturing this variation can be challenging. For one, fluctuating point charges based on the charge equilibration scheme have been developed for this purpose.[93,94] They have been mainly applied to simulate the structure and diffusivity of water, ions or amides in water but not for infrared (IR) spectroscopy. Recently, a charge equilibration scheme using a high dimensional NN to learn chemical hardness and electronegativities was presented.[95] Furthermore, generalizations to simulating larger molecules are difficult and no such extensions are available for multipolar charge models although their conformational dependence has been investigated.[96] For CO in myoglobin[97] and CN− in water[98] conformationally dependent MTPs have demonstrated to perform very well. For larger molecules than diatomics fluctuating charge models have been reported by using ML methods to predict point charges for an ensemble of structures to match their molecular dipole moment, e.g. within the PhysNet neural network approach.[73]

Dynamics and reactions on surfaces

The two established ways for computations involving surfaces are by periodic and cluster embedding models. The first model uses the periodicity of a unit cell to describe the electronic wave function of an infinite system. Periodic models allow meaningful computations for metal surfaces. Care has to be taken to minimize finite size effects such as self interaction of surface modifications or reactive sites. Increasing the size of the unit cell lowers the self interactions but leads to higher computational cost.[99] As realistic systems contain at least ∼102 atoms, simulations for reactions typically use generalized gradient approximation (GGA) functionals.[100] For non-reactive systems also meta-GGA and hybrid density functionals were already applied.[101,102] However, the accuracy can be increased by the specific reaction parameter (SRP) approach to density functional theory that optimizes the mix of exchange and correlation energies.[103] Surface systems are also modeled by a limited cluster surrounded by a sufficiently sized grid of point charges corresponding to the charge of cluster atoms or ions. Cluster embedding is usually restricted to insulators and semi conductors but is well suited for modelling surface impurities, supported surface clusters or reactants of low concentration.[104-106] Thus, with the cluster embedding method even coupled cluster calculations can be carried out.[105]

As for MD simulations in general, the quality of the potential energy surface but also the validity of the Born–Oppenheimer approximation for the surface process significantly determines the accuracy of the simulation.[107] The accuracy of potential energy surface is important for adiabatic processes such as elastic intramolecular energy redistribution within the adsorbate when interacting with the surface and inelastic energy transfer between adsorbate and surface. Experimentally, adsorbate–surface interaction can be probed by state-to-state scattering experiments[108,109] or time dependent desorption.[110-112]

Non-adiabatic effects on the dynamics must be considered for simulations especially for processes on metal surfaces. Here, one “weakly” non-adiabatic effect is called “electronic friction” that is the coupling between moving adsorbate atoms with a manifold of electronic states of the surface via electron–hole pair transitions. It is labelled “weakly” as the effect can be treated perturbatively as a damping force on the adsorbate movement.[107,113-115] Additionally a random force applies on the adsorbate for surfaces of finite temperatures. Effects on the dynamics by coupling between two or more PES of multiple electronic states are accounted for in the “strong” non-adiabatic regime.

Vibrational spectroscopy and dynamics in solution

Vibrational spectroscopy is particularly suitable for quantitative comparisons between experiments and simulations. This is primarily due to the precision with which such (laser-based) experiments can be carried out. In the context of experiments, atomistic simulations also yield positions and velocities of all atoms at all times. This provides necessary information to relate experimental observables, such as the IR spectrum, with specific structural features of a system. Even small peptides at ambient conditions can sample multiple conformations each of which exhibits potentially conformer-specific IR spectra. This applies to both, peptides in the gas phase and in solution.[116-118] Thus it is of interest to determine (a) what conformational substate leads to a particular spectroscopic response and (b) whether there exists a unique correspondence (“fingerprint”) between structure and spectroscopy. One possibility is to use extensive electronic structure calculations. This is, however, time consuming and usually only possible for gas-phase systems.[117,119] Alternatively, MD simulations with physics-based and improved force fields can be used to determine the underlying structural features by comparing computed and experimentally measured IR spectra.[120,121]

N-Methyl acetamide (NMA) is a topical example for which substantial experimental and computational work has been carried out. In a recent effort the frequency correlation function for NMA in water was determined from an energy function based on reproducing kernel Hilbert space representation for the [CONH] moiety and atomic multipoles up to quadrupole for the electrostatics.[15] Depending on the technique used to determine the amide-I stretch frequency (scanning along the –CO local mode or along the CONH normal mode) the three time scales τ1 to τ3 for the decay of the frequency–frequency correlation function (FFCF) were [0.02; 0.21; 1.00] ps or [0.02; 0.20; 0.81] ps compared with two time scales from experiments [(0.05 to 0.1); 1.6] ps and [0.01; 1.0] ps[122,123] and [0.06; 0.66] ps from simulations[123] from using a standard force field. The simulation results using a multipolar representation favour a larger value for the decay times which is more consistent with the experimental findings.

A slightly larger and equally well-studied system is trialanine (Ala3).[120,121,124-133] Most experiments were carried out under conditions that prefer the cationic species and agree that the conformational ensemble is dominated by the poly-proline II (ppII) structure, often together with some population of the β-sheet conformation and rare sampling of a right-handed α-helical structure. From 1-d and 2-d IR experiments in the amide-I region a notable study used spectroscopic data to refine the underlying conformational ensemble generated from MD simulations. For this, the experimentally measured IR spectra was reproduced best from Bayesian ensemble refinement[120] which effectively reweights a reference distribution with corresponding conformer-specific IR spectra. Interestingly, the final ensemble generated from such an approach was similar to the Ramachandran maps from explicit MD simulations using a multipolar representation which also correctly described the IR spectroscopy.[121]

IR spectroscopy can also be used to follow protein assembly and disassembly. Insulin binds to the insulin receptor in its monomeric form but is stored in the body as a zinc-bound hexamer each of which consists of three homodimers. Hence, the stability of the dimer to decay into two monomers is a physiologically relevant property. For human WT insulin the experimentally determined stabilization of the dimer with respect to two separated monomers is ΔG = −7.2 kcal mol−1.[134] Thermodynamic stabilities from free energy simulations are between −8.4 and −11.9 kcal mol−1 and simulations along the minimum energy path yield −12.4 kcal mol−1.[135-137] For pharmaceutical applications modified insulins have been synthesized but their dimerization free energies are unknown and challenging to be determined by standard experimental protocols. Hence, alternative means to assess the thermodynamic stability of insulin dimer are required. One possibility is to use infrared spectroscopy because the insulin dimerization interface involves breaking of several hydrogen bonds involving contacts with the –CO unit that give rise to amide-I spectra.[138,139] Atomistic simulations using multipolar force fields confirm that changes in the association state result in modified amide-I spectra.[140] Based on this, attempts can be made to relate changes in the spectroscopic response with the thermodynamic stability of insulin dimer.

The IR spectroscopy and reaction dynamics of diatomic ligands – including CO and NO – bound to and within Myoglobin (Mb) have been thoroughly investigated.[142] Most characteristically, the CO infrared spectrum is split with the two peaks separated by ∼10 cm−1. The spectroscopic signatures were associated with two distinct conformational substates but their structural assignment remained elusive despite dedicated efforts.[143,144] MD simulations with sufficiently detailed electrostatic models for photodissociated CO provided the necessary accuracy[145,146] and together with experimental mutation studies[147] the more red-shifted peak was associated with the Fe–OC orientation whereas the less red-shifted peak corresponds to the Fe–CO state.[148]Fig. 1 shows the structure of unligated CO in Mb together with the IR spectrum of the free CO ligand (dashed line). State-specific spectra for the two conformational substates (Fe–CO in green; Fe–OC in blue) are shown as solid lines. It is found that the total spectrum, which is the experimental observable, agrees favourably with the measured spectrum and that the substate-specific spectra allow assignment to an Fe–CO and Fe–OC motif. However, although the spectra provide an identification of the two states, the results also imply that while sampling the Fe–CO conformation the spectroscopy is still sensitive to the presence of the Fe–OC state and vice versa. This is due to the low isomerization barrier between the two states. They are separated by a barrier of ∼0.7 kcal mol−1 which is close to the experimentally reported barrier of 0.5 kcal mol−1.[149]

Fig. 1

Photodissociated CO and its infrared spectrum in myoglobin (Mb). Left panel: Mb secondary structure in green with heme (ball-and-stick) and CO (van der Waals spheres) and the heme-Fe (green). The black arrows indicate rotation of the photodissociated ligand in the active site. Right panel: The IR spectrum for photodissociated CO from simulations with a multipolar representation of the electrostatics. The spectrum from the entire trajectory (dashed line) is compared with the spectrum from parts of the trajectory sampling the Fe–CO (green) and the Fe–OC (blue) substates. The inset shows experimental spectra for 12CO in human Mb (glycerol/water), 13CO in human Mb (D2O), 13CO in sperm whale Mb (D2O), and 13CO in Hb (D2O).[141] Reproduced from ref. 141 with permission from American Institute of Physics, copyright 1995.

By their very nature, classical MD simulations neglect nuclear quantum effects such as zero point energy or tunneling. Rigorous treatment of quantum effects in the nuclear dynamics is only possible for small systems with ∼10 heavy atoms.[150] Even for such system sizes the calculations need to make approximations because already for 4 atoms using a product basis set with 10 basis functions per degree of freedom leads to a 106 × 106 dense matrix to be diagonalized. Approximate schemes for including quantum effects for spectroscopic applications include path integral approaches or ring polymer MD (RPMD)[151,152] or centroid MD (CMD) simulations.[153-155] Although computationally convenient, care has to be exercised that the computed spectra do not suffer from artifacts. For example, unusually large red shifts have been found from CMD simulations of liquid water which were ultimately traced back to the so called “curvature problem” in CMD.[156] Similarly, RPMD simulations have been applied to the infrared spectroscopy of a range of systems[157-160] but again care has to be taken to avoid and/or remedy artifacts.[161,162] Also, the infrared lineshapes from RPMD simulations have been found to be unusually broad.[163] In the words of a recent review: “[T]he difference between the various methods [for describing nuclear quantum effects] are comparable to those one should expect from performing simulations on an imperfect potential energy surface.”[164] Although quantum mechanical methods can obtain spectroscopically accurate results for small systems in the gas phase,[165] extending this to larger molecules and systems in solution is highly nontrivial. An interesting recent development are on-the-fly quantum dynamics simulations without precomputed PESs which have been applied, e.g., to the nonadiabatic dynamics of pyrazine.[166] Also, non-equilibrium RPMD and CMD simulations may provide effective routes to improved spectroscopic characterizations of systems in the condensed phase.[167,168]

Quantitative simulations are also possible for thermodynamic properties, such as pure liquid densities or hydration free energies. Such validated parametrizations are required for protein–ligand binding studies, for a molecular-level understanding of chromatography, or for reactions in solvents other than water. Recent progress has been made in the automated parametrization for a library of 430 molecules by focusing on the electrostatic, Lennard-Jones, torsional and 1–4 nonbonded interactions.[169] For the density of the pure liquid the experimental reference values were reproduced with an averaged unsigned error of 1.8% whereas for the enthalpy of vaporization it is 5.9% which both are considerable improvements over earlier parametrizations.[170,171] Nevertheless, to remove systematic differences (∼2 kcal mol−1) between experimental and computed hydration free energies (which were determined with the TIP3P model for water) required introducing an overall solute–solvent scaling factor to increase the solute–solvent by 15% on average.[169] Following a similar optimization protocol for the polarizable Drude model the errors on the pure liquid density and heats of vaporization were 2% and 6%, respectively, whereas for the hydration free energies the error decreased to 0.5 kcal mol−1 together with the SWM4-NDP water model.[172,173]

One recurring theme in assessing the quality of (empirical) energy functions is that a particular parametrization may yield favourable comparison for only one or a few experimentally determined properties (e.g. pure liquid density and vaporization energy) – even for a class or library of compounds – but not yield satisfactory results for other observables (hydration free energy).[172] Such transferability has been explicitly considered for cyanide in water.[98,174,175] Using one single parametrization based on fluctuating multipolar electrostatics[21] it was possible to obtain quantitative results for vibrational energy relaxation, the 1d- and 2d-IR spectroscopy, and the hydration free energy. The computed T1 times were T1 = 22 ± 2 ps and 68 ± 11 ps in H2 O and D2O, respectively, compared with 28 ± 7 ps and 71 ± 3 ps from experiments with T1(H2O)/T1(D2O) = 0.33 vs. T1(H2O)/T1(D2O) = 0.39.[174,176] For the 1d- and 2d-infrared spectroscopy the computed blue shift was 35 cm−1 compared with 44 cm−1 from experiments, the full width at half maximum of 13 cm−1vs. 15 cm−1 from experiment, and the tilt angle depending on waiting time agreed between experiment and simulations.[98,176] The computed hydration free energy of −76 kcal mol−1 compared with −72 to −77 kcal mol−1 from experiment.[175,177] Additional improvements can be obtained from using PES-morphing techniques.[178,179]

For pure water a range of force fields has been developed in the recent past. One of them is the AMOEBA model which uses atomic multipoles.[180,181] The iAMOEBA parametrization is very successful for a wide range of ∼30 properties[181] for which the majority of the computed values are within a few percent of the experimentally reported data. The E3B model follows a different strategy by adding explicit three-body terms, akin to a many body expansion.[182] Similarly, the HBB (Huang, Braams, Bowman) force field also uses a many-body expansion.[183] Finally, the most comprehensive water force field in various phases is probably the MB-Pol model which also builds on multipolar interactions and many-body polarization.[184]

Reaction dynamics

Reactions in the gas phase

One recent field which has witnessed quantitative simulations concerns small-molecule reactions relevant to hypersonics and atmospheric re-entry.[185-188] For this, a simulation strategy involving high-level electronic structure calculations, global and reactive RKHS-based PESs and QCT simulations has been successfully developed and used for a range of atom + diatom reactions.[68,189,190] The focus was primarily on computing thermal reaction rates, final state distributions and vibrational relaxation times over a wide temperature range, up to 20 000 K. A typical PES for one electronic state is based on 104 energies and the reference energies are typically reproduced within a few cm−1 by the RKHS interpolation. Computed thermal and vibrational relaxation rates agree to within a few percent with those measured experimentally which is a good basis for developing more coarse grained models.[191-193]

Malonaldehyde (MA), acetylacetone, and formic acid dimer (FAD) are topical systems for quantitative simulations of gas-phase spectroscopy and reactions. In particular MA and FAD have attracted considerable interest and quantitatively accurate results are available for tunneling splittings. The experimentally determined[194] splitting for MA is 21.6 cm−1 which compares with 23.8 cm−1 from MCTDH calculations and 21.6 cm−1 from Monte Carlo simulations on the same full dimensional PES.[195] For acetylacetone, which is related to MA by substituting hydrogen atoms by methyl-groups, the IR spectroscopy is of particular interest as the barrier for proton transfer is low and leads to signatures in the spectra.[196] Morphing[178,179] a parametrized PES suitable for following proton transfer and comparing the resulting infrared spectrum with that from experiments yields an estimated barrier for proton transfer of 2.35 kcal mol−1; see Fig. 2. Subsequent machine learning found a barrier of 3.25 kcal mol−1 from transfer learning to the PNO-LCCSDT(T)-F12 level of theory.[35]

Fig. 2

Experimental infrared (blue) and computed (green, black/orange, red) power spectrum in the region of the H-transfer mode of acetylacetone. For the computations PESs featuring different barrier heights (1.18, 2.35, 4.70) kcal mol−1 were used in a morphing-type approach[178] to assess the position of the proton transfer band. Best agreement between experimentally measured and computed IR spectra is for a barrier height of 2.35 kcal mol−1, compared with 3.2 kcal mol−1 from CCSD(T) calculations.[196] The signatures around 3000 cm−1 in the experimentally measured spectra are due to the CH stretch vibrations.

For FAD the gas phase infrared spectroscopy contained signatures for double proton transfer (DPT).[197] Comparing computed IR spectra from MD simulations on a morphed PES[178] with those measured experimentally yielded an estimated barrier for DPT of 7.2 kcal mol−1 [197] which compares with 7.3 kcal mol−1 from a subsequent analysis of microwave spectra[198] and 8.2 kcal mol−1 from fitting CCSD(T)-F12a energies calculated with the cc-pVTZ and aug-cc-pVTZ basis sets for H and C/O atoms (CCSD(T)-F12a/haTZ) to permutationally invariant polynomials.[199] A recent full dimensional PES which was transfer learned from the MP2 level of theory to CCSD(T) energies reported[200] a barrier for DPT of 7.92 kcal mol−1 and a dissociation energy for FAD in the gas phase from diffusion Monte Carlo simulations of D0 = −14.23 ± 0.08 kcal mol−1 in excellent agreement with an experimentally determined value of −14.22 ± 0.12 kcal mol−1.[201] These examples illustrate the benefit of accurate, reactive and full-dimensional PESs for quantitative simulations of gas phase processes.

QCT simulations were also applied to the photodissociation of formaldehyde (H2CO) following excitation of the 2143 band. The results from QCT simulations on a PIP-represented PES determined at the MRCI/cc-pVTZ level of theory compared with experiments yielded excellent results.[202] The properties considered included the speed distributions for state-selected CO, and the rotational and vibrational distributions of the CO and H2 products.[202] In a more recent study the kinetic energy release of the three-body breakup from Coulomb explosion experiments was compared with that from QCT simulations. For high kinetic energy particularly good agreement was found whereas for lower energies the agreement was rather more qualitative.[203]

Atomistic simulations can also be successfully used for determining vibrational relaxation rates. This was done[204] for the CO + CO collision based on an accurate 6-dimensional PES (frozen monomers) validated for 68 rovibrational levels with a root mean square error of 0.29 cm−1 between experiment and calculations.[205] The vibrational relaxation rates from QCT simulations for the v = 10 and v = 16 vibrational levels were 2.6 × 1010 cm3 s−1 and 3.1 × 1010 cm3 s−1 compared with 1.9 × 1010 cm3 s−1 and 5.0 × 1010 cm3 s−1 from experiments.[206] Such information is paramount for combustion processes and hypersonics and was also determined in near-quantitative agreement for the N + NO and O + CO systems.[190,207] Also, for the H + O3 → OH + O2 reaction, QCT simulations based on a PIP + NN representation of energies calculated at the MRCI level of theory reported final state vibrational distributions for the OH product peaking at v′ = 8, in near-quantitative agreement with available experiments.[208]

Finally, thermal rates for reactions have also been determined for bimolecular processes. As an example, the kinetics of the OH + HO2 reaction to form O2 + H2O was determined from QCT simulations using a full dimensional PES based on PIP + NN.[209] Over a temperature range between 250 K and 3000 K the computed rate is in good agreement with a range of experiments and reproduces in particular the negative temperature dependence reported at low temperatures.

Reactions in solution

For reactions in solution the additional complication is the presence of an environment that needs to be described by a separate energy function. The most common solvent, water, is particularly challenging to represent and despite immense effort to date no single parametrization outperforming all others is available.[210] In addition, many empirical energy functions which are a common starting point for investigating reactions in solution have been parametrized with one of the simpler water force fields such as TIP3P[211] or SPC/E.[212]

Despite the challenges, meaningful quantitative simulations for certain reaction types in solution are available. One of them concerns the Claisen rearrangement reaction.[213] Using multi state adiabatic reactive dynamics (MS-ARMD) simulations, the transformation from chorismate to prephenate and corresponding smaller molecular systems was investigated,[214] whereas the actual activation free energy did not reproduce the experimentally observed results. This can be achieved from using a dedicated parametrization as was done for EVB-based simulations of the same reaction.[215] With a model parametrized to reproduce the activation free energy in water the barrier height in the protein environment was reproduced to within 0.6 kcal mol−1 compared with experiment.

Another class of reactions that has been investigated in great detail in solution are SN2 reactions. One example is the SN2 reaction of haloalkane dehalogenase (DhlA) in which the nucleophile carboxylate from Asp124 replaces one of the halides of the CH2Cl–CH2Cl substrate, i.e. –COO− + (CH2Cl)2 → –OCO–CH2–CH2Cl + Cl−.[216] Following the empirical valence bond (EVB) approach,[14] first the reference reaction in water was parametrized to reproduce the experimentally observed activation free energy. This model was then used for simulations in the protein DhlA and the computed catalytic effect of 11.6 kcal mol−1 was in good agreement with the experimentally observed value of 11.7 kcal mol−1.[216] For the [Br–CH3–Cl]− reaction in solution parametrizations within the MS-VALBOND[217] and MS-ARMD[218] frameworks have been successfully used.[219] The catalytic effect in going from the gas phase to solution for the forward [Br–CH3 + Cl]− → [Br + CH3–Cl]− reaction is 17.4 kcal mol−1 using MS-ARMD at the MP2 level, compared with 15.2 kcal mol−1 from experiment.

Phosphate hydrolysis reactions have also been investigated at a quantitative level. For a number of substituted methyl phenyl phosphate diesters activation free energies were determined from EVB calculations and compared with experiment.[220] Typical differences between computed and measured activation free energies were below 1 kcal mol−1. From analysis of the MD trajectories it was concluded that phosphate transfer in these systems followed an associative rather than a dissociative pathway. For methyl transfer reactions, recent MS-ARMD simulations also found that the pathway is associative. For one of the systems, Pyr + MeBr experimentally determined barrier heights for acetonitrile and hexane as the solvent were 22.5 kcal mol−1 and 27.6 kcal mol−1 which compare with 23.2 kcal mol−1 and 28.1 kcal mol−1 from MS-ARMD simulations.[221]

Surface reactions

Molecular beam scattering experiments are an important tool to investigate the reaction dynamics at the gas–surface interface and yield information about reaction mechanisms, topological features and the gas–surface interaction potential. Their characterization and accurate prediction of suitable observables such as the final vibrational, rotational or translational energy distribution of the reaction products, reaction rates, or desorption energies from simulations is of great importance not only for understanding but also for designing new heterogeneous catalysts.[222-225]

One prominent process is the scattering of NO from a metal Au(111) surfaces that is focus of research for the past 20 years.[226] The scattering of initially highly vibrationally excited NO molecules (νini = 15) has shown vibrational relaxation of about 36 kcal mol−1 (150 kJ mol−1) or up to 10 vibrational quanta. This differs significantly for NO scattering from insulators such as LiF(001) for which negligible vibrational energy loss was observed which is indicative of non-adiabatic relaxation channels rather than purely mechanical ones.[226,227]

Early simulations by Tully and coworker based on the independent electron surface hopping (IESH) model allow metal-to-molecule charge transfer by including two electronic states in the model Hamiltonian.[228,229] Such an approach only qualitatively captures the large experimentally observed vibrational relaxation. As an example, simulations for the scattering of NO(νini = 16) with an incident energy of Eini = 0.5 eV lead to a most probable population of the νfin = 14 final vibrational state whereas experiments report pronounced relaxation to νfin,exp = 6.[109] Including electron friction in the model to perturbatively describe vibration–electron coupling yields relaxation to νfin = 13.[108,109,113,230] These differences between experiment and simulations were primarily attributed to inaccuracies in the PES which also may explain the observed overestimation of multibounce events in the simulations.[228,231]

Fig. 3a shows one of the most recent improvements towards a quantitative study of NO-vibrational relaxation using accurate energy functions in MD simulations including electron friction (MDEF, open green squares).[232] Vibrational relaxation of scattered NO(νini = 16, Eini = 0.52 eV) on Au(111) leads to a loss of 9 quanta to yield νfin = 7 which is close to the experimental peak at νfin,exp = 6.[226] This corresponds to a vibrational energy loss of ∼23 kcal mol−1 (∼1.0 eV) compared with the experimentally observed 36 kcal mol−1 (150 kJ mol−1). Even without electronic friction (red bars) QCT simulations with the embedded atom neural network (EANN) potential model already predicts a considerably more pronounced relaxation with a peak at νfin = 11. Conversely, when NO(νini = 12, Eini = 0.42 eV) scatters from LiF(001), the product state distribution P (νfin) strongly peaks at νfin = 12 and no vibrational relaxation takes place, see Fig. 3b. This finding primarily reflects the rather isotropic interaction between NO and the LiF surface that governs the dynamics.[232] Similarly, the predicted translational energy loss of scattered NO(νini = 1, Eini = 0.31 eV) of 〈ΔEtrans〉 = 4.8 kcal mol−1 (0.21 eV) agrees well with the experimentally determined value of 4.4 kcal mol−1 (0.19 eV).[233] Hence, a quantitative understanding of such experiments is possible.

Fig. 3

Final vibrational distribution of highly vibrationally excited NO scattered from (a) Au(111) and (b) LiF(001) from adiabatic (red) and non-adiabatic MD simulations including electronic friction (open green squares). Experimental data of are shown in black.[226,227] Reproduced from ref. 232 with permission from American Physical Society, copyright 2021.

These results emphasize the need for high-dimensional and accurate PESs to describe energy transfer between the surface and the impacting molecules together with including non-adiabatic effects especially for simulation on metal surfaces. In the present case the improvement is due to using a machine-learned, reactive PES that correctly captures the decrease in the NO dissociation barrier from ∼161 kcal mol−1 (7.0 eV) in the gas phase to 66 kcal mol−1 (2.86 eV)[232] on Au(111) which is below the vibrational energy at ν = 16 (corresponding to 78 kcal mol−1 or 3.4 eV).[109,234,235] The NO bond elongation up to 1.89 Å leads to softening of the molecular vibration and increased coupling with other surface degrees of freedom.

For CO on Au(111) a NN representation of the Behler–Parrinello type fitted to DFT reference data was used together with QCT simulations to investigate the scattering of vibrationally excited CO(vini = 2) from the surface. The simulation shows predominantly specular scattering with no significant deviation between incidence angle (9°) and scattering angle (peak at 10°), consistent with an experimentally observed scattering angle of 9°. More importantly, the experimentally observed CO(vfin = 1) product could be attributed to rapid nonadiabatic relaxation of chemisorbed CO.[236] The vibrationally relaxed CO can either desorb rapidly or transfer into a long living physisorbed state giving raise to the fast and slow component in the experimental time of flight spectra.[237]

Achieving quantitative simulation of surface reactions is still challenging but highly relevant for developing and improving heterogeneous catalysts.[100,238,239] Although the prediction of the reaction probability for the dissociative chemisorption of molecular H2 on metal surfaces is possible with high accuracy,[240] for practical application reliable prediction of heavier atoms are required. Neglecting surface atom motion, simulations of H2 on Cu(111) yield a reaction probability curve that is shifted by less than 1 kcal mol−1 along the collision energy range up to 19.1 kcal mol−1 (80 kJ mol−1).[241]

One limitation in such simulations concerns the level at which the electronic structure calculations can be carried out due to the system size and time scales required for meaningful simulations. The current state-of-the art for modeling molecule–metal reactions is still DFT with GGA functionals.[100] GGA density functionals have shown good performance for predicting sticking probabilities of dissociative adsorptions if the work function of the metal surface exceeds 7 eV.[242] One successful example is the match between ab initio MD simulations of the dissociative adsorption of CH4 on Ni(111) surfaces at 550 K with molecular beam experiments. The SRP approach on density functionals was used for an accurate interpolation of the dissociation barrier height.[103,243] They achieved chemical accuracy for low incident energies ranging from 24.1 kcal mol−1 up to 28.7 kcal mol−1 (101–120 kJ mol−1).[243] At higher incident energies the differences between simulation and experiment increased. It was argued that this is due to a deficiency in the classical dynamics description to account for quasi-resonant energy transfer. However, quantitative agreement was found with initially excited CH-stretch vibrations (ν = 1) at incident energies up to 38.2 kcal mol−1 (160 kJ mol−1). Subsequently, the development of a 15-dimensional PIP-based PES for CH4 on Ni(111) fitted to nearly 200 000 reference data points has opened up the possibility for statistically quantitative description.[244] Chemical accuracy for the dissociative adsorption of CHD3 on Pt(111) by ab initio MD simulations was also achieved.[245-247]

For metal surface systems with low work functions pure density functionals perform less well for reaction barriers and sticking probabilities of dissociatively adsorbed molecules.[242] One well-studied example system is the adsorption of O2 on Al(111) surfaces. Here, experiments have shown small initial sticking probabilities of 1% at low incident energies (∼0.7 kcal mol−1, 30 meV) and an increase to ∼90% at high incident energies (>13.8 kcal mol−1, 600 meV).[248] Adiabatic MD simulations on PESs from GGA density functionals yield a barrier-less O2 dissociation on the Al(111) surface and predict a 100% sticking probability even for the lowest simulated incident energy of 0.6 kcal mol−1 (25 meV). Only energy surfaces from spin-constrained density functional computations with O2 in the electronic triplet state yield qualitative agreement for the sticking probabilities.[249] Nearly quantitative agreement has been achieved by using a PES determined from embedded correlated wave function theory[250] that treats a molecule–surface embedded subsystem with more accurate electron correlation methods and the extended surface at the DFT level. MD simulations were carried out on a fitted 6-dimensional London–Eyring–Polanyi–Sato PES[251,252] using 700 reference points for O2 on a 10 to 14 atoms embedded Al(111) slab in a periodic supercell on MP2 level. The results yield a simulated probability sticking curve that is about 1.1 kcal mol−1 (0.05 eV) shifted towards higher incident energies.[253]

In another recent effort[254] the vibrationally induced isomerization of CO on NaCl surfaces was investigated from QCT simulations. Experiments have reported the isomerization of OC–NaCl to CO–NaCl for highly vibrationally excited CO. For the largest cluster considered, containing 13 CO molecules adsorbed onto the NaCl surface, isomerization starts with v = 22 quanta in the CO stretch, consistent with experiment. Even after isomerization, the CO adsorbate remains vibrationally excited which is also what experiments report.

Quantitative MD simulations have also been possible for chromatographic systems. The dynamics of single molecules at the solid/liquid interface plays important roles in surface science (adsorption/desorption), material science, and interfacial chemistry. The molecular mechanisms underlying retention in RPLC depends on both, specific and non-specific interactions between the analyte and its environment.[255-260] The non-polar stationary phase is often a functionalized silica surface to which alkyl chains of various lengths are tethered, and the typical mobile phase is a water/methanol or water/acetonitrile mixture. Despite the superficially “simple” chemical composition of such systems, a molecular understanding underlying the separation process is challenging as the system is highly dynamical, heterogeneous and disordered.[261-264] Experimentally determined capacity factors which can be converted into retention times and free energies of binding are available for toluene, benzene, chlorobenzene, and phenol.[265] Using multipolar force fields that reproduce the hydration free energies for the probe molecules,[266] free energy differences ΔΔGPhCl→PhR for mutating –Cl in chlorobenzene to –H (benzene), –OH (phenol), and –CH3 (toluene) were determined inside the chromatographic column. With a 20 : 80 MeOH/H2O solvent mixture these computed free energy changes for toluene, benzene, and phenol are [0.05; 0.57; 1.25] kcal mol−1 compared with [0.06; 0.60; 1.43] kcal mol−1 from experiment and for the 50 : 50 MeOH/H2O mixture they are [0.08; 0.51; 1.01] kcal mol−1 compared with [0.20; 0.61; 1.16] kcal mol−1.[267]

Outlook

In this outlook a number of possible future developments are described. They include conceptual and practical aspects of force fields and further refinements to better capture the underlying physics of the molecules considered.

The majority of atomistic simulations in solution are carried out with rather empirical energy functions for the solvent. This is particularly obvious for protein dynamics using TIP3P[211] or SPC/E[212] models for the solvent which are convenient but not quantitative. One of the main reasons for this is the fact that the corresponding energy functions for the proteins have been parametrized with such water models early on. Only replacing the underlying water model will lead to an imbalance between the water–water and water–protein interactions which is not desirable. Hence, only a rigorous “bottom-up” reparametrization of protein force fields using improved but still computationally efficient water energy functions can address this. Recently, an effort has been made to emulate a wide range of rigid, fixed-charge and polarizable water parametrizations using the minimally distributed charge model[88] within the same implementation.[268] This can serve as a proxy for a unified parametrization platform for small molecules and protein fragments in solution based on improved energy functions for water. Together with ongoing efforts to improve energy functions for small molecules and protein building blocks[90,172] such approaches hold promise to make atomistic simulations quantitative in a broader sense and lead to computing platforms that even allow predictive simulations for large systems.

For chemical reactions in solution the reorganization of the electrons makes an important contribution to the solute–solvent interaction. In the gas phase such effects can be absorbed in the total energy if a global representation of the total energy function is possible. However, in solution, the charge reorganization needs to be included explicitly as a function of geometry because the solvent–solute interaction is described by explicit electrostatics in one or the other way. A recent example concerns the recombination dynamics of CO and oxygen atoms. For the process in gas phase all necessary interactions are accounted for in one global energy function[190] whereas for the process on amorphous solid water the reorganization of the electrons between the CO + O reactant and the CO2 product needs to be included explicitly.[45,269]

Conformationally dependent electrostatics has already been explored to some extent in the past.[270] However, the applications were often restricted to water[93] with the exception of the fluctuating charge implementation in CHARMM.[271] The model itself is based on charge equilibration between bonded atoms based on their electronegativity and chemical hardness rather than first principles electrostatic potentials from electronic structure calculations. For additional accuracy, multipolar energy functions can include conformational dependence of the MTPs which has, however, only been done for diatomics.[97,98]

As an example for the influence of conformationally dependent charges, Fig. 4 reports the CO2 formation probability for different initial conditions of the reactive MD simulations. Recombination simulations were initiated from a CO(center of mass)–O separation of (i) R = 4.0 Å for a previously used[45] fixed point charge model in the collinear OC–O (θ = 180°) conformation and (ii) for separations ranging from R = 4.0 Å to R = 8.0 Å with a fluctuating charge model. With fixed charges[45,269] the oxygen atom interacts more strongly with the water hydrogen atoms due to its charge of qO = −0.2 e. Although some charge transfer between the water surface and the oxygen atom is expected from electronic structure calculations[45] the magnitude of the charge is probably too high. On the other hand, with the fluctuating charge model it is possible to correctly describe the change in the charges of the CO and O moieties as they approach one another. In this model, the magnitude of the point charges depends on the OC–O distance. Charges from reference electronic structure calculations were found to change in a sigmoidal fashion between the CO2 product geometry and the CO + O reactant state.

Fig. 4

Normalized reaction time distributions for CO + O recombination from 500 simulations, for reactions initiated from θ = 180° and R = 4.0, 6.0, 8.0 Å using fixed[45,269] and fluctuating point charge (FPC) models. The inset reports the distribution for R = 4.0 Å and almost all recombination takes place within 1 ps. The reaction probability found with R = 4.0 Å is 32% and 99.2% with fixed and fluctuating charges respectively.

As is shown in the inset of Fig. 4, for the shortest separation the normalized reaction time distributions for the two charge models are close but the total product formation probability differs considerably. It is 32% for fixed charges and 99% for fluctuating charges from running 500 simulations. As in the reactant geometry the CO and O fragments are close to electrically neutral they interact weakly with the surface and diffusion is facile (diffusional barrier of 0.2 kcal mol−1 to 2 kcal mol−1 for atomic oxygen)[272] which increases the probability for reactive encounters. Hence, recombination from short distances R is almost 100% and recombination from further away is still possible on the 100 ps time scale, see Fig. 4. The desorption energy for atomic oxygen with qO = −0.2 e is ∼7 kcal mol−1 which also explains its slow surface diffusion with fixed charges. Experimentally, the desorption energy ranges from 2.4 to 3.5 kcal mol−1 [273] compared with 2.2 kcal mol−1 from simulations with neutral oxygen (qO = 0 e).[272] For molecular oxygen (O2) atomistic simulations found[274] desorption energies of 1.5 to 2.0 kcal mol−1 which compares favourably with experiments that report a value of 1.8 kcal mol−1.[275] In summary, with fluctuating charges correctly describing the asymptotic states of the OC + O → CO2 reaction the diffusional and desorption barriers for the oxygen atom are in quantitative agreement with experiment and allow to realistically model recombination on amorphous solid water.

For energy transfer and spectroscopy the coupling between internal degrees of freedom is important. However, in most empirical energy functions this mechanical coupling is implicit and primarily governed by coordinate transformations between Cartesian coordinates which are used for following the molecular dynamics and the internal coordinates in which the energy function itself is described. There are examples for force fields such as the Merck Molecular Force Field (MMFF),[276] which include cross terms, e.g. between neighboring bonds, bonds and valence angles. However, the number of parameters to be determined increases considerably and the additional number of terms affects the computational performance. ML-based techniques provide an opportunity to include such couplings from rigorous electronic structure data and first examples demonstrate their accuracy for chemical reactions and spectroscopy.[200,277]

Dynamical simulations and the chemical understanding of surface processes are a crucial factor for innovations in heterogeneous catalysis. The difficulties and large effort in performing accurate, high-level quantum electronic computations together with incorporating non-adiabatic effects during molecule–surface interactions has often hindered quantitative agreement between simulations and experimental observations. Increased understanding in surface effects, sophisticated models for accurate potential energy surfaces and ML potentials that allow the simulation of adsorbate and surface atom dynamics have lead to increasing success of reproducing experimental results. However, there are examples like the dissociative chemisorption of HCl on Au(111) for which quantitative molecular dynamics simulation could not successfully achieved yet.[37,278-281]

Quantitative atomistic simulations provide considerable scope for molecular-level understanding of complex chemical and biological systems. For this, computationally efficient, versatile and implementations of sufficiently accurate representations of the total energy of the systems are required. The “quantitative” nature of such simulations is ultimately judged from comparison with experiments which in itself has measurement errors associated with it. An integrated approach combining knowledge from experiment and simulation together with a realistic assessment of the uncertainties involved in both will be of particular interest to arrive at a comprehensive description and understanding of complex systems. In this regard, ML-based techniques will contribute to both, the representation of the energy function and the molecular simulations themselves, and the quantification of the underlying uncertainties.

Conflicts of interest

There are no conflicts to declare.