Search

In computational physics, chemistry, and biology, the implementation of new techniques in a shared and open source software lowers barriers to entry and promotes rapid scientific progress. However, effectively training new software users presents several challenges. Common methods like direct knowledge transfer and in-person workshops are limited in reach and comprehensiveness. Furthermore, while the COVID-19 pandemic highlighted the benefits of online training, traditional online tutorials can quickly become outdated and may not cover all the software's functionalities. To address these issues, here we introduce ``PLUMED Tutorials'', a collaborative model for developing, sharing, and updating online tutorials. This initiative utilizes repository management and continuous integration to ensure compatibility with software updates. Moreover, the tutorials are interconnected to form a structured learning path and are enriched with automatic annotations to provide broader context. This paper illustrates the development, features, and advantages of PLUMED Tutorials, aiming to foster an open community for creating and sharing educational resources., Comment: 26 pages, 5 figures

Many biological processes occur on time scales longer than those accessible to molecular dynamics simulations. Identifying collective variables (CVs) and introducing an external potential to accelerate them is a popular approach to address this problem. In particular, $\texttt{PLUMED}$ is a community-developed library that implements several methods for CV-based enhanced sampling. This chapter discusses two recent developments that have gained popularity in recent years. The first is the On-the-fly Probability Enhanced Sampling (OPES) method as a biasing scheme. This provides a unified approach to enhanced sampling able to cover many different scenarios: from free energy convergence to the discovery of metastable states, from rate calculation to generalized ensemble simulation. The second development concerns the use of machine learning (ML) approaches to determine CVs by learning the relevant variables directly from simulation data. The construction of these variables is facilitated by the $\texttt{mlcolvar}$ library, which allows them to be optimized in Python and then used to enhance sampling thanks to a native interface inside $\texttt{PLUMED}$. For each of these methods, in addition to a brief introduction, we provide guidelines, practical suggestions and point to examples from the literature to facilitate their use in the study of the process of interest.

A first principles machine learning model has been developed aimed at studying the formation of calcium carbonate from aqueous solution using molecular dynamics simulations. The model, dubbed SCAN-ML, reproduces accurately the potential energy surface derived from ab initio density-functional theory within the SCAN approximation for the exchange and correlation functional. A broad range of properties have been calculated relevant to ions in solution, solid phases, and the calcite/water interface. Careful comparison with results from experiments and semi-empirical force fields shows that SCAN-ML provides an excellent description of this system, surpassing state-of-the-art force fields for many properties, while providing a benchmark for many quantities that are currently beyond the reach of direct ab initio molecular dynamics. A key feature of SCAN-ML is its ability to capture chemical reactions, including the carbonate to bicarbonate transformation. Our model thus paves the way for the study of reactive crystallization pathways in calcium carbonate, which are currently poorly understood., Comment: 13 pages, 6 figures

The electrical double layer (EDL) at aqueous solution-metal oxide interfaces critically affects many fundamental processes in electrochemistry, geology and biology, yet understanding its microscopic structure is challenging for both theory and experiments. Here we employ ab initio-based machine learning potentials including long-range electrostatics in large-scale atomistic simulations of the EDL at the TiO2-electrolyte interface. Our simulations provide a molecular-scale picture of the EDL that demonstrates the limitations of standard mean-field models. We further develop a method to accurately calculate the electrostatic potential drop at the interface. The computed capacitance originating from the adsorbed charges and the potential drop agrees with experiments, supporting the reliability of our description of the EDL. The larger interfacial capacitance of basic relative to acidic solutions originates from the higher affinity of the cations for the oxide surface and gives rise to distinct charging mechanisms on negative and positive surfaces.

Understanding the condensed-phase behavior of chiral molecules is important in biology, as well as in a range of technological applications, such as the manufacture of pharmaceuticals. Here, we use molecular dynamics simulations to study a chiral four-site molecular model that exhibits a second-order symmetry-breaking phase transition from a supercritical racemic liquid, into subcritical D-rich and L-rich liquids. We determine the infinite-size critical temperature using the fourth-order Binder cumulant, and we show that the finite-size scaling behavior of the order parameter is compatible with the 3D Ising universality class. We also study the spontaneous D-rich to L-rich transition at a slightly subcritical temperature $T\approx0.985 T_c$ and our findings indicate that the free energy barrier for this transformation increases with system size as $N^{2/3}$ where $N$ is the number of molecules, consistent with a surface-dominated phenomenon. The critical behavior observed herein suggests a mechanism for chirality selection in which a liquid of chiral molecules spontaneously forms a phase enriched in one of the two enantiomers as the temperature is lowered below the critical point. Furthermore, the increasing free energy barrier with system size indicates that fluctuations between the L-rich and D-rich phases are suppressed as the size of the system increases, trapping it in one of the two enantiomerically-enriched phases. Such a process could provide the basis for an alternative explanation for the origin of biological homochirality. We also conjecture the possibility of observing nucleation at subcritical temperatures under the action of a suitable chiral external field., Comment: 9 pages, 7 figures

The formation of ice in the atmosphere affects precipitation and cloud properties, and plays a key role in the climate of our planet. Although ice can form directly from liquid water at deeply supercooled conditions, the presence of foreign particles can aid ice formation at much warmer temperatures. Over the past decade, experiments have highlighted the remarkable efficiency of feldspar minerals as ice nuclei compared to other particles present in the atmosphere. However, the exact mechanism of ice formation on feldspar surfaces has yet to be fully understood. Here, we develop a first-principles machine-learning model for the potential energy surface aimed at studying ice nucleation at microcline feldspar surfaces. The model is able to reproduce with high fidelity the energies and forces derived from density-functional theory (DFT) based on the SCAN exchange and correlation functional. We apply the machine-learning force field to study different fully-hydroxylated terminations of the (100), (010), and (001) surfaces of microcline exposed to vacuum. Our calculations suggest that terminations that do not minimize the number of broken bonds are preferred in vacuum. We also study the structure of supercooled liquid water in contact with microcline surfaces, and find that water density correlations extend up to around 1 nm from the surfaces. Finally, we show that the force field maintains a high accuracy during the simulation of ice formation at microcline surfaces, even for large systems of around 30,000 atoms. Future work will be directed towards the calculation of nucleation free energy barriers and rates using the force field developed herein, and understanding the role of different microcline surfaces on ice nucleation., Comment: 12 pages, 6 figures

The possible existence of a liquid-liquid critical point in deeply supercooled water has been a subject of debate in part due to the challenges associated with providing definitive experimental evidence. Pioneering work by Mishima and Stanley [Nature 392, 164 (1998) and Phys. Rev. Lett. 85, 334 (2000)] sought to shed light on this problem by studying the melting curves of different ice polymorphs and their metastable continuation in the vicinity of the expected location of the liquid-liquid transition and its associated critical point. Based on the continuous or discontinuous changes in slope of the melting curves, Mishima suggested that the liquid-liquid critical point lies between the melting curves of ice III and ice V. Here, we explore this conjecture using molecular dynamics simulations with a purely-predictive machine learning model based on ab initio quantum-mechanical calculations. We study the melting curves of ices III, IV, V, VI, and XIII using this model and find that the melting lines of all the studied ice polymorphs are supercritical and do not intersect the liquid-liquid transition locus. We also find a pronounced, yet continuous, change in slope of the melting lines upon crossing of the locus of maximum compressibility of the liquid. Finally, we analyze critically the literature in light of our findings, and conclude that the scenario in which melting curves are supercritical is favored by the most recent computational and experimental evidence. Thus, although the preponderance of experimental and computational evidence is consistent with the existence of a second critical point in water, the behavior of the melting lines of ice polymorphs does not provide strong evidence in support of this viewpoint, according to our calculations., Comment: 21 pages, 4 figures, supplementary information

A longstanding question in water research is the possibility that supercooled liquid water can undergo a liquid-liquid phase transition (LLT) into high- and low-density liquids. We used several complementary molecular simulation techniques to evaluate the possibility of an LLT in an ab initio neural network model of water trained on density functional theory calculations with the SCAN exchange correlation functional. We conclusively show the existence of a first-order LLT and an associated critical point in the SCAN description of water, representing the first definitive computational evidence for an LLT in water from first principles., Comment: 14 pages, 3 figures. Supplemental material contains 11 pages, 8 figures

Molecular simulations have provided valuable insight into the microscopic mechanisms underlying homogeneous ice nucleation. While empirical models have been used extensively to study this phenomenon, simulations based on first-principles calculations have so far proven prohibitively expensive. Here, we circumvent this difficulty by using an efficient machine learning model trained on density-functional theory (DFT) energies and forces. We compute nucleation rates at atmospheric pressure, over a broad range of supercoolings, using the seeding technique and systems of up to hundreds of thousands of atoms simulated with ab initio accuracy. The key quantity provided by the seeding technique is the size of the critical cluster (i.e., a size such that the cluster has equal probabilities of growing or melting at the given supersaturation) which is used together with the equations of classical nucleation theory to compute nucleation rates. We find that nucleation rates for our model at moderate supercoolings are in good agreement with experimental measurements within the error of our calculation. We also study the impact of properties such as the thermodynamic driving force, interfacial free energy, and stacking disorder on the calculated rates., Comment: 20 pages, 5 figures

Ice Ih, the common form of ice in the biosphere, contains proton disorder. Its proton-ordered counterpart, ice XI, is thermodynamically stable below 72 K. However, even below this temperature the formation of ice XI is kinetically hindered and experimentally it is obtained by doping ice with KOH. Doping creates ionic defects that promote the migration of protons and the associated change in proton configuration. In this article, we mimic the effect of doping in molecular dynamics simulations using a bias potential that enhances the formation of ionic defects. The recombination of the ions thus formed proceeds through fast migration of the hydroxide and results in the jump of protons along a hydrogen bond loop. This provides a physical and expedite way to change the proton configuration, and to accelerate diffusion in proton configuration space. A key ingredient of this approach is a machine learning potential trained with density functional theory data and capable of modeling molecular dissociation. We exemplify the usefulness of this idea by studying the order-disorder transition using an appropriate order parameter to distinguish the proton environments in ice Ih and XI. We calculate the changes in free energy, enthalpy, and entropy associated with the transition. Our estimated entropy agrees with experiment within the error bars of our calculation., Comment: 17 pages, 10 figures

Machine learning models are rapidly becoming widely used to simulate complex physicochemical phenomena with ab initio accuracy. Here, we use one such model as well as direct density functional theory (DFT) calculations to investigate the phase equilibrium of water, hexagonal ice (Ih), and cubic ice (Ic), with an eye towards studying ice nucleation. The machine learning model is based on deep neural networks and has been trained on DFT data obtained using the SCAN exchange and correlation functional. We use this model to drive enhanced sampling simulations aimed at calculating a number of complex properties that are out of reach of DFT-driven simulations and then employ an appropriate reweighting procedure to compute the corresponding properties for the SCAN functional. This approach allows us to calculate the melting temperature of both ice polymorphs, the driving force for nucleation, the heat of fusion, the densities at the melting temperature, the relative stability of ice Ih and Ic, and other properties. We find a correct qualitative prediction of all properties of interest. In some cases, quantitative agreement with experiment is better than for state-of-the-art semiempirical potentials for water. Our results also show that SCAN correctly predicts that ice Ih is more stable than ice Ic., Comment: 20 pages, 9 figures

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.

We study the phase equilibrium between liquid water and ice Ih modeled by the TIP4P/Ice interatomic potential using enhanced sampling molecular dynamics simulations. Our approach is based on the calculation of ice Ih-liquid free energy differences from simulations that visit reversibly both phases. The reversible interconversion is achieved by introducing a static bias potential as a function of an order parameter. The order parameter was tailored to crystallize the hexagonal diamond structure of oxygen in ice Ih. We analyze the effect of the system size on the ice Ih-liquid free energy differences and we obtain a melting temperature of 270 K in the thermodynamic limit. This result is in agreement with estimates from thermodynamic integration (272 K) and coexistence simulations (270 K). Since the order parameter does not include information about the coordinates of the protons, the spontaneously formed solid configurations contain proton disorder as expected for ice Ih., Comment: 9 pages, 6 figures

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.

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.

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

We present a method for performing multithermal-multibaric molecular dynamics simulations that sample entire regions of the temperature-pressure (TP) phase diagram. The method uses a variational principle [Valsson and Parrinello, Phys. Rev. Lett. 113, 090601 (2014)] in order to construct a bias that leads to a uniform sampling in energy and volume. The intervals of temperature and pressure are taken as inputs and the relevant energy and volume regions are determined on the fly. In this way the method guarantees adequate statistics for the chosen TP region. We show that our multithermal-multibaric simulations can be used to calculate all static physical quantities for all temperatures and pressures in the targeted region of the TP plane. We illustrate our approach by studying the density anomaly of TIP4P/Ice water., Comment: 6 pages, 6 figures

In the ongoing intense quest to increase the photoconversion efficiencies of lead halide perovskites, it has become evident that optimizing the morphology of the material is essential to achieve high peformance. Despite the fact that nucleation plays a key role in controlling the crystal morphology, very little is known about the nucleation and crystal growth processes. Here, we perform metadynamics simulations of nucleation of methylammonium lead triiodide (MAPI) in order to unravel the atomistic details of perovskite crystallization from a $\gamma$-butyrolactone solution. The metadynamics trajectories show that the nucleation process takes place in several stages. Initially, dense amorphous clusters mainly consisting of lead and iodide appear from the homogeneous solution. These clusters evolve into lead iodide (PbI$_{2}$) like structures. Subsequently, methylammonium (MA$^{+}$) ions diffuse into this PbI$_{2}$-like aggregates triggering the transformation into a perovskite crystal through a solid-solid transformation. Demonstrating the crucial role of the monovalent cations in crystallization, our simulations provide key insights into the evolution of the perovskite microstructure which is essential to make high-quality perovskite based solar cells and optoelectronics.

An important characteristic that determines the behavior of a solute in water is whether it is hydrophobic or hydrophilic. The traditional classification is based on chemical experience and heuristics. However, this does not reveal how the local environment modulates this important property. We present a local fingerprint for hydrophobicity and hydrophilicity inspired by the two body contribution to the entropy. This fingerprint is an inexpensive, quantitative and physically meaningful way of studying hydrophilicity and hydrophobicity that only requires as input the water-solute radial distribution functions. We apply our fingerprint to octanol, benzene and the 20 proteinogenic amino acids. Our measure of hydrophilicity is coherent with chemical experience and, moreover, it also shows how the character of an atom can change as its environment is changed. Lastly, we use the fingerprint as collective variable in a funnel metadynamics simulation of a host-guest system. The fingerprint serves as a desolvation collective variable that enhances transitions between the bound and unbound states.

We study by computer simulation the nucleation of a supersaturated Lennard-Jones vapor into the liquid phase. The large free energy barriers to transition make the time scale of this process impossible to study by ordinary molecular dynamics simulations.Therefore we use a recently developed enhanced sampling method [Valsson and Parrinello, Phys. Rev. Lett. 113, 090601 (2014)] based on the variational determination of a bias potential. We differ from previous applications of this method in that the bias is constructed on the basis of the physical model provided by the classical theory of nucleation. We examine the technical problems associated with this approach. Our results are very satisfactory and will pave the way for calculating the nucleation rates in many systems., Comment: 3 figures

We introduce a computational method to discover polymorphs in molecular crystals at finite temperature. The method is based on reproducing the crystallization process starting from the liquid and letting the system discover the relevant polymorphs. This idea, however, conflicts with the fact that crystallization has a time scale much longer than that of molecular simulations. In order to bring the process within affordable simulation time, we enhance the fluctuations of a collective variable by constructing a bias potential with well tempered metadynamics. We use as collective variable an entropy surrogate based on an extended pair correlation function that includes the correlation between the orientation of pairs of molecules. We also propose a similarity metric between configurations based on the extended pair correlation function and a generalized Kullback-Leibler divergence. In this way, we automatically classify the configurations as belonging to a given polymorph using our metric and a hierarchical clustering algorithm. We find all relevant polymorphs for both substances and we predict new polymorphs. One of them is stabilized at finite temperature by entropic effects., Comment: 7 pages, 4 figures

We introduce a new fingerprint that allows distinguishing between liquid-like and solid-like atomic environments. This fingerprint is based on an approximate expression for the entropy projected on individual atoms. When combined with a local enthalpy, this fingerprint acquires an even finer resolution and it is capable of discriminating between different crystal structures., Comment: 6 pages, 4 figures

Crystallization is a process of great practical relevance in which rare but crucial fluctuations lead to the formation of a solid phase starting from the liquid. Like in all first order first transitions there is an interplay between enthalpy and entropy. Based on this idea, to drive crystallization in molecular simulations, we introduce two collective variables, one enthalpic and the other entropic. Defined in this way, these collective variables do not prejudge the structure the system is going to crystallize into. We show the usefulness of this approach by studying the case of sodium and aluminum that crystallize in the bcc and fcc crystalline structure, respectively. Using these two generic collective variables, we perform variationally enhanced sampling and well tempered metadynamics simulations, and find that the systems transform spontaneously and reversibly between the liquid and the solid phases., Comment: 4 pages, 2 figures