Non-Equilibrium Modeling of Concentration-Driven processes with Constant Chemical Potential Molecular Dynamics Simulations.

ConspectusConcentration-driven processes in solution, i.e., phenomena that are sustained by persistent concentration gradients, such as crystallization and surface adsorption, are fundamental chemical processes. Understanding such phenomena is crucial for countless applications, from pharmaceuticals to biotechnology. Molecular dynamics (MD), both in- and out-of-equilibrium, plays an essential role in the current understanding of concentration-driven processes. Computational costs, however, impose drastic limitations on the accessible scale of simulated systems, hampering the effective study of such phenomena. In particular, due to these size limitations, closed system MD of concentration-driven processes is affected by solution depletion/enrichment that unavoidably impacts the dynamics of the chemical phenomena under study. As a notable example, in simulations of crystallization from solution, the transfer of monomers between the liquid and crystal phases results in a gradual depletion/enrichment of solution concentration, altering the driving force for phase transition. In contrast, this effect is negligible in experiments, given the macroscopic size of the solution volume. Because of these limitations, accurate MD characterization of concentration-driven phenomena has proven to be a long-standing simulation challenge. While disparate equilibrium and nonequilibrium simulation strategies have been proposed to address the study of such processes, the methodologies are in continuous development.In this context, a novel simulation technique named constant chemical potential molecular dynamics (CµMD) was recently proposed. CµMD employs properly designed, concentration-dependent external forces that regulate the flux of solute species between selected subregions of the simulation volume. This enables simulations of systems under a constant chemical drive in an efficient and straightforward way. The CµMD scheme was originally applied to the case of crystal growth from solution and then extended to the simulation of various physicochemical processes, resulting in new variants of the method. This Account illustrates the CµMD method and the key advances enabled by it in the framework of in silico chemistry. We review results obtained in crystallization studies, where CµMD allows growth rate calculations and equilibrium shape predictions, and in adsorption studies, where adsorption thermodynamics on porous or solid surfaces was correctly characterized via CµMD. Furthermore, we will discuss the application of CµMD variants to simulate permeation through porous materials, solution separation, and nucleation upon fixed concentration gradients. While presenting the numerous applications of the method, we provide an original and comprehensive assessment of concentration-driven simulations using CµMD. To this end, we also shed light on the theoretical and technical foundations of CµMD, underlining the novelty and specificity of the method with respect to existing techniques while stressing its current limitations. Overall, the application of CµMD to a diverse range of fields provides new insight into many physicochemical processes, the in silico study of which has been hitherto limited by finite-size effects. In this context, CµMD stands out as a general-purpose method that promises to be an invaluable simulation tool for studying molecular-scale concentration-driven phenomena.

Properties of aqueous electrolyte solutions at carbon electrodes: effects of concentration and surface charge on solution structure, ion clustering and thermodynamics in the electric double layer.

Surfaces are able to control physical-chemical processes in multi-component solution systems and, as such, find application in a wide range of technological devices. Understanding the structure, dynamics and thermodynamics of non-ideal solutions at surfaces, however, is particularly challenging. Here, we use Constant Chemical Potential Molecular Dynamics (CµMD) simulations to gain insight into aqueous NaCl solutions in contact with graphite surfaces at high concentrations and under the effect of applied surface charges: conditions where mean-field theories describing interfaces cannot (typically) be reliably applied. We discover an asymmetric effect of surface charge on the electric double layer structure and resulting thermodynamic properties, which can be explained by considering the affinity of the surface for cations and anions and the cooperative adsorption of ions that occurs at higher concentrations. We characterise how the sign of the surface charge affects ion densities and water structure in the double layer and how the capacitance of the interface-a function of the electric potential drop across the double layer-is largely insensitive to the bulk solution concentration. Notably, we find that negatively charged graphite surfaces induce an increase in the size and concentration of extended liquid-like ion clusters confined to the double layer. Finally, we discuss how concentration and surface charge affect the activity coefficients of ions and water at the interface, demonstrating how electric fields in this region should be explicitly considered when characterising the thermodynamics of both solute and solvent at the solid/liquid interface.

Exploring the Binding Pathway of Novel Nonpeptidomimetic Plasmepsin V Inhibitors.

Predicting the interaction modes and binding affinities of virtual compound libraries is of great interest in drug development. It reduces the cost and time of lead compound identification and selection. Here we apply path-based metadynamics simulations to characterize the binding of potential inhibitors to the Plasmodium falciparum aspartic protease plasmepsin V (plm V), a validated antimalarial drug target that has a highly mobile binding site. The potential plm V binders were identified in a high-throughput virtual screening (HTVS) campaign and were experimentally verified in a fluorescence resonance energy transfer (FRET) assay. Our simulations allowed us to estimate compound binding energies and revealed relevant states along binding/unbinding pathways in atomistic resolution. We believe that the method described allows the prioritization of compounds for synthesis and enables rational structure-based drug design for targets that undergo considerable conformational changes upon inhibitor binding.

Constant chemical potential-quantum mechanical-molecular dynamics simulations of the graphene-electrolyte double layer.

We present the coupling of two frameworks-the pseudo-open boundary simulation method known as constant potential molecular dynamics simulations (CµMD), combined with quantum mechanics/molecular dynamics (QMMD) calculations-to describe the properties of graphene electrodes in contact with electrolytes. The resulting CµQMMD model was then applied to three ionic solutions (LiCl, NaCl, and KCl in water) at bulk solution concentrations ranging from 0.5 M to 6 M in contact with a charged graphene electrode. The new approach we are describing here provides a simulation protocol to control the concentration of electrolyte solutions while including the effects of a fully polarizable electrode surface. Thanks to this coupling, we are able to accurately model both the electrode and solution side of the double layer and provide a thorough analysis of the properties of electrolytes at charged interfaces, such as the screening ability of the electrolyte and the electrostatic potential profile. We also report the calculation of the integral electrochemical double layer capacitance in the whole range of concentrations analyzed for each ionic species, while the quantum mechanical simulations provide access to the differential and integral quantum capacitance. We highlight how subtle features, such as the adsorption of potassium graphene or the tendency of the ions to form clusters contribute to the ability of graphene to store charge, and suggest implications for desalination.

Understanding Metal-Organic Framework Nucleation from a Solution with Evolving Graphs.

A mechanistic understanding of metal-organic framework (MOF) synthesis and scale-up remains underexplored due to the complex nature of the interactions of their building blocks. In this work, we investigate the collective assembly of building units at the early stages of MOF nucleation, using MIL-101(Cr) as a prototypical example. Using large-scale molecular dynamics simulations, we observe that the choice of solvent (water and N,N-dimethylformamide), the introduction of ions (Na+ and F-) and the relative populations of MIL-101(Cr) half-secondary building unit (half-SBU) isomers have a strong influence on the cluster formation process. Additionally, the shape, size, nucleation and growth rates, crystallinity, and short and long-range order largely vary depending on the synthesis conditions. We evaluate these properties as they naturally emerge when interpreting the self-assembly of MOF nuclei as the time evolution of an undirected graph. Solution-induced conformational complexity and ionic concentration have a dramatic effect on the morphology of clusters emerging during assembly. While pure solvents lead to the rapid formation of a small number of large clusters, the presence of ions in aqueous solutions results in smaller clusters and slower nucleation. This diversity is captured by the key features of the graph representation. Principle component analysis on graph properties reveals that only a small number of molecular descriptors is needed to deconvolute MOF self-assembly. Descriptors such as the average coordination number between half-SBUs and fractal dimension are of particulalr interest as they can be can be followed experimentally by techniques like by time-resolved spectroscopy. Ultimately, graph theory emerges as an approach that can be used to understand complex processes revealing molecular descriptors accessible by both simulation and experiment.

Droplet Microfluidics for the Label-Free Extraction of Complete Phase Diagrams and Kinetics of Liquid-Liquid Phase Separation in Finite Volumes.

Liquid-liquid phase separation of polymer and protein solutions is central in many areas of biology and material sciences. Here, an experimental and theoretical framework is provided to investigate the thermodynamics and kinetics of liquid-liquid phase separation in volumes comparable to cells. The strategy leverages droplet microfluidics to accurately measure the volume of the dense phase generated by liquid-liquid phase separation of solutions confined in micro-sized compartments. It is shown that the measurement of the volume fraction of the dense phase at different temperatures allows the evaluation of the binodal lines that determine the coexistence region of the two phases in the temperature-concentration phase diagram. By applying a thermodynamic model of phase separation in finite volumes, it is further shown that the platform can predict and validate kinetic barriers associated with the formation of a dense droplet in a parent dilute phase, therefore connecting thermodynamics and kinetics of liquid-liquid phase separation.

Multiple pathways in NaCl homogeneous crystal nucleation.

NaCl crystal nucleation from metastable solutions has long been considered to occur according to a single-step mechanism where the growth in the size and crystalline order of the emerging nuclei is simultaneous. Recent experimental observations suggest that significant ion-ion correlations occur in solution and that NaCl crystals can emerge from disordered intermediates which is seemingly at odds with this established view. Here, we performed biased and unbiased molecular dynamics simulations to analyse and characterise the pathways to crystalline phases from solutions far into the metastable region. We find that large liquid-like NaCl clusters emerge as the solution concentration is increased and a wide distribution of crystallisation pathways are observed with two-step nucleation pathways-where crystalline order emerges in dense liquid NaCl regions-being more dominant than one-step pathways to phase separation far into the metastable region. Analyses of cluster size populations and the ion pair association constant show that these clusters are transient, unlike the thermodynamically stable prenucleation cluster solute species that were suggested in other mineralising systems. A Markov state model was developed to analyse the mechanisms and timescales for nucleation from unbiased molecular dynamics trajectories in a reaction coordinate space characterising the dense regions in clusters and crystalline order. This allowed calculation of the committor probabilities for the system to relax to the solution or crystal states and to estimate the rate of nucleation, which shows excellent agreement with literature values. From a fundamental nucleation perspective, our work highlights the need to extend the attribute 'critical' to an ensemble of clusters which can display a broad range of structures and include sizeable disordered domains depending upon the reaction conditions. Moreover, our recent simulation studies demonstrated that carbon surfaces catalyse the formation of liquid-like NaCl networks which, combined with the observations here, suggests that alternative pathways beyond the single-step mechanism can be exploited to control the crystallisation of NaCl.

Investigating the quasi-liquid layer on ice surfaces: a comparison of order parameters.

Ice surfaces are characterized by pre-melted quasi-liquid layers (QLLs), which mediate both crystal growth processes and interactions with external agents. Understanding QLLs at the molecular level is necessary to unravel the mechanisms of ice crystal formation. Computational studies of the QLLs heavily rely on the accuracy of the methods employed for identifying the local molecular environment and arrangements, discriminating between solid-like and liquid-like water molecules. Here we compare the results obtained using different order parameters to characterize the QLLs on hexagonal ice (Ih) and cubic ice (Ic) model surfaces investigated with molecular dynamics (MD) simulations in a range of temperatures. For the classification task, in addition to the traditional Steinhardt order parameters in different flavours, we select an entropy fingerprint and a deep learning neural network approach (DeepIce), which are conceptually different methodologies. We find that all the analysis methods give qualitatively similar trends for the behaviours of the QLLs on ice surfaces with temperature, with some subtle differences in the classification sensitivity limited to the solid-liquid interface. The thickness of QLLs on the ice surface increases gradually as the temperature increases. The trends of the QLL size and of the values of the order parameters as a function of temperature for the different facets may be linked to surface growth rates which, in turn, affect crystal morphologies at lower vapour pressure. The choice of the order parameter can be therefore informed by computational convenience except in cases where a very accurate determination of the liquid-solid interface is important.

Identifying Conformational Isomers of Organic Molecules in Solution via Unsupervised Clustering.

We present a systematic approach for the identification of statistically relevant conformational macrostates of organic molecules from molecular dynamics trajectories. The approach applies to molecules characterized by an arbitrary number of torsional degrees of freedom and enables the transferability of the macrostates definition across different environments. We formulate a dissimilarity measure between molecular configurations that incorporates information on the characteristic energetic cost associated with transitions along all relevant torsional degrees of freedom. Such metric is employed to perform unsupervised clustering of molecular configurations based on the Fast Search and Find of Density Peaks algorithm. We apply this method to investigate the equilibrium conformational ensemble of Sildenafil, a conformationally complex pharmaceutical compound, in different environments including the crystal bulk, the gas phase, and three different solvents (acetonitrile, 1-butanol, and toluene). We demonstrate that while Sildenafil can adopt more than 100 metastable conformational configurations, only 12 are significantly populated across all of the environments investigated. Despite the complexity of the conformational space, we find that the most abundant conformers in solution are the closest to the conformers found in the most common Sildenafil crystal phase.

Molecular Level Understanding of the Free Energy Landscape in Early Stages of Metal-Organic Framework Nucleation.

The assembly mechanism of  Metal-Organic Frameworks (MOFs) is controlled by the choice of solvent and the presence of spectator ions. In this paper, we apply  enhanced sampling molecular dynamics methods to investigate the role of solvent and ions in the early stages of the synthesis of MIL-101(Cr). Microsecond-long well-tempered metadynamics simulations uncover a  rich  structural free energy landscape, with secondary building units (SBUs) adopting distinct crystal and noncrystal like configurations. In the presence of ions (Na+, F-), we observe a complex effect on the crystallinity of SBUs. By  modulating the interactions between terephthalate linkers and Cr atoms, ions affect the abundance of crystal-like SBUs, consequently controlling the percentage of defects. Solvent effects are assessed by comparing water with   N, N-dimethylformamide, in which SBU adducts are appreciably more stable and compact. These results shed light on how solvent and ionic strength impact the free energy of the assembly phenomena that ultimately control material synthesis.