Estimating Free-Energy Surfaces and Their Convergence from Multiple, Independent Static and History-Dependent Biased Molecular-Dynamics Simulations with Mean Force Integration

Addressing the sampling problem is central to obtaining quantitative insight from molecular dynamics simulations. Adaptive biased sampling methods, such as metadynamics, tackle this issue by perturbing the Hamiltonian of a system with a history-dependent bias potential, enhancing the exploration of the ensemble of configurations and estimating the corresponding free energy surface (FES). Nevertheless, efficiently assessing and systematically improving their convergence remains an open problem. Here, building on mean force integration (MFI), we develop and test a metric for estimating the convergence of FESs obtained by combining asynchronous, independent simulations subject to diverse biasing protocols, including static biases, different variants of metadynamics, and various combinations of static and history-dependent biases. The developed metric and the ability to combine independent simulations granted by MFI enable us to devise strategies to systematically improve the quality of FES estimates. We demonstrate our approach by computing FES and assessing the convergence of a range of systems of increasing complexity, including one- and two-dimensional analytical FESs, alanine dipeptide, a Lennard-Jones supersaturated vapor undergoing liquid droplet nucleation, and the model of a colloidal system crystallizing via a two-step mechanism. The methods presented here can be generally applied to biased simulations and are implemented in pyMFI, a publicly accessible, open-source Python library

A variational approach to assess reaction coordinates for two-step crystallization

Molecule- and particle-based simulations provide the tools to test, in microscopic detail, the validity of classical nucleation theory. In this endeavor, determining nucleation mechanisms and rates for phase separation requires an appropriately defined reaction coordinate to describe the transformation of an out-of-equilibrium parent phase for which myriad options are available to the simulator. In this article, we describe the application of the variational approach to Markov processes to quantify the suitability of reaction coordinates to study crystallization from supersaturated colloid suspensions. Our analysis indicates that collective variables (CVs) that correlate with the number of particles in the condensed phase, the system potential energy, and approximate configurational entropy often feature as the most appropriate order parameters to quantitatively describe the crystallization process. We apply time-lagged independent component analysis to reduce high-dimensional reaction coordinates constructed from these CVs to build Markov State Models (MSMs), which indicate that two barriers separate a supersaturated fluid phase from crystals in the simulated environment. The MSMs provide consistent estimates for crystal nucleation rates, regardless of the dimensionality of the order parameter space adopted; however, the two-step mechanism is only consistently evident from spectral clustering of the MSMs in higher dimensions. As the method is general and easily transferable, the variational approach we adopt could provide a useful framework to study controls for crystal nucleation

Into the dynamics of a supramolecular polymer at submolecular resolution

To rationally design supramolecular polymers capable of self-healing or reconfiguring their structure in a dynamically controlled way, it is imperative to gain access into the intrinsic dynamics of the supramolecular polymer (dynamic exchange of monomers) while maintaining a high-resolution description of the monomer structure. But this is prohibitively difficult at experimental level. Here we show atomistic, coarse-grained modelling combined with advanced simulation approaches to characterize the molecular mechanisms and relative kinetics of monomer exchange in structural variants of a synthetic supramolecular polymer in different conditions. We can capture differences in supramolecular dynamics consistent with the experimental observations, revealing that monomer exchange in and out the fibres originates from the defects present in their supramolecular structure. At the same time, the submolecular resolution of this approach offers a molecular-level insight into the dynamics of these bioinspired materials, and a flexible tool to obtain structure-dynamics relationships for a variety of polymeric assemblies

Communication: Role of explicit water models in the helix folding/unfolding processes

In the last years, it has become evident that computer simulations can assume a relevant role in modelling protein dynamical motions for their ability to provide a full atomistic image of the processes under investigation. The ability of the current protein force-fields in reproducing the correct thermodynamics and kinetics systems behaviour is thus an essential ingredient to improve our understanding of many relevant biological functionalities. In this work, employing the last developments of the metadynamics framework, we compare the ability of state-of-the-art all-atom empirical functions and water models to consistently reproduce the folding and unfolding of a helix turn motif in a model peptide. This theoretical study puts in evidence that the choice of the water models can influence the thermodynamic and the kinetics of the system under investigation, and for this reason cannot be considered trivial

Reducing crystal structure overprediction of ibuprofen with large scale molecular dynamics simulations

The control of the crystal form is a central issue in the pharmaceutical industry. The identification of putative polymorphs through Crystal Structure Prediction (CSP) methods is based on lattice energy calculations, which are known to significantly over-predict the number of plausible crystal structures. A valuable tool to reduce overprediction is to employ physics-based, dynamic simulations to coalesce lattice energy minima separated by small barriers into a smaller number of more stable geometries once thermal effects are introduced. Molecular dynamics simulations and enhanced sampling methods can be employed in this context to simulate crystal structures at finite temperature and pressure. Here we demonstrate the applicability of approaches based on molecular dynamics to systematically process realistic CSP datasets containing several hundreds of crystal structures. The system investigated is ibuprofen, a conformationally flexible active pharmaceutical ingredient that crystallises both in enantiopure forms and as a racemic mixture. By introducing a hierarchical approach in the analysis of finite-temperature supercell configurations, we can post-process a dataset of 555 crystal structures, identifying 65% of the initial structures as labile, while maintaining all the experimentally known crystal structures in the final, reduced set. Moreover, the extensive nature of the initial dataset allows one to gain quantitative insight into the persistence and the propensity to transform of crystal structures containing common hydrogen-bonded intermolecular interaction motifs

Solvent Dynamics and Thermodynamics at the Crystal-Solution Interface of Ibuprofen

The choice of solvent is key in the manufacturing of solution-grown crystals due to the critical effect it can exert on their morphology. Here we set out to investigate the dynamics and thermodynamics of solvent molecules at the crystal-solution interface for the morphologically dominant crystal faces of ibuprofen. In particular, we evaluate how thermodynamically favourable the desorption of a solvent molecule is and estimate the rate of exchange of adsorbed solvent molecules with molecules from the bulk solution. This analysis is carried out for all four morphologically dominant crystal faces of ibuprofen {100}, {002}, {011} and {110}, and ten solvents, i.e. water, 1-butanol, toluene, cyclohexanone, cyclohexane, acetonitrile, trichloromethane, methanol, ethyl acetate and ethanol. Our work reveals that the structure of the solution and the exchange dynamics can be strongly dependent on both the crystal face and the solvent, i.e. the same solvent can show radically different structure when in contact with different faces, alternatively the same face can induce different structuring in different solvents. Moreover, we find particularly strong surface-solvent interactions for the {002} and {100} crystal faces in several of the solvents examined. We conclude that the role of desolvation in the growth process is solvent- and face-specific, and therefore it has the potential of impacting the crystal shape anisotropy. We provide a framework to rationalise this effect based on molecular simulations of the crystal/solution interface

Electrochemistry, ion adsorption and dynamics in the double layer: a study of NaCl(aq) on graphite

Graphite and related sp2 carbons are ubiquitous electrode materials with particular promise for use in e.g., energy storage and desalination devices, but very little is known about the properties of the carbon–electrolyte double layer at technologically relevant concentrations. Here, the (electrified) graphite–NaCl(aq) interface was examined using constant chemical potential molecular dynamics (CμMD) simulations; this approach avoids ion depletion (due to surface adsorption) and maintains a constant concentration, electroneutral bulk solution beyond the surface. Specific Na+ adsorption at the graphite basal surface causes charging of the interface in the absence of an applied potential. At moderate bulk concentrations, this leads to accumulation of counter-ions in a diffuse layer to balance the effective surface charge, consistent with established models of the electrical double layer. Beyond ∼0.6 M, however, a combination of over-screening and ion crowding in the double layer results in alternating compact layers of charge density perpendicular to the interface. The transition to this regime is marked by an increasing double layer size and anomalous negative shifts to the potential of zero charge with incremental changes to the bulk concentration. Our observations are supported by changes to the position of the differential capacitance minimum measured by electrochemical impedance spectroscopy, and are explained in terms of the screening behaviour and asymmetric ion adsorption. Furthermore, a striking level of agreement between the differential capacitance from solution evaluated in simulations and measured in experiments allows us to critically assess electrochemical capacitance measurements which have previously been considered to report simply on the density of states of the graphite material at the potential of zero charge. Our work shows that the solution side of the double layer provides the more dominant contribution to the overall measured capacitance. Finally, ion crowding at the highest concentrations (beyond ∼5 M) leads to the formation of liquid-like NaCl clusters confined to highly non-ideal regions of the double layer, where ion diffusion is up to five times slower than in the bulk. The implications of changes to the speciation of ions on reactive events in the double layer are discussed

Explicit temperature coupling in phase-field crystal models of solidification

We present a phase-field crystal (PFC) model for solidification that accounts for thermal transport and a temperature-dependent lattice parameter. Elasticity effects are characterized through the continuous elastic field computed from the microscopic density field. We showcase the model capabilities via selected numerical investigations which focus on the prototypical growth of two-dimensional crystals from the melt, resulting in faceted shapes and dendrites. This work sets the grounds for a comprehensive mesoscale model of solidification including thermal expansion

Improved time integration for phase-field crystal models of solidification

We optimize a numerical time-stabilization routine for the phase-field crystal (PFC) models of solidification. By numerical experiments, we showcase that our approach can improve the accuracy of underlying time integration schemes by a few orders of magnitude. We investigate different time integration schemes. Moreover, as a prototypical example for applications, we extend our numerical approach to a PFC model of solidification with an explicit temperature coupling.Comment: 5 pages, 3 figure

Speciation of Substituted Benzoic Acids in Solution: Evaluation of Spectroscopic and Computational Methods for the Identification of Associates and Their Role in Crystallization

Self-association of four benzoic acid derivatives 2-chloro-4-nitrobenzoic acid, 2-methyl-4-nitrobenzoic acid, 3-hydroxybenzoic acid, and 2,6-dimethoxybenzoic acid in solution was investigated using spectroscopic measurements (Fourier transform infrared and 1H and 13C NMR spectroscopy) and molecular simulation methods. Based on the formation of hydrogen bonds, solvents can be divided into two groups: apolar solvents or solvents with a low hydrogen bond acceptor propensity, in which the benzoic acid derivatives form hydrogen-bonded dimers, and solvents with hydrogen bond acceptor propensity β > 0.3, interacting with the carboxylic group of benzoic acid, thus screening its interaction in the formation of self-associates. The formation propensity and structure of self-associates stabilized by weak interactions, such as π···π stacking and CH3···π interactions, however, are determined by the substituents in the benzene ring. Despite all the studied compounds being polymorphic, in none of the cases, an unequivocal structural link between self-associates present in the solution and the crystal form was observed