Systematic improvement of empirical energy functions in the era of machine learning.

The impact of targeted replacement of individual terms in empirical force fields is quantitatively assessed for pure water, dichloromethane (CH 2 Cl 2 ), and solvated K + and Cl - ions. For the electrostatic interactions, point charges (PCs) and machine learning (ML)-based minimally distributed charges (MDCM) fitted to the molecular electrostatic potential are evaluated together with electrostatics based on the Coulomb integral. The impact of explicitly including second-order terms is investigated by adding a fragment molecular orbital (FMO)-derived polarization energy to an existing force field, in this case CHARMM. It is demonstrated that anisotropic electrostatics reduce the RMSE for water (by 1.4 kcal/mol), CH 2 Cl 2 (by 0.8 kcal/mol) and for solvated Cl - clusters (by 0.4 kcal/mol). An additional polarization term can be neglected for CH 2 Cl 2 but further improves the models for pure water (by ∼ 1.0 kcal/mol) and hydrated Cl - (by 0.4 kcal/mol), and is key for solvated K + , reducing the RMSE by 2.3 kcal/mol. A 12-6 Lennard-Jones functional form performs satisfactorily with PC and MDCM electrostatics, but is not appropriate for descriptions that account for the electrostatic penetration energy. The importance of many-body contributions is assessed by comparing a strictly 2-body approach with self-consistent reference data. Two-body interactions suffice for CH 2 Cl 2 whereas water and solvated K + and Cl - ions require explicit many-body corrections. Finally, a many-body-corrected dimer potential energy surface exceeds the accuracy attained using a conventional empirical force field, potentially reaching that of an FMO calculation. The present work systematically quantifies which terms improve the performance of an existing force field and what reference data to use for parametrizing these terms in a tractable fashion for ML fitting of pure and heterogeneous systems.

OH-Formation following vibrationally induced reaction dynamics of H2COO.

The reaction dynamics of H2COO to form HCOOH and dioxirane as first steps for OH-elimination is quantitatively investigated. Using a machine learned potential energy surface (PES) at the CASPT2/aug-cc-pVTZ level of theory vibrational excitation along the CH-normal mode νCH with energies up to 40.0 kcal mol-1 (∼5νCH) leads almost exclusively to HCOOH which further decomposes into OH + HCO. Although the barrier to form dioxirane is only 21.4 kcal mol-1 the reaction probability to form dioxirane is two orders of magnitude lower if the CH-stretch mode is excited. Following the dioxirane-formation pathway is facile, however, if the COO-bend vibration is excited together with energies equivalent to ∼2νCH or ∼3νCOO. For OH-formation in the atmosphere the pathway through HCOOH is probably most relevant because the alternative pathways (through dioxirane or formic acid) involve several intermediates that can de-excite through collisions, relax via internal vibrational relaxation (IVR), or pass through loose and vulnerable transition states (formic acid). This work demonstrates how, by selectively exciting particular vibrational modes, it is possible to dial into desired reaction channels with a high degree of specificity.

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.

Low-temperature kinetics for the N + NO reaction: experiment guides the way.

The reaction N(4S) + NO(X2Π) → O(3P) + N2(X1Σ+g) plays a pivotal role in the conversion of atomic to molecular nitrogen in dense interstellar clouds and in the atmosphere. Here we report a joint experimental and computational investigation of the N + NO reaction with the aim of providing improved constraints on its low temperature reactivity. Thermal rates were measured over the 50 to 296 K range in a continuous supersonic flow reactor coupled with pulsed laser photolysis and laser induced fluorescence for the production and detection of N(4S) atoms, respectively. With decreasing temperature, the experimentally measured reaction rate was found to monotonously increase up to a value of (6.6 ± 1.3) × 10-11 cm3 s-1 at 50 K. To confirm this finding, quasi-classical trajectory simulations were carried out on a previously validated, full-dimensional potential energy surface (PES). However, around 50 K the computed rates decreased which required re-evaluation of the reactive PES in the long-range part due to a small spurious barrier with a height of ∼40 K in the entrance channel. By exploring different correction schemes the measured thermal rates can be adequately reproduced, displaying a clear negative temperature dependence over the entire temperature range. The possible astrochemical implications of an increased reaction rate at low temperature are also discussed.

Conformational and state-specific effects in reactions of 2,3-dibromobutadiene with Coulomb-crystallized calcium ions.

Recent advances in experimental methodology enabled studies of the quantum-state- and conformational dependence of chemical reactions under precisely controlled conditions in the gas phase. Here, we generated samples of selected gauche and s-trans 2,3-dibromobutadiene (DBB) by electrostatic deflection in a molecular beam and studied their reaction with Coulomb crystals of laser-cooled Ca+ ions in an ion trap. The rate coefficients for the total reaction were found to strongly depend on both the conformation of DBB and the electronic state of Ca+. In the (4p)2P1/2 and (3d)2D3/2 excited states of Ca+, the reaction is capture-limited and faster for the gauche conformer due to long-range ion-dipole interactions. In the (4s)2S1/2 ground state of Ca+, the reaction rate for s-trans DBB still conforms with the capture limit, while that for gauche DBB is strongly suppressed. The experimental observations were analysed with the help of adiabatic capture theory, ab initio calculations and reactive molecular dynamics simulations on a machine-learned full-dimensional potential energy surface of the system. The theory yields near-quantitative agreement for s-trans-DBB, but overestimates the reactivity of the gauche-conformer compared to the experiment. The present study points to the important role of molecular geometry even in strongly reactive exothermic systems and illustrates striking differences in the reactivity of individual conformers in gas-phase ion-molecule reactions.

The first HyDRA challenge for computational vibrational spectroscopy.

Vibrational spectroscopy in supersonic jet expansions is a powerful tool to assess molecular aggregates in close to ideal conditions for the benchmarking of quantum chemical approaches. The low temperatures achieved as well as the absence of environment effects allow for a direct comparison between computed and experimental spectra. This provides potential benchmarking data which can be revisited to hone different computational techniques, and it allows for the critical analysis of procedures under the setting of a blind challenge. In the latter case, the final result is unknown to modellers, providing an unbiased testing opportunity for quantum chemical models. In this work, we present the spectroscopic and computational results for the first HyDRA blind challenge. The latter deals with the prediction of water donor stretching vibrations in monohydrates of organic molecules. This edition features a test set of 10 systems. Experimental water donor OH vibrational wavenumbers for the vacuum-isolated monohydrates of formaldehyde, tetrahydrofuran, pyridine, tetrahydrothiophene, trifluoroethanol, methyl lactate, dimethylimidazolidinone, cyclooctanone, trifluoroacetophenone and 1-phenylcyclohexane-cis-1,2-diol are provided. The results of the challenge show promising predictive properties in both purely quantum mechanical approaches as well as regression and other machine learning strategies.

Experimental and Theoretical Studies of Hyperthermal N + O2 Collisions.

The dynamics of hyperthermal N(4S) + O2 collisions were investigated both experimentally and theoretically. Crossed molecular beams experiments were performed at an average center-of-mass (c.m.) collision energy of ⟨Ecoll⟩ = 77.5 kcal mol-1, with velocity- and angle-resolved product detection by a rotatable mass spectrometer detector. Nonreactive (N + O2) and reactive (NO + O) product channels were identified. In the c.m. reference frame, the nonreactively scattered N atoms and reactively scattered NO molecules were both directed into the forward direction with respect to the initial direction of the reagent N atoms. On average, more than 90% of the available energy (⟨Eavl⟩ = 77.5 kcal mol-1) was retained in translation of the nonreactive products (N + O2), whereas a much smaller fraction of the available energy for the reactive pathway (⟨Eavl⟩ = 109.5 kcal mol-1) went into translation of the NO + O products, and the distribution of translational energies for this channel was broad, indicating extensive internal excitation in the nascent NO molecules. The experimentally derived c.m. translational energy and angular distributions of the reactive products suggested at least two dynamical pathways to the formation of NO + O. Quasiclassical trajectory (QCT) calculations were performed with a collision energy of Ecoll = 77 kcal mol-1 using two sets of potential energy surfaces, denoted as PES-I and PES-II, and these theoretical results were compared to each other and to the experimental results. PES-I is a reproducing kernel Hilbert space (RKHS) representation of multireference configurational interaction (MRCI) energies, while PES-II is a many-body permutation invariant polynomial (MB-PIP) fit of complete active space second order perturbation (CASPT2) points. The theoretical investigations were both consistent with the experimental suggestion of two dynamical pathways to produce NO + O, where reactive collisions may proceed on the doublet (12A') and quartet (14A') surfaces. When analyzed with this theoretical insight, the experimental c.m. translational energy and angular distributions were in reasonably good agreement with those predicted by the QCT calculations, although minor differences were observed which are discussed. Theoretical translational energy and angular distributions for the nonreactive N + O2 products matched the experimental translational energy and angular distributions almost quantitatively. Finally, relative yields for the nonreactive and reactive scattering channels were determined from the experiment and from both theoretical methods, and all results are in reasonable agreement.

Transfer-learned potential energy surfaces: Toward microsecond-scale molecular dynamics simulations in the gas phase at CCSD(T) quality.

The rise of machine learning has greatly influenced the field of computational chemistry and atomistic molecular dynamics simulations in particular. One of its most exciting prospects is the development of accurate, full-dimensional potential energy surfaces (PESs) for molecules and clusters, which, however, often require thousands to tens of thousands of ab initio data points restricting the community to medium sized molecules and/or lower levels of theory (e.g., density functional theory). Transfer learning, which improves a global PES from a lower to a higher level of theory, offers a data efficient alternative requiring only a fraction of the high-level data (on the order of 100 are found to be sufficient for malonaldehyde). This work demonstrates that even with Hartree-Fock theory and a double-zeta basis set as the lower level model, transfer learning yields coupled-cluster single double triple [CCSD(T)]-level quality for H-transfer barrier energies, harmonic frequencies, and H-transfer tunneling splittings. Most importantly, finite-temperature molecular dynamics simulations on the sub-µs time scale in the gas phase are possible and the infrared spectra determined from the transfer-learned PESs are in good agreement with the experiment. It is concluded that routine, long-time atomistic simulations on PESs fulfilling CCSD(T)-standards become possible.

Interaction at a distance: Xenon migration in Mb.

The transport of ligands, such as NO or O2, through internal cavities is essential for the function of globular proteins, including hemoglobin, myoglobin (Mb), neuroglobin, truncated hemoglobins, or cytoglobin. For Mb, several internal cavities (Xe1 through Xe4) were observed experimentally and they were linked to ligand storage. The present work determines barriers for xenon diffusion and relative stabilization energies for the ligand in the initial and final pocket, linking a transition depending on the occupancy state of the remaining pockets from both biased and unbiased molecular dynamics simulations. It is found that the energetics of a particular ligand migration pathway may depend on the direction in which the transition is followed and the occupancy state of the other cavities. Furthermore, the barrier height for a particular transition can depend in a non-additive fashion on the occupancy of either cavity A or B or simultaneous population of both cavities, A and B. Multiple repeats for the Xe1 → Xe2 transition reveal that the activation barrier is a distribution of barrier heights rather than one single value, which is confirmed by a distribution of transition times for the same transition from unbiased simulations. Dynamic cross correlation maps demonstrate that correlated motions occur between adjacent residues or through space, residue Phe138 is found to be a gate for the Xe1 → Xe2 transition, and the volumes of the internal cavities vary along the diffusion pathway, indicating that there is dynamic communication between the ligand and the protein. These findings suggest that Mb is an allosteric protein.

PhysNet meets CHARMM: A framework for routine machine learning/molecular mechanics simulations.

Full-dimensional potential energy surfaces (PESs) based on machine learning (ML) techniques provide a means for accurate and efficient molecular simulations in the gas and condensed phase for various experimental observables ranging from spectroscopy to reaction dynamics. Here, the MLpot extension with PhysNet as the ML-based model for a PES is introduced into the newly developed pyCHARMM application programming interface. To illustrate the conception, validation, refining, and use of a typical workflow, para-chloro-phenol is considered as an example. The main focus is on how to approach a concrete problem from a practical perspective and applications to spectroscopic observables and the free energy for the -OH torsion in solution are discussed in detail. For the computed IR spectra in the fingerprint region, the computations for para-chloro-phenol in water are in good qualitative agreement with experiment carried out in CCl4. Moreover, relative intensities are largely consistent with experimental findings. The barrier for rotation of the -OH group increases from ∼3.5 kcal/mol in the gas phase to ∼4.1 kcal/mol from simulations in water due to favorable H-bonding interactions of the -OH group with surrounding water molecules.

Water dynamics around T0 vs R4 of hemoglobin from local hydrophobicity analysis.

The local hydration around tetrameric hemoglobin (Hb) in its T0 and R4 conformational substates is analyzed based on molecular dynamics simulations. Analysis of the local hydrophobicity (LH) for all residues at the α1ß2 and α2ß1 interfaces, responsible for the quaternary T → R transition, which is encoded in the Monod-Wyman-Changeux model, as well as comparison with earlier computations of the solvent accessible surface area, makes clear that the two quantities measure different aspects of hydration. Local hydrophobicity quantifies the presence and structure of water molecules at the interface, whereas "buried surface" reports on the available space for solvent. For simulations with Hb frozen in its T0 and R4 states, the correlation coefficient between LH and buried surface is 0.36 and 0.44, respectively, but it increases considerably if the 95% confidence interval is used. The LH with Hb frozen and flexible changes little for most residues at the interfaces but is significantly altered for a few select ones: Thr41α, Tyr42α, Tyr140α, Trp37ß, Glu101ß (for T0) and Thr38α, Tyr42α, Tyr140α (for R4). The number of water molecules at the interface is found to increase by ∼25% for T0 → R4, which is consistent with earlier measurements. Since hydration is found to be essential to protein function, it is clear that hydration also plays an essential role in allostery.

Molecular-level understanding of the rovibrational spectra of N2O in gaseous, supercritical, and liquid SF6 and Xe.

The transition between the gas-, supercritical-, and liquid-phase behavior is a fascinating topic, which still lacks molecular-level understanding. Recent ultrafast two-dimensional infrared spectroscopy experiments suggested that the vibrational spectroscopy of N2O embedded in xenon and SF6 as solvents provides an avenue to characterize the transitions between different phases as the concentration (or density) of the solvent increases. The present work demonstrates that classical molecular dynamics (MD) simulations together with accurate interaction potentials allows us to (semi-)quantitatively describe the transition in rotational vibrational infrared spectra from the P-/R-branch line shape for the stretch vibrations of N2O at low solvent densities to the Q-branch-like line shapes at high densities. The results are interpreted within the classical theory of rigid-body rotation in more/less constraining environments at high/low solvent densities or based on phenomenological models for the orientational relaxation of rotational motion. It is concluded that classical MD simulations provide a powerful approach to characterize and interpret the ultrafast motion of solutes in low to high density solvents at a molecular level.

Structure, Organization, and Heterogeneity of Water-Containing Deep Eutectic Solvents.

The spectroscopy and structural dynamics of a deep eutectic mixture (KSCN/acetamide) with varying water content is investigated from 2D IR (with the C-N stretch vibration of the SCN- anions as the reporter) and THz spectroscopy. Molecular dynamics simulations correctly describe the nontrivial dependence of both spectroscopic signatures depending on water content. For the 2D IR spectra, the MD simulations relate the steep increase in the cross-relaxation rate at high water content to the parallel alignment of packed SCN- anions. Conversely, the nonlinear increase of the THz absorption with increasing water content is mainly attributed to the formation of larger water clusters. The results demonstrate that a combination of structure-sensitive spectroscopies and molecular dynamics simulations provides molecular-level insights into the emergence of heterogeneity of such mixtures by modulating their composition.

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

Correction for 'Transfer learned potential energy surfaces: accurate anharmonic vibrational dynamics and dissociation energies for the formic acid monomer and dimer' by Silvan Käser et al., Phys. Chem. Chem. Phys., 2022, 24, 5269-5281, https://doi.org/10.1039/D1CP04393E.

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

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

Double proton transfer in hydrated formic acid dimer: Interplay of spatial symmetry and solvent-generated force on reactivity.

The double proton transfer (DPT) reaction in the hydrated formic acid dimer (FAD) is investigated at molecular-level detail. For this, a global and reactive machine learned (ML) potential energy surface (PES) is developed to run extensive (more than 100 ns) mixed ML/MM molecular dynamics (MD) simulations in explicit molecular mechanics (MM) solvent at MP2-quality for the solute. Simulations with fixed - as in a conventional empirical force field - and conformationally fluctuating - as available from the ML-based PES - charge models for FAD show a significant impact on the competition between DPT and dissociation of FAD into two formic acid monomers. With increasing temperature the barrier height for DPT in solution changes by about 10% (∼1 kcal mol-1) between 300 K and 600 K. The rate for DPT is largest, ∼1 ns-1, at 350 K and decreases for higher temperatures due to destabilisation and increased probability for dissociation of FAD. The water solvent is found to promote the first proton transfer by exerting a favourable solvent-induced Coulomb force along the O-H⋯O hydrogen bond whereas the second proton transfer is significantly controlled by the O-O separation and other conformational degrees of freedom. Double proton transfer in hydrated FAD is found to involve a subtle interplay and balance between structural and electrostatic factors.

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.

Hydration dynamics and IR spectroscopy of 4-fluorophenol.

Halogenated groups are relevant in pharmaceutical applications and potentially useful spectroscopic probes for infrared spectroscopy. In this work, the structural dynamics and infrared spectroscopy of para-fluorophenol (F-PhOH) and phenol (PhOH) is investigated in the gas phase and in water using a combination of experiment and molecular dynamics (MD) simulations. The gas phase and solvent dynamics around F-PhOH and PhOH is characterized from atomistic simulations using empirical energy functions with point charges or multipoles for the electrostatics, Machine Learning (ML) based parametrizations and with full ab initio (QM) and mixed Quantum Mechanical/Molecular Mechanics (QM/MM) simulations with a particular focus on the CF- and OH-stretch region. The CF-stretch band is heavily mixed with other modes whereas the OH-stretch in solution displays a characteristic high-frequency peak around 3600 cm-1 most likely associated with the -OH group of PhOH and F-PhOH together with a characteristic progression below 3000 cm-1 due to coupling with water modes which is also reproduced by several of the simulations. Solvent and radial distribution functions indicate that the CF-site is largely hydrophobic except for simulations using point charges which renders them unsuited for correctly describing hydration and dynamics around fluorinated sites. The hydrophobic character of the CF-group is particularly relevant for applications in pharmaceutical chemistry with a focus on local hydration and interaction with the surrounding protein.

Quantum and quasi-classical dynamics of the C(3P) + O2(3Σ-g) → CO(1Σ+) + O(1D) reaction on its electronic ground state.

The dynamics of the C(3P) + O2(3Σ-g) → CO(1Σ+) + O(1D) reaction on its electronic ground state is investigated by using time-dependent wave packet propagation (TDWP) and quasi-classical trajectory (QCT) simulations. For the moderate collision energies considered (Ec = 0.001 to 0.4 eV, corresponding to a range from 10 K to 4600 K) the total reaction probabilities from the two different treatments of the nuclear dynamics agree very favourably. The undulations present in P(E) from the quantum mechanical treatment can be related to stabilization of the intermediate CO2 complex with lifetimes on the 0.05 ps time scale. This is also confirmed from direct analysis of the TDWP simulations and QCT trajectories. Product diatom vibrational and rotational level resolved state-to-state reaction probabilities from TDWP and QCT simulations agree well except for the highest product vibrational states (v' ≥ 15) and for the lowest product rotational states (j' ≤ 10). Opening of the product vibrational level CO(v' = 17) requires ∼0.2 eV from QCT and TDWP simulations with O2(j = 0) and decreases to 0.04 eV if all initial rotational states are included in the QCT analysis, compared with Ec > 0.04 eV obtained from experiments. It is thus concluded that QCT simulations are suitable for investigating and realistically describe the C(3P) + O2(3Σ-g) → CO(1Σ+) + O(1D) reaction down to low collision energies when compared with results from a quantum mechanical treatment using TDWPs.

Combining Machine Learning and Spectroscopy to Model Reactive Atom + Diatom Collisions.

The prediction of product translational, vibrational, and rotational energy distributions for arbitrary initial conditions for reactive atom + diatom collisions is of considerable practical interest in atmospheric re-entry. Because of the large number of accessible states, determination of the necessary information from explicit (quasi-classical or quantum) dynamics studies is impractical. Here, a machine-learned (ML) model based on translational energy and product vibrational states assigned from a spectroscopic, ro-vibrational coupled energy expression based on the Dunham expansion is developed and tested quantitatively. All models considered in this work reproduce final state distributions determined from quasi-classical trajectory (QCT) simulations with R2 ∼ 0.98. As a further validation, thermal rates determined from the machine-learned models agree with those from explicit QCT simulations and demonstrate that the atomistic details are retained by the machine learning which makes them suitable for applications in more coarse-grained simulations. More generally, it is found that ML is suitable for designing robust and accurate models from mixed computational/experimental data which may also be of interest in other areas of the physical sciences.