Biomolecular dynamics with machine-learned quantum-mechanical force fields trained on diverse chemical fragments.

The GEMS method enables molecular dynamics simulations of large heterogeneous systems at ab initio quality.

Automatic identification of chemical moieties.

In recent years, the prediction of quantum mechanical observables with machine learning methods has become increasingly popular. Message-passing neural networks (MPNNs) solve this task by constructing atomic representations, from which the properties of interest are predicted. Here, we introduce a method to automatically identify chemical moieties (molecular building blocks) from such representations, enabling a variety of applications beyond property prediction, which otherwise rely on expert knowledge. The required representation can either be provided by a pretrained MPNN, or be learned from scratch using only structural information. Beyond the data-driven design of molecular fingerprints, the versatility of our approach is demonstrated by enabling the selection of representative entries in chemical databases, the automatic construction of coarse-grained force fields, as well as the identification of reaction coordinates.

Accurate global machine learning force fields for molecules with hundreds of atoms.

Global machine learning force fields, with the capacity to capture collective interactions in molecular systems, now scale up to a few dozen atoms due to considerable growth of model complexity with system size. For larger molecules, locality assumptions are introduced, with the consequence that nonlocal interactions are not described. Here, we develop an exact iterative approach to train global symmetric gradient domain machine learning (sGDML) force fields (FFs) for several hundred atoms, without resorting to any potentially uncontrolled approximations. All atomic degrees of freedom remain correlated in the global sGDML FF, allowing the accurate description of complex molecules and materials that present phenomena with far-reaching characteristic correlation lengths. We assess the accuracy and efficiency of sGDML on a newly developed MD22 benchmark dataset containing molecules from 42 to 370 atoms. The robustness of our approach is demonstrated in nanosecond path-integral molecular dynamics simulations for supramolecular complexes in the MD22 dataset.

SpookyNet: Learning force fields with electronic degrees of freedom and nonlocal effects.

Machine-learned force fields combine the accuracy of ab initio methods with the efficiency of conventional force fields. However, current machine-learned force fields typically ignore electronic degrees of freedom, such as the total charge or spin state, and assume chemical locality, which is problematic when molecules have inconsistent electronic states, or when nonlocal effects play a significant role. This work introduces SpookyNet, a deep neural network for constructing machine-learned force fields with explicit treatment of electronic degrees of freedom and nonlocality, modeled via self-attention in a transformer architecture. Chemically meaningful inductive biases and analytical corrections built into the network architecture allow it to properly model physical limits. SpookyNet improves upon the current state-of-the-art (or achieves similar performance) on popular quantum chemistry data sets. Notably, it is able to generalize across chemical and conformational space and can leverage the learned chemical insights, e.g. by predicting unknown spin states, thus helping to close a further important remaining gap for today's machine learning models in quantum chemistry.

Impact of the Characteristics of Quantum Chemical Databases on Machine Learning Prediction of Tautomerization Energies.

An essential aspect for adequate predictions of chemical properties by machine learning models is the database used for training them. However, studies that analyze how the content and structure of the databases used for training impact the prediction quality are scarce. In this work, we analyze and quantify the relationships learned by a machine learning model (Neural Network) trained on five different reference databases (QM9, PC9, ANI-1E, ANI-1, and ANI-1x) to predict tautomerization energies from molecules in Tautobase. For this, characteristics such as the number of heavy atoms in a molecule, number of atoms of a given element, bond composition, or initial geometry on the quality of the predictions are considered. The results indicate that training on a chemically diverse database is crucial for obtaining good results and also that conformational sampling can partly compensate for limited coverage of chemical diversity. The overall best-performing reference database (ANI-1x) performs on average by 1 kcal/mol better than PC9, which, however, contains about 2 orders of magnitude fewer reference structures. On the other hand, PC9 is chemically more diverse by a factor of ∼5 as quantified by the number of atom-in-molecule-based fragments (amons) it contains compared with the ANI family of databases. A quantitative measure for deficiencies is the Kullback-Leibler divergence between reference and target distributions. It is explicitly demonstrated that when certain types of bonds need to be covered in the target database (Tautobase) but are undersampled in the reference databases, the resulting predictions are poor. Examples of this include the poor performance of all databases analyzed to predict C(sp2)-C(sp2) double bonds close to heteroatoms and azoles containing N-N and N-O bonds. Analysis of the results with a Tree MAP algorithm provides deeper understanding of specific deficiencies in predicting tautomerization energies by the reference datasets due to inadequate coverage of chemical space. Capitalizing on this information can be used to either improve existing databases or generate new databases of sufficient diversity for a range of machine learning (ML) applications in chemistry.

Machine Learning Force Fields.

In recent years, the use of machine learning (ML) in computational chemistry has enabled numerous advances previously out of reach due to the computational complexity of traditional electronic-structure methods. One of the most promising applications is the construction of ML-based force fields (FFs), with the aim to narrow the gap between the accuracy of ab initio methods and the efficiency of classical FFs. The key idea is to learn the statistical relation between chemical structure and potential energy without relying on a preconceived notion of fixed chemical bonds or knowledge about the relevant interactions. Such universal ML approximations are in principle only limited by the quality and quantity of the reference data used to train them. This review gives an overview of applications of ML-FFs and the chemical insights that can be obtained from them. The core concepts underlying ML-FFs are described in detail, and a step-by-step guide for constructing and testing them from scratch is given. The text concludes with a discussion of the challenges that remain to be overcome by the next generation of ML-FFs.

Isomerization and decomposition reactions of acetaldehyde relevant to atmospheric processes from dynamics simulations on neural network-based potential energy surfaces.

Acetaldehyde (AA) isomerization [to vinylalcohol (VA)] and decomposition (into either CO + CH4 or H2 + C2H2O) are studied using a fully dimensional, reactive potential energy surface represented as a neural network (NN). The NN, trained on 432 399 reference structures from MP2/aug-cc-pVTZ calculations, has a mean absolute error of 0.0453 kcal/mol and a root mean squared error of 1.186 kcal mol-1 for a test set of 27 399 structures. For the isomerization process AA → VA, the minimum dynamical path implies that the C-H vibration and the C-C-H (with H being the transferring H-atom) and the C-C-O angles are involved to surmount the 68.2 kcal/mol barrier. Using an excess energy of 93.6 kcal/mol-the typical energy available in the solar spectrum and sufficient to excite to the first electronically excited state-to initialize the molecular dynamics, no isomerization to VA is observed on the 500 ns time scale. Only with excess energies of ∼127.6 kcal/mol (including the zero point energy of the AA molecule), isomerization occurs on the nanosecond time scale. Given that collisional quenching times under tropospheric conditions are ∼1 ns, it is concluded that formation of VA following photoexcitation of AA from actinic photons is unlikely. This also limits the relevance of this reaction pathway to be a source for formic acid.

Thermal activation of methane by MgO+: temperature dependent kinetics, reactive molecular dynamics simulations and statistical modeling.

The kinetics of MgO+ + CH4 was studied experimentally using the variable ion source, temperature adjustable selected ion flow tube (VISTA-SIFT) apparatus from 300-600 K and computationally by running and analyzing reactive atomistic simulations. Rate coefficients and product branching fractions were determined as a function of temperature. The reaction proceeded with a rate of k = 5.9 ± 1.5 × 10-10(T/300 K)-0.5±0.2 cm3 s-1. MgOH+ was the dominant product at all temperatures, but Mg+, the co-product of oxygen-atom transfer to form methanol, was observed with a product branching fraction of 0.08 ± 0.03(T/300 K)-0.8±0.7. Reactive molecular dynamics simulations using a reactive force field, as well as a neural network trained on thousands of structures yield rate coefficients about one order of magnitude lower. This underestimation of the rates is traced back to the multireference character of the transition state [MgOCH4]+. Statistical modeling of the temperature-dependent kinetics provides further insight into the reactive potential surface. The rate limiting step was found to be consistent with a four-centered activation of the C-H bond, in agreement with previous calculations. The product branching was modeled as a competition between dissociation of an insertion intermediate directly after the rate-limiting transition state, and traversing a transition state corresponding to a methyl migration leading to a Mg-CH3OH+ complex, though only if this transition state is stabilized significantly relative to the dissociated MgOH+ + CH3 product channel. An alternative, non-statistical mechanism is discussed, whereby a post-transition state bifurcation in the potential surface could allow the reaction to proceed directly from the four-centered TS to the Mg-CH3OH+ complex thereby allowing a more robust competition between the product channels.

Reactive atomistic simulations of Diels-Alder reactions: The importance of molecular rotations.

The Diels-Alder reaction between 2,3-dibromo-1,3-butadiene and maleic anhydride has been studied by means of multisurface adiabatic reactive molecular dynamics and the PhysNet neural network architecture. This system is used as a prototype to explore the concertedness, synchronicity, and possible ways of promotion of Diels-Alder reactions. Analysis of the minimum dynamic path indicates that rotational energy is crucial (∼65%) to drive the system toward the transition state in addition to collision energy (∼20%). Comparison with the reaction of butadiene and maleic anhydride shows that the presence of bromine substituents in the diene accentuates the importance of rotational excitation to promote the reaction. At the high total energies at which reactive events are recorded, the reaction is found to be direct and mostly synchronous.

Exhaustive state-to-state cross sections for reactive molecular collisions from importance sampling simulation and a neural network representation.

High-temperature, reactive gas flow is inherently nonequilibrium in terms of energy and state population distributions. Modeling such conditions is challenging even for the smallest molecular systems due to the extremely large number of accessible states and transitions between them. Here, neural networks (NNs) trained on explicitly simulated data are constructed and shown to provide quantitatively realistic descriptions which can be used in mesoscale simulation approaches such as Direct Simulation Monte Carlo to model gas flow at the hypersonic regime. As an example, the state-to-state cross sections for N(4S) + NO(2Π) → O(3P) + N2(X1Σg +) are computed from quasiclassical trajectory (QCT) simulations. By training NNs on a sparsely sampled noisy set of state-to-state cross sections, it is demonstrated that independently generated reference data are predicted with high accuracy. State-specific and total reaction rates as a function of temperature from the NN are in quantitative agreement with explicit QCT simulations and confirm earlier simulations, and the final state distributions of the vibrational and rotational energies agree as well. Thus, NNs trained on physical reference data can provide a viable alternative to computationally demanding explicit evaluation of the microscopic information at run time. This will considerably advance the ability to realistically model nonequilibrium ensembles for network-based simulations.

PhysNet: A Neural Network for Predicting Energies, Forces, Dipole Moments, and Partial Charges.

In recent years, machine learning (ML) methods have become increasingly popular in computational chemistry. After being trained on appropriate ab initio reference data, these methods allow for accurately predicting the properties of chemical systems, circumventing the need for explicitly solving the electronic Schrödinger equation. Because of their computational efficiency and scalability to large data sets, deep neural networks (DNNs) are a particularly promising ML algorithm for chemical applications. This work introduces PhysNet, a DNN architecture designed for predicting energies, forces, and dipole moments of chemical systems. PhysNet achieves state-of-the-art performance on the QM9, MD17, and ISO17 benchmarks. Further, two new data sets are generated in order to probe the performance of ML models for describing chemical reactions, long-range interactions, and condensed phase systems. It is shown that explicitly including electrostatics in energy predictions is crucial for a qualitatively correct description of the asymptotic regions of a potential energy surface (PES). PhysNet models trained on a systematically constructed set of small peptide fragments (at most eight heavy atoms) are able to generalize to considerably larger proteins like deca-alanine (Ala10): The optimized geometry of helical Ala10 predicted by PhysNet is virtually identical to ab initio results (RMSD = 0.21 Å). By running unbiased molecular dynamics (MD) simulations of Ala10 on the PhysNet-PES in gas phase, it is found that instead of a helical structure, Ala10 folds into a "wreath-shaped" configuration, which is more stable than the helical form by 0.46 kcal mol-1 according to the reference ab initio calculations.

Sampling reactive regions in phase space by following the minimum dynamic path.

Understanding mechanistic aspects of reactivity lies at the heart of chemistry. Once the potential energy surface (PES) for a system of interest is known, reactions can be studied by computational means. While the minimum energy path (MEP) between two minima of the PES can give some insight into the topological changes required for a reaction to occur, it lacks dynamical information and is an unrealistic depiction of the reactive process. For a more realistic view, molecular dynamics (MD) simulations are required. However, this usually involves generating thousands of trajectories in order to sample a few reactive events and is therefore much more computationally expensive than calculating the MEP. In this work, it is shown that a "minimum dynamic path" (MDP) can be constructed, which, contrary to the MEP, provides insight into the reaction dynamics. It is shown that the underlying concepts can be extended to directly sample reactive regions in phase space. The sampling method and the MDP are demonstrated on the well-known 2-dimensional Müller-Brown PES and for a realistic 12-dimensional reactive PES for sulfurochloridic acid, a proxy molecule used to study vibrationally induced photodissociation of sulfuric acid.

A reactive, scalable, and transferable model for molecular energies from a neural network approach based on local information.

Despite the ever-increasing computer power, accurate ab initio calculations for large systems (thousands to millions of atoms) remain infeasible. Instead, approximate empirical energy functions are used. Most current approaches are either transferable between different chemical systems, but not particularly accurate, or they are fine-tuned to a specific application. In this work, a data-driven method to construct a potential energy surface based on neural networks is presented. Since the total energy is decomposed into local atomic contributions, the evaluation is easily parallelizable and scales linearly with system size. With prediction errors below 0.5 kcal mol-1 for both unknown molecules and configurations, the method is accurate across chemical and configurational space, which is demonstrated by applying it to datasets from nonreactive and reactive molecular dynamics simulations and a diverse database of equilibrium structures. The possibility to use small molecules as reference data to predict larger structures is also explored. Since the descriptor only uses local information, high-level ab initio methods, which are computationally too expensive for large molecules, become feasible for generating the necessary reference data used to train the neural network.

Molecular Oxygen Formation in Interstellar Ices Does Not Require Tunneling.

The formation of molecular oxygen in and on amorphous ice in the interstellar medium requires oxygen diffusion to take place. Recent experiments suggest that this process involves quantum tunneling of the oxygen atoms at sufficiently low temperatures. Fitting experimental diffusion rates between 6 and 25 K to an expression that accounts for the roughness of the surface yields excellent agreement. The molecular dynamics of adsorbed oxygen is characterized by rapid intrasite dynamics, followed by intersite transitions over distances of ∼10 Å. Explicit simulations using a realistic free-energy surface for oxygen diffusion on amorphous ice down to 10 K show that quantum tunneling is not required for mobility of adsorbed oxygen. This is confirmed by comparing quantum and classical simulations using the same free-energy surface. The ratio of diffusional and desorption energy Edif/ Edes = 275/1082 ≈ 0.3 is at the lower end of typically used values but is still consistent with the assumptions made in models for interstellar chemistry.

Minimal distributed charges: Multipolar quality at the cost of point charge electrostatics.

Most empirical force fields use atom-centered point charges (PCs) to represent the electrostatic potential (ESP) around molecules. While such PC models are computationally efficient, they are unable to capture anisotropic electronic features, such as σ holes or lone pairs. These features are better described using atomic multipole (MTP) moments, which significantly improve the quality of the resulting ESP. However, the improvement comes at the expense of a considerably increased computational complexity and cost for calculating the interaction energies and forces. In the present work, a novel minimal distributed charge model (MDCM) based on off-centered point charges is presented and the quality of the resulting ESP is compared to the performance of MTPs and atom-centered PC models for several test molecules. All three models are fitted using the same algorithm based on differential evolution, which is available as a Fortran90 program from the authors upon request. We show that the MDCM is capable of approximating the reference ab initio ESP with an accuracy as good as, or better than, MTPs without the need for computationally expensive higher order multipoles. Further it is demonstrated that the MDCM is numerically stable in molecular dynamics simulations and is able to reproduce electrostatic interaction energies and thermodynamic quantities with the same accuracy as MTPs at reduced computational cost.

Quantum and quasiclassical trajectory studies of rotational relaxation in Ar-N2+ collisions.

The collision of N2+ with Ar is studied using quantum and classical methods. The dynamics was followed on a new potential energy surface based on ab initio energies computed at the UCCSD(T)-F12a/aug-cc-pVTZ level, using the correct analytical long range behaviour and a reproducing kernel representation. Comparison with multi-reference MRCI+Q calculations establish that UCCSD(T)-F12a is a sufficiently high level of theory for this problem. Results from quantum close coupling and quasiclassical trajectory calculations agree favourably with each other and the rates for inelastic collisions are lower than those from Langevin theory. This differs from previous calculations on a zero point-corrected potential energy surface (PES) and indicates that such corrections, although potentially useful, should not be applied in the present case. Despite the rather large differences between the potential energy surfaces, the computed rates are within one order of magnitude of one another which suggests that the quality of the PES is not the main reason for the remaining disagreement between computation and experiment. Also, the fraction of inelastic rotational collisions exceeds 20% in all cases irrespective of whether quantum or classical dynamics is used. Previous experimental rate coefficients for N2+(ν = 0, j = 6) colliding with Ar suggest that the rotational quantum number is largely conserved. This can not be confirmed from any of the simulations and calls for new single molecule experiments.

Toolkit for the Construction of Reproducing Kernel-Based Representations of Data: Application to Multidimensional Potential Energy Surfaces.

In the early days of computation, slow processor speeds limited the amount of data that could be generated and used for scientific purposes. In the age of big data, the limiting factor usually is the method with which large amounts of data are analyzed and useful information is extracted. A typical example from chemistry are high-level ab initio calculations for small systems, which have nowadays become feasible even if energies at many different geometries are required. Molecular dynamics simulations often require several thousand distinct trajectories to be run. Under such circumstances suitable analytical representations of potential energy surfaces (PESs) based on ab initio calculations are required to propagate the dynamics at an acceptable cost. In this work we introduce a toolkit which allows the automatic construction of multidimensional PESs from gridded ab initio data based on reproducing kernel Hilbert space (RKHS) theory. The resulting representations require no tuning of parameters and allow energy and force evaluations at ab initio quality at the same cost as empirical force fields. Although the toolkit is primarily intended for constructing multidimensional potential energy surfaces for molecular systems, it can also be used for general machine learning purposes. The software is published under the MIT license and can be downloaded, modified, and used in other projects for free.

Migration of small ligands in globins: Xe diffusion in truncated hemoglobin N.

In heme proteins, the efficient transport of ligands such as NO or O2 to the binding site is achieved via ligand migration networks. A quantitative assessment of ligand diffusion in these networks is thus essential for a better understanding of the function of these proteins. For this, Xe migration in truncated hemoglobin N (trHbN) of Mycobacterium Tuberculosis was studied using molecular dynamics simulations. Transitions between pockets of the migration network and intra-pocket relaxation occur on similar time scales (10 ps and 20 ps), consistent with low free energy barriers (1-2 kcal/mol). Depending on the pocket from where Xe enters a particular transition, the conformation of the side chains lining the transition region differs which highlights the coupling between ligand and protein degrees of freedom. Furthermore, comparison of transition probabilities shows that Xe migration in trHbN is a non-Markovian process. Memory effects arise due to protein rearrangements and coupled dynamics as Xe moves through it.

Nonadiabatic effects in electronic and nuclear dynamics.

Due to their very nature, ultrafast phenomena are often accompanied by the occurrence of nonadiabatic effects. From a theoretical perspective, the treatment of nonadiabatic processes makes it necessary to go beyond the (quasi) static picture provided by the time-independent Schrödinger equation within the Born-Oppenheimer approximation and to find ways to tackle instead the full time-dependent electronic and nuclear quantum problem. In this review, we give an overview of different nonadiabatic processes that manifest themselves in electronic and nuclear dynamics ranging from the nonadiabatic phenomena taking place during tunnel ionization of atoms in strong laser fields to the radiationless relaxation through conical intersections and the nonadiabatic coupling of vibrational modes and discuss the computational approaches that have been developed to describe such phenomena. These methods range from the full solution of the combined nuclear-electronic quantum problem to a hierarchy of semiclassical approaches and even purely classical frameworks. The power of these simulation tools is illustrated by representative applications and the direct confrontation with experimental measurements performed in the National Centre of Competence for Molecular Ultrafast Science and Technology.

Collision-induced rotational excitation in N2 (+)((2)Σg (+),v=0)-Ar: Comparison of computations and experiment.

The collisional dynamics of N2 (+)((2)Σg (+)) cations with Ar atoms is studied using quasi-classical simulations. N2 (+)-Ar is a proxy to study cooling of molecular ions and interesting in its own right for molecule-to-atom charge transfer reactions. An accurate potential energy surface (PES) is constructed from a reproducing kernel Hilbert space (RKHS) interpolation based on high-level ab initio data. The global PES including the asymptotics is fully treated within the realm of RKHS. From several ten thousand trajectories, the final state distribution of the rotational quantum number of N2 (+) after collision with Ar is determined. Contrary to the interpretation of previous experiments which indicate that up to 98% of collisions are elastic and conserve the quantum state, the present simulations find a considerably larger number of inelastic collisions which supports more recent findings.