Search

This letter illustrates the opinion of the molecular dynamics (MD) community on the need to adopt a new FAIR paradigm for the use of molecular simulations. It highlights the necessity of a collaborative effort to create, establish, and sustain a database that allows findability, accessibility, interoperability, and reusability of molecular dynamics simulation data. Such a development would democratize the field and significantly improve the impact of MD simulations on life science research. This will transform our working paradigm, pushing the field to a new frontier. We invite you to support our initiative at the MDDB community (https://mddbr.eu/community/)

We investigate learning the eigenfunctions of evolution operators for time-reversal invariant stochastic processes, a prime example being the Langevin equation used in molecular dynamics. Many physical or chemical processes described by this equation involve transitions between metastable states separated by high potential barriers that can hardly be crossed during a simulation. To overcome this bottleneck, data are collected via biased simulations that explore the state space more rapidly. We propose a framework for learning from biased simulations rooted in the infinitesimal generator of the process and the associated resolvent operator. We contrast our approach to more common ones based on the transfer operator, showing that it can provably learn the spectral properties of the unbiased system from biased data. In experiments, we highlight the advantages of our method over transfer operator approaches and recent developments based on generator learning, demonstrating its effectiveness in estimating eigenfunctions and eigenvalues. Importantly, we show that even with datasets containing only a few relevant transitions due to sub-optimal biasing, our approach recovers relevant information about the transition mechanism.

Determining the kinetic bottlenecks that make transitions between metastable states difficult is key to understanding important physical problems like crystallization, chemical reactions, or protein folding. In all these phenomena, the system spends a considerable amount of time in one metastable state before making a rare but important transition to a new state. The rarity of these events makes their direct simulation challenging, if not impossible. We propose a method to explore the distribution of configurations that the system passes as it translocates from one metastable basin to another. We shall refer to this set of configurations as the transition state ensemble. We base our method on the committor function and the variational principle to which it obeys. We find the minimum of the variational principle via a self-consistent procedure that does not require any input besides the knowledge of the initial and final state. Right from the start, our procedure focuses on sampling the transition state ensemble and allows harnessing a large number of such configurations. With the help of the variational principle, we perform a detailed analysis of the transition state ensemble, ranking quantitatively the degrees of freedom mostly involved in the transition and opening the way for a systematic approach for the interpretation of simulation results and the construction of collective variables.

The understanding of prototropic tautomerism in water and the characterization of solvent effects on protomeric equilibrium pose significant challenges. Using molecular dynamics simulations based on state-of-the-art deep learning potential and enhanced sampling methods, we provide a comprehensive description of all configurational transformations in glycine solvated in water and determine accurate free energy profiles of these processes. We observe that the tautomerism between the neutral and zwitterionic forms of solvated glycine can occur by both intramolecular proton transfer in glycine and intermolecular proton transfer in the contact ion pair (anionic glycine and hydronium ion) or the separated ion pair (cationic glycine and hydroxide ion).

A variety of enhanced sampling methods predict multidimensional free energy landscapes associated with biological and other molecular processes as a function of a few selected collective variables (CVs). The accuracy of these methods is crucially dependent on the ability of the chosen CVs to capture the relevant slow degrees of freedom of the system. For complex processes, finding such CVs is the real challenge. Machine learning (ML) CVs offer, in principle, a solution to handle this problem. However, these methods rely on the availability of high-quality datasets -- ideally incorporating information about physical pathways and transition states -- which are difficult to access, therefore greatly limiting their domain of application. Here, we demonstrate how these datasets can be generated by means of enhanced sampling simulations in trajectory space via the metadynamics of paths [arXiv:2002.09281] algorithm. The approach is expected to provide a general and efficient way to generate efficient ML-based CVs for the fast prediction of free energy landscapes in enhanced sampling simulations. We demonstrate our approach with two numerical examples, a two-dimensional model potential and the isomerization of alanine dipeptide, using deep targeted discriminant analysis as our ML-based CV of choice.

The anomalous $\lambda$-transition of liquid sulfur, which is supposed to be related to the transformation of eight-membered sulfur rings into long polymeric chains, has attracted considerable attention. However, a detailed description of the underlying dynamical polymerization process is still missing. Here, we study the structures and the mechanism of the polymerization processes of liquid sulfur across the $\lambda$-transition as well as its reverse process of formation of the rings. We do so by performing ab-initio-quality molecular dynamics simulations thanks to a combination of machine learning potentials and state-of-the-art enhanced sampling techniques. With our approach, we obtain structural results that are in good agreement with the experiments and we report precious dynamical insights into the mechanisms involved in the process.

Studying the pathways of ligand-receptor binding is essential to understand the mechanism of target recognition by small molecules. The binding free energy and kinetics of protein-ligand complexes can be computed using molecular dynamics (MD) simulations, often in quantitative agreement with experiments. However, only a qualitative picture of the ligand binding/unbinding paths can be obtained through a conventional analysis of the MD trajectories. Besides, the higher degree of manual effort involved in analyzing pathways limits its applicability in large-scale drug discovery. Here we address this limitation by introducing an automated approach for analyzing molecular transition paths with a particular focus on protein-ligand dissociation. Our method is based on the dynamic time-warping (DTW) algorithm, originally designed for speech recognition. We accurately classified molecular trajectories using a very generic descriptor set of contacts or distances. Our approach outperforms manual classification by distinguishing between parallel dissociation channels, within the pathways identified by visual inspection. Most notably, we could compute exit-path-specific ligand-dissociation kinetics. The unbinding timescale along the fastest path agrees with the experimental residence time, providing a physical interpretation to our entirely data-driven protocol. In combination with appropriate enhanced sampling algorithms, this technique can be used for the initial exploration of ligand-dissociation pathways as well as for calculating path-specific thermodynamic and kinetic properties.

This work examines the conformational ensemble involved in $\beta$-hairpin folding by means of advanced molecular dynamics simulations and dimensionality reduction. A fully atomistic description of the protein and the surrounding solvent molecules is used and this complex energy landscape is sampled by means of parallel tempering metadynamics simulations. The ensemble of configurations explored is analysed using the recently proposed sketch-map algorithm. Further simulations allow us to probe how mutations affect the structures adopted by this protein. We find that many of the configurations adopted by a mutant are the same as those adopted by the wild type protein. Furthermore, certain mutations destabilize secondary structure containing configurations by preventing the formation of hydrogen bonds or by promoting the formation of new intramolecular contacts. Our analysis demonstrates that machine-learning techniques can be used to study the energy landscapes of complex molecules and that the visualizations that are generated in this way provide a natural basis for examining how the stabilities of particular configurations of the molecule are affected by factors such as temperature or structural mutations., Comment: Pre-print of doi:10.1021/ct500950z, reproduced with permission from the ACS

Interatomic potentials learned using machine learning methods have been successfully applied to atomistic simulations. However, accurate models require large training datasets, while generating reference calculations is computationally demanding. To bypass this difficulty, we propose a transfer learning algorithm that leverages the ability of graph neural networks (GNNs) to represent chemical environments together with kernel mean embeddings. We extract a feature map from GNNs pre-trained on the OC20 dataset and use it to learn the potential energy surface from system-specific datasets of catalytic processes. Our method is further enhanced by incorporating into the kernel the chemical species information, resulting in improved performance and interpretability. We test our approach on a series of realistic datasets of increasing complexity, showing excellent generalization and transferability performance, and improving on methods that rely on GNNs or ridge regression alone, as well as similar fine-tuning approaches., Comment: 20 pages, 4 figures, 7 tables, published in NeurIPS 2023

Identifying a reduced set of collective variables is critical for understanding atomistic simulations and accelerating them through enhanced sampling techniques. Recently, several methods have been proposed to learn these variables directly from atomistic data. Depending on the type of data available, the learning process can be framed as dimensionality reduction, classification of metastable states or identification of slow modes. Here we present $\texttt{mlcolvar}$, a Python library that simplifies the construction of these variables and their use in the context of enhanced sampling through a contributed interface to the PLUMED software. The library is organized modularly to facilitate the extension and cross-contamination of these methodologies. In this spirit, we developed a general multi-task learning framework in which multiple objective functions and data from different simulations can be combined to improve the collective variables. The library's versatility is demonstrated through simple examples that are prototypical of realistic scenarios.

The phase change compound Ge$_2$Sb$_2$Te$_5$ (GST225) is exploited in advanced non-volatile electronic memories and in neuromorphic devices which both rely on a fast and reversible transition between the crystalline and amorphous phases induced by Joule heating. The crystallization kinetics of GST225 is a key functional feature for the operation of these devices. We report here on the development of a machine-learned interatomic potential for GST225 that allowed us to perform large scale molecular dynamics simulations (over 10000 atoms for over 100 ns) to uncover the details of the crystallization kinetics in a wide range of temperatures of interest for the programming of the devices. The potential is obtained by fitting with a deep neural network (NN) scheme a large quantum-mechanical database generated within Density Functional Theory. The availability of a highly efficient and yet highly accurate NN potential opens the possibility to simulate phase change materials at the length and time scales of the real devices.

The study of the rare transitions that take place between long lived metastable states is a major challenge in molecular dynamics simulations. Many of the methods suggested to address this problem rely on the identification of the slow modes of the system which are referred to as collective variables. Recently machine learning methods have been used to learn the collective variables as functions of a large number of physical descriptors. Among many such methods Deep Targeted Discriminant Analysis has proven to be useful. This collective variable is built from data harvested in short unbiased simulation in the two basins. Here we enrich the set of data on which the Deep Targeted Discriminant Analysis collective variable is built by adding data coming from the transition path ensemble. These are collected from a number of reactive trajectories obtained using the On-the-fly Probability Enhanced Sampling Flooding method. The collective variables thus trained, lead to a more accurate sampling and faster convergence. The performance of these new collective variables is tested on a number of representative examples.

We present a new method to compute free energies at a quantum mechanical (QM) level of theory from molecular simulations using cheap reference potential energy functions, such as force fields. To overcome the poor overlap between the reference and target distributions, we generalize targeted free energy perturbation (TFEP) to employ multiple configuration maps. While TFEP maps have been obtained before from an expensive training of a normalizing flow neural network (NN), our multimap estimator allows us to use the same set of QM calculations to both optimize the maps and estimate the free energy, thus removing almost completely the overhead due to training. A multimap extension of the multistate Bennett acceptance ratio estimator is also derived for cases where samples from two or more states are available. Furthermore, we propose a one-epoch learning policy that can be used to efficiently avoid overfitting when computing the loss function is expensive compared to generating data. Finally, we show how our multimap approach can be combined with enhanced sampling strategies to overcome the pervasive problem of poor convergence due to slow degrees of freedom. We test our method on the HiPen dataset of drug-like molecules and fragments, and we show that it can accelerate the calculation of the free energy difference of switching from a force field to a DFTB3 potential by about 3 orders of magnitude compared to standard FEP and by a factor of about 8 compared to previously published nonequilibrium calculations., Comment: Main Text: 12 pages, 5 figures, 7 equations. Supplemental Material: 17 pages, 5 figures, 22 equations

A key step in the development of a new drug is the design of drug-excipient complexes that lead to an optimal drug release kinetics. Computational chemistry, specifically enhanced sampling molecular dynamics methods, can play a key role in this context by minimizing the need for expensive experiments and reducing cost and time. Here, we show that recent advances in enhanced sampling methodologies can be brought to fruition in this area. We demonstrate the potential of these methodologies by simulating the release kinetics of the complex praziquantel-montmorillonite in water. Praziquantel finds promising applications in the treatment of schistosomiasis, but its biopharmaceutical profile needs to be improved and a cheap material like the montmorillonite clay would be a very convenient excipient. We simulate the drug release both from surface and interlayer space and find that the time it takes to the drug to diffuse from the interlayer space to the solution controls the release time.

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.

Present-day atomistic simulations generate long trajectories of ever more complex systems. Analyzing these data, discovering metastable states, and uncovering their nature is becoming increasingly challenging. In this paper, we first use the variational approach to conformation dynamics to discover the slowest dynamical modes of the simulations. This allows the different metastable states of the system to be located and organized hierarchically. The physical descriptors that characterize metastable states are discovered by means of a machine learning method. We show in the cases of two proteins, Chignolin and Bovine Pancreatic Trypsin Inhibitor, how such analysis can be effortlessly performed in a matter of seconds. Another strength of our approach is that it can be applied to the analysis of both unbiased and biased simulations., Comment: Main text: 10 pages, 4 figures. Supplementary Info: 4 pages, 5, figures

The process of ligand-protein unbinding is crucial in biophysics. Water is an essential part of any biological system and yet, many aspects of its role remain elusive. Here, we simulate with state-of-the-art enhanced sampling techniques the binding of Benzamidine to Trypsin which is a much studied and paradigmatic ligand-protein system. We use machine learning methods to determine efficient collective coordinates for the complex non-local network of water. These coordinates are used to perform On-the-fly Probability Enhanced Sampling simulations, which we adapt to calculate also the ligand residence time. Our results, both static and dynamic, are in good agreement with experiments. We find that the presence of a water molecule located at the bottom of the binding pocket allows via a network of hydrogen bonds the ligand to be released into the solution. On a finer scale, even when unbinding is allowed, another water molecule further modulates the exit time.

We combine molecular dynamics simulations with experiments to estimate solubilities of organic salts in complex growth environments. We predict the solubility by simulations of the growth and dissolution of ions at the crystal surface kink sites at different solution concentrations. Thereby, the solubility is identified as the solution's salt concentration, where the energy of the ion pair dissolved in solution equals the energy of the ion pair crystallized at the kink sites. The simulation methodology is demonstrated for the case of anhydrous sodium acetate crystallized from various solvent-antisolvent mixtures. To validate the predicted solubilities, we have measured the solubilities of sodium acetate in-house, using an experimental setup and measurement protocol that guarantees moisture-free conditions, which is key for a hygroscopic compound like sodium acetate. We observe excellent agreement between the experimental and the computationally evaluated solubilities for sodium acetate in different solvent-antisolvent mixtures. Given the agreement and the rich data the simulations produce, we can use them to complement experimental tasks which in turn will reduce time and capital in the design of complicated industrial crystallization processes of organic salts.

In adaptive-bias enhanced sampling methods, a bias potential is added to the system to drive transitions between metastable states. The bias potential is a function of a few collective variables and is gradually modified according to the underlying free energy surface. We show that when the collective variables are suboptimal, there is an exploration-convergence tradeoff, and one must choose between a quickly converging bias that will lead to fewer transitions, or a slower to converge bias that can explore the phase space more efficiently but might require a much longer time to produce an accurate free energy estimate. The recently proposed On-the-fly Probability Enhanced Sampling (OPES) method focuses on fast convergence, but there are cases where fast exploration is preferred instead. For this reason, we introduce a new variant of the OPES method that focuses on quickly escaping metastable states, at the expense of convergence speed. We illustrate the benefits of this approach on prototypical systems and show that it outperforms the popular metadynamics method.

The long timescale associated with ligand residence times renders their computation challenging. Therefore, the influence of factors like solvation and steric hindrance on residence times are not fully understood. Here, we demonstrate in a set of model host-guest systems that the recently developed Gaussian Mixture Based Enhanced Sampling allows residence times to be computed and enables understanding their unbinding mechanism. We observe that guest unbinding often proceeds via a series of intermediate states that can be labelled by the number of water molecules present in the binding cavity. And in several cases the residence time is correlated to the water trapping times in the cavity., Comment: 11 pages, including Supplementary Information, 5 figures

The determination of efficient collective variables is crucial to the success of many enhanced sampling methods. As inspired by previous discrimination approaches, we first collect a set of data from the different metastable basins. The data are then projected with the help of a neural network into a low-dimensional manifold in which data from different basins are well discriminated. This is here guaranteed by imposing that the projected data follows a preassigned distribution. The collective variables thus obtained lead to an efficient sampling and often allow reducing the number of collective variables in a multi-basin scenario. We first check the validity of the method in two-state systems. We then move to multi-step chemical processes. In the latter case, at variance with previous approaches, one single collective variable suffices, leading not only to computational efficiency but to a very clear representation of the reaction free energy profile.

The development of enhanced sampling methods has greatly extended the scope of atomistic simulations, allowing long-time phenomena to be studied with accessible computational resources. Many such methods rely on the identification of an appropriate set of collective variables. These are meant to describe the system's modes that most slowly approach equilibrium. Once identified, the equilibration of these modes is accelerated by the enhanced sampling method of choice. An attractive way of determining the collective variables is to relate them to the eigenfunctions and eigenvalues of the transfer operator. Unfortunately, this requires knowing the long-term dynamics of the system beforehand, which is generally not available. However, we have recently shown that it is indeed possible to determine efficient collective variables starting from biased simulations. In this paper, we bring the power of machine learning and the efficiency of the recently developed on-the-fly probability enhanced sampling method to bear on this approach. The result is a powerful and robust algorithm that, given an initial enhanced sampling simulation performed with trial collective variables or generalized ensembles, extracts transfer operator eigenfunctions using a neural network ansatz and then accelerates them to promote sampling of rare events. To illustrate the generality of this approach we apply it to several systems, ranging from the conformational transition of a small molecule to the folding of a mini-protein and the study of materials crystallization.

We present an approach that extends the theory of targeted free energy perturbation (TFEP) to calculate free energy differences and free energy surfaces at an accurate quantum mechanical level of theory from a cheaper reference potential. The convergence is accelerated by a mapping function that increases the overlap between the target and the reference distributions. Building on recent work, we show that this map can be learned with a normalizing flow neural network, without requiring simulations with the expensive target potential but only a small number of single-point calculations, and, crucially, avoiding the systematic error that was found previously. We validate the method by numerically evaluating the free energy difference in a system with a double-well potential and by describing the free energy landscape of a simple chemical reaction in the gas phase., Comment: Main Text: 7 pages, 2 figures, 18 equations. Supplemental Material: 4 pages, 26 equations

We present a modified version of the nudged elastic band (NEB) algorithm to find minimum energy paths con-necting two known configurations. We show that replacing the harmonic band-energy term with a discretized version of the Onsager-Machlup action leads to a NEB algorithm with adaptive spring lengths that automatically increase the resolution of the minimum energy path around the saddle point of the potential energy surface. The method has the same computational cost per optimization step of the standard NEB algorithm and does not introduce additional parameters. We present applications to the isomerization of alanine dipeptide, the elimination of hydrogen from ethane and the healing of a 5-77-5 defect in graphene., Comment: Main text: 9 pages, 7 Figures, 4 tables; Supplementary Materials: 9 pages, 15 Figures

The study of liquid-liquid phase transition has attracted considerable attention. One interesting example of such phenomenon is phosphorus for which the existence a first-order phase transition between a low density insulating molecular phase and a conducting polymeric phase has been experimentally established. In this paper, we model this transition by an ab-initio quality molecular dynamics simulation and explore a large portion of the liquid section of the phase diagram. We draw the liquid-liquid coexistence curve and determine that it terminates into a second-order critical point. Close to the critical point, large coupled structure and electronic structure fluctuations are observed., Comment: 6 pages and 5 figures

We present a molecular dynamics simulation method for the computation of the solubility of organic crystals in solution. The solubility is calculated based on the equilibrium free energy difference between the solvated solute and its crystallized state at the crystal surface kink site. In order to efficiently sample the growth and dissolution process, we have carried out well-tempered Metadynamics simulations with a collective variable that captures the slow degrees of freedom, namely the solute diffusion to and adsorption at the kink site together with the desolvation of the kink site. Simulations were performed at different solution concentrations using constant chemical potential molecular dynamics and the solubility was identified at the concentration at which the free energy values between the grown and dissolved kink states were equal. The effectiveness of this method is demonstrated by its success in reproducing the experimental trends of solubility of urea and naphthalene in a variety of solvents.

Supersolid is a mysterious and puzzling state of matter whose possible existence has stirred a vigorous debate among physicists for over 60 years. Its elusive nature stems from the coexistence of two seemingly contradicting properties, long-range order and superfluidity. We report computational evidence of a supersolid phase of deuterium under high pressure ($p >800$ GPa) and low temperature (T $<$ 1.0 K). In our simulations, that are based on bosonic path integral molecular dynamics, we observe a highly concerted exchange of atoms while the system preserves its crystalline order. The exchange processes are favoured by the soft core interactions between deuterium atoms that form a densely packed metallic solid. At the zero temperature limit, Bose-Einstein condensation is observed as the permutation probability of $N$ deuterium atoms approaches $1/N$ with a finite superfluid fraction. Our study provides concrete evidence for the existence of a supersolid phase in high-pressure deuterium and could provide insights on the future investigation of supersolid phases in real materials., Comment: 5 pages, 3 figures

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

The study of chemical reactions in aqueous media is very important for its implications in several fields of science, from biology to industrial processes. Modelling these reactions is however difficult when water directly participates in the reaction. Since it requires a fully quantum mechanical description of the system, $\textit{ab-initio}$ molecular dynamics is the ideal candidate to shed light on these processes. However, its scope is limited by a high computational cost. A popular alternative is to perform molecular dynamics simulations powered by machine learning potentials, trained on an extensive set of quantum mechanical calculations. Doing so reliably for reactive processes is difficult because it requires including very many intermediate and transition state configurations. In this study, we used an active learning procedure accelerated by enhanced sampling to harvest such structures and to build a neural-network potential to study the urea decomposition process in water. This allowed us to obtain the free energy profiles of this important reaction in a wide range of temperatures, to discover a number of novel metastable states and to improve the accuracy of the kinetic rates calculations. Furthermore, we found that the formation of the zwitterionic intermediate has the same probability of occurring via an acidic or a basic pathway, which could be the cause of the insensitivity of reaction rates to the pH solution.

The sampling problem lies at the heart of atomistic simulations and over the years many different enhanced sampling methods have been suggested towards its solution. These methods are often grouped into two broad families. On the one hand methods such as umbrella sampling and metadynamics that build a bias potential based on few order parameters or collective variables. On the other hand, tempering methods such as replica exchange that combine different thermodynamic ensembles in one single expanded ensemble. We instead adopt a unifying perspective, focusing on the target probability distribution sampled by the different methods. This allows us to introduce a new class of collective-variables-based bias potentials that can be used to sample any of the expanded ensembles normally sampled via replica exchange. We also provide a practical implementation, by properly adapting the iterative scheme of the recently developed on-the-fly probability enhanced sampling method [Invernizzi and Parrinello, J. Phys. Chem. Lett. 11.7 (2020)], which was originally introduced for metadynamics-like sampling. The resulting method is very general and can be used to achieve different types of enhanced sampling. It is also reliable and simple to use, since it presents only few and robust external parameters and has a straightforward reweighting scheme. Furthermore, it can be used with any number of parallel replicas. We show the versatility of our approach with applications to multicanonical and multithermal-multibaric simulations, thermodynamic integration, umbrella sampling, and combinations thereof.

One of the main applications of atomistic computer simulations is the calculation of ligand binding energies. The accuracy of these calculations depends on the force field quality and on the thoroughness of configuration sampling. Sampling is an obstacle in modern 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 of host-guest systems from the SAMPL5 challenge. We obtain highly accurate binding energies and good agreement with experiments. The role of water during the binding process is then analysed in some detail.

Accurate prediction of a gas solubility in a liquid is crucial in many areas of chemistry, and a detailed understanding of the molecular mechanism of the gas solvation continues to be an active area of research. Here, we extend the idea of constant chemical potential molecular dynamics (C{\mu}MD) approach to the calculation of the gas solubility in the liquid under constant gas chemical potential conditions. As a representative example, we utilize this method to calculate the isothermal solubility of carbon dioxide in water. Additionally, we provide microscopic insight into the mechanism of solvation that preferentially occurs in areas of the surface where the hydrogen network is broken.

We present an ab initio molecular dynamics (MD) investigation of the tautomeric equilibrium for aqueous solutions of glycine and acetone at realistic experimental conditions. Metadynamics is used to accelerate proton migration among tautomeric centers. Due to the formation of complex water-ion structures involved the proton dynamics in the aqueous environment, standard enhanced sampling approaches may face severe limitations in providing a general description of the phenomenon. Recently, we developed a set of Collective Variables (CVs) designed to study protons transfer reactions in complex condensed systems [Grifoni et al. PNAS, 2019, 116(10), 4054-4057]. In this work we applied this approach to study proton dissociation dynamics leading to tautomeric interconversion of biologically and chemically relevant prototypical systems, namely glycine and acetone in water. Although relatively simple from a chemical point of view, the results show that even for these small systems complex reaction pathways and non-trivial conversion dynamics are observed. The generality of our method allows obtaining these results without providing any prior information on the dissociation dynamics but only the atomic species that can exchange protons in the process. Our results agree with literature estimates and demonstrate the general applicability of this method in the study of tautomeric reactions.

CO2 has attracted considerable attention in the recent years due to its role in the greenhouse effect and environmental management. While its reaction with water has been studied extensively, the same cannot be said for reactivity in supercritical CO2 phase, where the conjugate acid/base equilibria proceed through different mechanisms and activation barriers. In spite of the apparent simplicity of the CO2 + H2O reaction, the collective effect of different environments has drastic influence on the free energy profile. Enhanced sampling techniques and well-tailored collective variables provide a detailed picture of the enthalpic and entropic drivers underscoring the differences in the formation mechanism of carbonic acid in the gas, aqueous and supercritical CO2 phases.

We present a method for performing path integral molecular dynamics (PIMD) simulations for fermions and address its sign problem. PIMD simulations are widely used for studying many-body quantum systems at thermal equilibrium. However, they assume that the particles are distinguishable and neglect bosonic and fermionic exchange effects. Interacting fermions play a key role in many chemical and physical systems, such as electrons in quantum dots and ultracold trapped atoms. A direct sampling of the fermionic partition function is impossible using PIMD since its integrand is not positive definite. We show that PIMD simulations for fermions are feasible by employing our recently developed method for bosonic PIMD and reweighting the results to obtain fermionic expectation values. The approach is tested against path integral Monte Carlo (PIMC) simulations for up to 7 electrons in a two-dimensional quantum dot for a range of interaction strengths. However, like PIMC, the method suffers from the sign problem at low temperatures. We propose a simple approach for alleviating it by simulating an auxiliary system with a larger average sign and obtaining an upper bound to the energy of the original system using the Bogoliubov inequality. This allows fermions to be studied at temperatures lower than would otherwise have been feasible using PIMD, as demonstrated in the case of a three-electron quantum dot. Our results extend the boundaries of PIMD simulations of fermions and will hopefully stimulate the development of new approaches for tackling the sign problem.

We present a method to sample reactive pathways via biased molecular dynamics simulations in trajectory space. We show that the use of enhanced sampling techniques enables unconstrained exploration of multiple reaction routes. Time correlation functions are conveniently computed via reweighted averages along a single trajectory and kinetic rates are accessed at no additional cost. These abilities are illustrated analyzing a model potential and the umbrella inversion of NH$_3$ in water. The algorithm allows a parallel implementation and promises to be a powerful tool for the study of rare events., Comment: Main Text: 5 pages, 5 figures, 7 equations. Supplemental Material: 6 pages, 11 figures, 11 equations

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.

Many processes in chemistry and physics take place on timescales that cannot be explored using standard molecular dynamics simulations. This renders the use of enhanced sampling mandatory. Here we introduce an enhanced sampling method that is based on constructing a model probability density from which a bias potential is derived. The model relies on the fact that in a physical system most of the configurations visited can be grouped into isolated metastable islands. To each island we associate a distribution that is fitted to a Gaussian mixture. The different distributions are linearly combined together with coefficients that are computed self consistently. Remarkably, from this biased dynamics, rates of transition between different metastable states can be straightforwardly computed.

Metadynamics is an enhanced sampling method of great popularity, based on the on-the-fly construction of a bias potential that is function of a selected number of collective variables. We propose here a change in perspective that shifts the focus from the bias to the probability distribution reconstruction, while keeping some of the key characteristics of metadynamics, such as the flexible on-the-fly adjustments to the free energy estimate. The result is an enhanced sampling method that presents a drastic improvement in convergence speed, especially when dealing with suboptimal and/or multidimensional sets of collective variables. The method is especially robust and easy to use, in fact it requires only few simple parameters to be set, and it has a straightforward reweighting scheme to recover the statistics of the unbiased ensemble. Furthermore it gives more control on the desired exploration of the phase space, since the deposited bias is not allowed to grow indefinitely and it does not push the simulation to uninteresting high free energy regions. We demonstrate the performance of the method in a number of representative examples.

A widely spread method of crystal preparation is to precipitate it from a supersaturated solution. In such a process, control of solution concentration is of paramount importance. Nucleation process, polymorph selection, and crystal habits depend crucially on this thermodynamic parameter. When performing simulations in the canonical ensemble as the crystalline phase is deposited the solution is depleted of solutes. This unavoidable modification of the thermodynamic conditions leads to significant artifact. Here we adopt the idea of the constant chemical potential molecular dynamics approach of Perego et al. [J. Chem. Phys. 2015, 142, 144113] to the study of nucleation. Our method allows determining the crystal nucleus size and nucleation rates at constant supersaturation. As an example we study the homogeneous nucleation of sodium chloride from its supersaturated aqueous solution.

Trapped Bosons exhibit fundamental physical phenomena and are potentially useful for quantum technologies. We present a method for simulating Bosons using path integral molecular dynamics. A main challenge for simulations is including all permutations due to exchange symmetry. We show that evaluation of the potential can be done recursively, avoiding explicit enumeration of permutations, and scales cubically with system size. The method is applied to Bosons in a 2D trap and agrees with essentially exact results. An analysis of the role of exchange with decreasing temperature is also presented.

Aqueous solutions provide a medium for many important reactions in chemical synthesis, industrial processes, environmental chemistry, and biological functions. It is an accepted fact that aqueous solvents can be direct participants to the reaction process and not act only as simple passive dielectrics. Assisting water molecules and proton wires are thus essential for the efficiency of many reactions. Here we study the decomposition of urea into ammonia and isocyanic acid by means of enhanced ab initio molecular dynamics simulations. We highlight the role of the solvent molecules and their interactions with the reactants providing a proper description of the reaction mechanism and how the water hydrogen-bond network affects the reaction dynamics. Reaction free energy and rates have been calculated taking into account this important effect.

Determining the complete set of ligands' binding/unbinding pathways is important for drug discovery and to rationally interpret mutation data. Here we have developed a metadynamics-based technique that addressed this issue and allows estimating affinities in the presence of multiple escape pathways. Our approach is shown on a Lysozyme T4 variant in complex with the benzene molecule. The calculated binding free energy is in agreement with experimental data. Remarkably, not only we were able to find all the previously identified ligand binding pathways, but also we uncovered 3 new ones. This results were obtained at a small computational cost, making this approach valuable for practical applications, such as screening of small compounds libraries.

We used molecular dynamics simulations to predict the steady state crystal shape of naphthalene grown from ethanol solution. The simulations were performed at constant supersaturation by utilizing a recently proposed algorithm [Perego et al., J. Chem. Phys., 142, 2015, 144113]. To bring the crystal growth within the timescale of a molecular dynamics simulation we applied Well-Tempered Metadynamics with a spatially constrained collective variable, which focuses the sampling on the growing layer. We estimated that the resulting steady state crystal shape corresponds to a rhombic prism, which is in line with experiments. Further, we observed that at the investigated supersaturations, the $\{00\bar{1}\}$ face grows in a two step two dimensional nucleation mechanism while the considerably faster growing faces $\{1\bar{1}0\}$ and $\{20\bar{1}\}$ grow new layers with a one step two dimensional nucleation mechanism.

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 not to 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.

From the Ising model and the Lennard-Jones fluid, to water and the iron-carbon system, phase diagrams are an indispensable tool to understand phase equilibria. In spite of the effort of the simulation community the calculation of a large portion of a phase diagram using computer simulation is still today a significant challenge. Here we propose a method to calculate phase diagrams involving liquid and solid phases by the reversible transformation of the liquid and the solid. To this end we introduce an order parameter that breaks the rotational symmetry and we leverage our recently introduced method to sample the multithermal-multibaric ensemble. In this way in a single molecular dynamics simulation we are able to compute the liquid-solid coexistence line for entire regions of the temperature and pressure phase diagram. We apply our approach to the bcc-liquid phase diagram of sodium and the fcc-bcc-liquid phase diagram of aluminum., Comment: 9 pages, 4 figures

Free energy sampling methods allow studying the full dynamics of activated processes. Unfortunately, the affordable accuracy of the potential describing the energy and forces of the system is usually rather low. Here we introduce a new method that combining metadynamics and free energy perturbation allows calculating accurate quantum chemical free energies for chemical reactions. To prove the effectiveness of this new approach we study the S N 2 reaction of CH3F + Cl- -> CH3Cl + F- in vacuo and solvated by water. Comparisons are made with harmonic transition state theory to show how this method could provide accurate equilibrium and rate constants for complex systems.