Search

We introduce a novel enhanced sampling approach named OPES flooding for calculating the kinetics of rare events from atomistic molecular dynamics simulation. This method is derived from the On-the-fly-Probability-Enhanced-Sampling (OPES) approach [Invernizzi and Parrinello, JPC Lett. 2020], which has been recently developed for calculating converged free energy surfaces for complex systems. In this paper, we describe the theoretical details of the OPES flooding technique and demonstrate the application on three systems of increasing complexity: barrier crossing in a two-dimensional double well potential, conformational transition in the alanine dipeptide in gas phase, and the folding and unfolding of the chignolin polypeptide in aqueous environment. From extensive tests, we show that the calculation of accurate kinetics not only requires the transition state to be bias-free, but the amount of bias deposited should also not exceed the effective barrier height measured along the chosen collective variables. In this vein, the possibility of computing rates from biasing suboptimal order parameters has also been explored. Furthermore, we describe the choice of optimum parameter combinations for obtaining accurate results from limited computational effort.

Enhanced sampling techniques have revolutionised molecular dynamics (MD) simulations, enabling the study of rare events and the calculation of free energy differences in complex systems. One of the main families of enhanced sampling techniques uses physical degrees of freedom called collective variables (CVs) to accelerate a system’s dynamics and recover the original system’s statistics. However, encoding all the relevant degrees of freedom in a limited number of CVs is challenging, particularly in large biophysical systems. Another category of techniques, such as parallel tempering, simulates multiple replicas of the system in parallel, without requiring CVs. However, these methods may explore less relevant high-energy portions of the phase space and become computationally expensive for large systems. To overcome the limitations of both approaches, we propose a replica exchange method called OneOPES that combines the power of multi-replica simulations and CV-based enhanced sampling. This method efficiently accelerates the crossing of both enthalpic and entropic barriers without the need for optimal CVs definition, parameters fine tuning nor the use of a large number of replicas, as demonstrated by its successful applications to protein-ligand binding and protein folding benchmark systems. Our approach shows promise as a new direction in the development of enhanced sampling techniques for molecular dynamics simulations, providing an efficient and robust framework for the study of complex and unexplored problems.

The main protease from SARS-CoV-2 is a homodimer. Yet, a recent 0.1 ms long molecular dynamics simulation shows that it readily undergoes a symmetry breaking event on passing from the solid state to the aqueous solution. As a result, the subunits present distinct conformations of the binding pocket. By analysing this long time simulation, here we uncover a previously unrecognised role of water molecules in triggering the transition. Interestingly, each subunit presents a different collection of long-lived water molecules. Enhanced sampling methods performed here, along with machine learning approaches, further establish that the transition to the asymmetric state is essentially irreversible.

One of the main applications of atomistic computer simulations is the calculation of ligand binding free energies. The accuracy of these calculations depends on the force field quality and on the thoroughness of configuration sampling. Sampling is an obstacle in simulations due to the frequent appearance of kinetic bottlenecks in the free energy landscape. Very often this difficulty is circumvented by enhanced sampling techniques. Typically, these techniques depend on the introduction of appropriate collective variables that are meant to capture the system’s degrees of freedom. In ligand binding, water has long been known to play a key role, but its complex behaviour has proven difficult to fully capture. In this paper we combine machine learning with physical intuition to build a non-local and highly efficient water-describing collective variable. We use it to study a set of host-guest systems from the SAMPL5 challenge. We obtain highly accurate binding free energies and good agreement with experiments. The role of water during the binding process is then analysed in some detail. © 2021, The Author(s)., Nature Communications, 12 (1), ISSN:2041-1723

Over the last few decades enhanced sampling methods have made great strides. Here, we exploit this progress and propose a modular workflow for blind reaction discovery and characterization of reaction paths. Central to our strategy is the use of the recently developed explore variant of the on-the-fly probability enhanced sampling method. Like metadynamics, this method is based on the identification of appropriate collective variables. Our first step is the discovery of new chemical reactions and it is performed biasing a one dimensional collective variable derived from spectral graph theory. Once new reaction pathways are detected, we construct ad-hoc tailored neural-network based collective variables to improve sampling of specific reactions and finally we refine the results using free energy perturbation theory. Our workflow has been successfully applied to both intramolecular and intermolecular reactions. Without any chemical hypothesis, we discovered several possible products, computed the free energy surface at semiempirical level, and finally refined it with a more accurate Hamiltonian. Our workflow requires minimal user input, and thanks to its modularity and flexibility, can extend the scope of ab initio molecular dynamics for the exploration and characterization of reaction space.

Disentangling the mechanistic details of a chemical reaction pathway is a hard problem that often requires a considerable amount of chemical intuition and a component of luck. Experiments struggle in observing short-life metastable intermediates, while computer simulations often rely upon a good initial guess. In this work, we propose a method that, from the simulations of a reactant and a product state, searches for reaction mechanisms connecting the two by exploring the configuration space through metadynamics, a well-known enhanced molecular dynamics method. The key quantity underlying this search is based on the use of an approach called harmonic linear discriminant analysis which allows a systematic construction of collective variables. Given the reactant and product states, we choose a set of descriptors capable of discriminating between the two states. In order to not prejudge the results, generic descriptors are introduced. The fluctuations of the descriptors in the two states are used to construct collective variables. We use metadynamics in an exploratory mode to discover the intermediates and the transition states that lead from reactant to product. The search is at first conducted at a low theory level. The calculation is then refined, and the energy of the intermediates and transition states discovered during metadynamics is computed again using a higher level of theory. The method's aim is to offer a simple reaction mechanism search procedure that helps in saving time and is able to find unexpected mechanisms that defy well established chemical paradigms. We apply it to two reactions, showing that a high level of complexity can be hidden even in seemingly trivial and small systems. The method can be applied to larger systems, such as reactions in solution or catalysis.

The phenomenon of solidification of a substance from its liquid phase is of the greatest practical and theoretical importance, and atomistic simulations can provide precious information towards its understanding and control. Unfortunately, the time scale for crystallization is much larger than what can be explored in standard simulations. Enhanced sampling methods can overcome this time scale hurdle. Here we employ the on-the-fly probability enhanced sampling method that is a recent evolution of metadynamics. This method, like many others, relies on the definition of appropriate collective variables able to capture the slow degrees of freedom. To this effect we introduce collective coordinates of general applicability to crystallization simulations. They are based on the peaks of the three-dimensional structure factor that are combined non-linearly via the Deep Linear Discriminant Analysis machine learning method. We apply the method to the study of crystallization of a multicomponent system, Sodium Chloride and a molecular system, Carbon Dioxide., Comment: 23 pages, 3 figures

Designing an appropriate set of collective variables is crucial to the success of several enhanced sampling methods. Here we focus on how to obtain such variables from information limited to the metastable states. We characterize these states by a large set of descriptors and employ neural networks to compress this information in a lower-dimensional space, using Fisher's linear discriminant as an objective function to maximize the discriminative power of the network. We test this method on alanine dipeptide, using the non-linearly separable dataset composed by atomic distances. We then study an intermolecular aldol reaction characterized by a concerted mechanism. The resulting variables are able to promote sampling by drawing non-linear paths in the physical space connecting the fluctuations between metastable basins. Lastly, we interpret the behavior of the neural network by studying its relation to the physical variables. Through the identification of its most relevant features, we are able to gain chemical insight into the process.

Trapped bosons exhibit fundamental physical phenomena and are at the core of emerging quantum technologies. We present a method for simulating bosons using path integral molecular dynamics. The main difficulty in performing such simulations is enumerating all ring-polymer configurations, which arise due to permutations of identical particles. We show that the potential and forces at each time step can be evaluated by using a recurrence relation which avoids enumerating all permutations, while providing the correct thermal expectation values. The resulting algorithm scales cubically with system size. The method is tested and applied to bosons in a 2-dimensional (2D) trap and agrees with analytical results and numerical diagonalization of the many-body Hamiltonian. An analysis of the role of exchange effects at different temperatures, through the relative probability of different ring-polymer configurations, is also presented.

This chapter presents a new methodological development that represents an improvement over ECEID. Through the use of a canonical transformation on the ECEID Hamiltonian, it is possible to include the dynamics of an extra classical degree of freedom in the model, a semiclassical oscillator position. Physically, we can think of it as a representation of the oscillator centroid, but it enters the derivation just as a time-dependent reference position.

Thermalization between electronic and vibrational degrees of freedom arises in a range of physical situations spanning widely different time and length scales.

Heat transport can be seen microscopically as a sequence of collisions between energy carriers. The concepts of mean free path and mean free time play a key role in heat transport problems, the former being the average distance that a carrier covers between two collisions and the latter the average time, taken over a large number of collision events.

In this chapter we present the derivation of the method that we developed and used in our simulations: Effective Correlated Electron-Ion Dynamics (ECEID). We start from the model Hamiltonian and, after some approximations, we obtain a set of equations of motion (EOM) for electronic and vibrational quantities. We also discuss total energy conservation, the Open-Boundaries implementation and the Many-Body to One-Body projection of the EOM. We then outline the implementation of the method in a code and perform scaling tests with a varying number of oscillators and electronic sites. In Chap. 9 we will present a recent and more general reformulation of the method.

In this chapter we discuss explicit methods that tackle electron-phonon problems at the microscopic level. To capture the exact dynamics of coupled systems of interacting electrons and phonon modes is an insurmountable task for computer simulations. The tremendous memory footprint of the full problem makes simulations unfeasible even for systems including just a few electrons. Approximations must be employed and some information has to be discarded. Depending on the problem under study, it is common to focus only on certain aspects of the full time-dependent coupled problem.

When high-energy radiation penetrates living cells, it ionizes molecules along its path and can cause cell death by damaging DNA. Radiation events involve a sequence of processes that ultimately require a clear microscopic understanding [1]. Only about one third of the cellular damage is produced by direct interaction of the ionizing radiation with DNA, while the rest is due to secondary species, produced in the first hundreds to thousands of femtoseconds following the primary irradiation of the system [2].

Irradiation of biological matter triggers a cascade of secondary particles that interact with their surroundings, resulting in damage. Low-energy electrons are one of the main secondary species and electron-phonon interaction plays a fundamental role in their dynamics. We have developed a method to capture the electron-phonon inelastic energy exchange in real time and have used it to inject electrons into a simple system that models a biological environment, a water chain. We simulated both an incoming electron pulse and a steady stream of electrons and found that electrons with energies just outside bands of excited molecular states can enter the chain through phonon emission or absorption. Furthermore, this phonon-assisted dynamical behaviour shows great sensitivity to the vibrational temperature, highlighting a crucial controlling factor for the injection and propagation of electrons in water.

We present a quantum simulation method that follows the dynamics of out-of-equilibrium many-body systems of electrons and oscillators in real time. Its cost is linear in the number of oscillators and it can probe time scales from attoseconds to hundreds of picoseconds. Contrary to Ehrenfest dynamics, it can thermalize starting from a variety of initial conditions, including electronic population inversion. While an electronic temperature can be defined in terms of a nonequilibrium entropy, a Fermi-Dirac distribution in general emerges only after thermalization. These results can be used to construct a kinetic model of electron-phonon equilibration based on the explicit quantum dynamics.

The discovery of effective hydrogen storage materials is fundamental for the progress of a clean energy economy. Ammonia borane ($\mathrm{H_3BNH_3}$) has attracted great interest as a promising candidate but the reaction path that leads from its solid phase to hydrogen release is not yet fully understood. To address the need for insights in the atomistic details of such a complex solid state process, in this work we use \textit{ab-initio} molecular dynamics and metadynamics to study the early stages of AB dehydrogenation. We show that the formation of ammonia diborane ($\mathrm{H_3NBH_2(}$$\mu$$\mathrm{-H)BH_3}$) leads to the release of $\mathrm{NH_4^+}$, which in turn triggers an autocatalytic $\mathrm{H_2}$ production cycle. Our calculations provide a model for how complex solid state reactions can be theoretically investigated and rely upon the presence of multiple ammonia borane molecules, as substantiated by standard quantum-mechanical simulations on a cluster.