Consensus on the solubility of NaCl in water from computer simulations using the chemical potential route

Structural and dynamical properties of aqueous NaCl brines confined in kaolinite nanopores

Quantifying thermodynamics, structural, and dynamical properties of brine confined in clay pores is critical for a variety of geo-energy applications, including underground hydrogen storage (UHS) and carbon capture and sequestration (CCS). Atomistic molecular dynamics simulations are applied here to study aqueous NaCl brines within 10-Å kaolinite slit pores. NaCl concentrations are chosen at 5, 10, 12.5, and 15 wt. %, all below the solubility limit and high enough to provide statistically relevant information. The distribution of the ions within the nanopores is found not to be homogeneous. Explicitly, Na+ cations, preferentially attracted to the siloxane surface, accumulate in regions with low water density, whereas Cl- anions, attracted to the gibbsite surface of kaolinite, are found within the hydration layers. Confinement affects the properties of ions, with ion pairing being more pronounced within the pore than in bulk aqueous solutions at similar temperatures, pressures, and compositions. Conversely, the ions affect the properties of confined water. For example, the lifetime of water-water hydrogen bonds in confinement is shortened within the hydration shells; increasing salinity from 5 to 12.5 wt. % reduces the likelihood of water density fluctuations near the kaolinite surfaces, although when the NaCl concentration rises from 12.5 to 15 wt. %, Cl- anions enhance the likelihood of density fluctuations for the hydration layer near the gibbsite surface. The simulated molecular trajectories are studied further to extract diffusion coefficients. While confinement in the kaolinite nanopore reduces the mobility of all species, non-monotonic trends are observed as a function of salt concentration. The trends seem associated with the likelihood of ion pairing. Furthermore, the diffusion coefficients for the cations are predicted to be higher than those for the anions, which is contrary to what is typically observed in bulk brines. Because density fluctuations are correlated with properties such as the solubility of gases in confined water, our observations may have important implications for geo-energy applications such as UHS and CCS.

Extreme pH Tolerance in Peptide Coacervates Mediated by Multivalent Hydrogen Bonds for Enzyme-Triggered Oral Drug Delivery

Biopolymer-based complex coacervates hold promising prospects in the field of biomedicine. However, their low stability in environments with extreme pH and high salt concentrations, largely due to weakly charged biomacromolecules and insufficient understanding of their assembly processes, has hindered their practical applications in oral drug delivery. Here, we have developed Dopa-containing peptide-based complex coacervates that are stable across a wide range of pH (1-11) and salt concentrations. Large-scale all-atom molecular dynamics simulations reveal that multivalent hydrogen bonds control the assembly pathway of the coacervates and boost their stability. Systematic point mutations reveal that various multivalent molecular interactions can synergistically tune the properties of complex coacervates. Such peptide coacervates show high drug encapsulation efficacy and trypsin-triggered release, presenting great potential for oral drug delivery applications. Our multivalent hydrogen bond-mediated peptide coacervates provide new design principles of engineering functional coacervates for diverse applications.

Chemically Informed Coarse-Graining of Electrostatic Forces in Charge-Rich Biomolecular Condensates

Biomolecular condensates composed of highly charged biomolecules, such as DNA, RNA, chromatin, and nucleic-acid binding proteins, are ubiquitous in the cell nucleus. The biophysical properties of these charge-rich condensates are largely regulated by electrostatic interactions. Residue-resolution coarse-grained models that describe solvent and ions implicitly are widely used to gain mechanistic insights into the biophysical properties of condensates, offering transferability, computational efficiency, and accurate predictions for multiple systems. However, their predictive accuracy diminishes for charge-rich condensates due to the implicit treatment of solvent and ions. Here, we present Mpipi-Recharged, a residue-resolution coarse-grained model that improves the description of charge effects in biomolecular condensates containing disordered proteins, multidomain proteins, and/or disordered single-stranded RNAs. Mpipi-Recharged introduces a pair-specific asymmetric Yukawa electrostatic potential, informed by atomistic simulations. We show that this asymmetric coarse-graining of electrostatic forces captures intricate effects, such as charge blockiness, stoichiometry variations in complex coacervates, and modulation of salt concentration, without requiring explicit solvation. Mpipi-Recharged provides excellent agreement with experiments in predicting the phase behavior of highly charged condensates. Overall, Mpipi-Recharged improves the computational tools available to investigate the physicochemical mechanisms regulating biomolecular condensates, enhancing the scope of computer simulations in this field.

Ion-Ion Structural Correlation and Dynamics of Water in Aqueous NaCl Solutions with a Wide Range of Concentrations

The behavior of water in concentrated ionic solutions, including supersaturated conditions, is crucial for numerous material and energy conversion processes and fundamental research. All electrolytes whether they "structure-make" or "structure-break" the water structure lead to slower water motion. This study investigates the structure and dynamics of aqueous NaCl solutions across a wide range of concentrations. On the structural side, the primary focus is on ion-ion correlations. In terms of dynamics, we demonstrate that the slowing down of water dynamics continues even beyond the saturated state. We identify three distinct types of dynamics at large concentrations: ballistic, trapped, and diffusive. The van Hove correlation function exhibits no signs of relaxation within a time interval where particle motion is effectively halted. The system displays dynamical heterogeneities, confirmed by evaluating non-Gaussian parameters for the self-part of the van Hove function and identifying the mobile particles. These particles form clusters, with the largest sizes occurring when the non-Gaussian parameters are at their maximum. Additionally, we discuss the relaxation times associated with these systems using the incoherent intermediate scattering function and establish a connection with the mode-coupling theory.

Solution Thermodynamics of l-Glutamic Acid Polymorphs from Finite-Sized Molecular Dynamics Simulations

Efficiently obtaining atomic-scale thermodynamic parameters characterizing crystallization from solution is key to developing the modeling strategies needed in the quest for digital design strategies for industrial crystallization processes. Based on the thermodynamics of crystal nucleation in confined solutions, we develop a simulation framework to efficiently estimate the solubility and surface tension of organic crystals in solution from a few unbiased molecular dynamics simulations at a reference temperature. We then show that such a result can be extended with minimal computational overhead to capture the solubility curve. This enables an efficient and self-consistent estimate of the solubility and limit of solution stability associated with crystal nucleation in molecular systems from equilibrium molecular dynamics without the need for sophisticated free energy calculations. We apply our analysis to investigate the relative thermodynamic stability and aqueous solubility of the α and β polymorphs of l-glutamic acid. Our analysis enables an efficient appraisal of emergent ensemble properties associated with the thermodynamics of nucleation from solutions against experimental data, demonstrating that while the absolute solubility is still far from being quantitatively captured by an off-the-shelf point charge transferable force field, the relative polymorphic stability and solubility obtained from finite temperature simulation are consistent with the experimentally available information on glutamic acid. We foresee the ability to efficiently obtain solubility information from a limited number of computational experiments as a critical component of high-throughput polymorph screenings.

Crumbling crystals: on the dissolution mechanism of NaCl in water

Dissolution of ionic salts in water is ubiquitous, particularly for NaCl. However, an atomistic scale understanding of the process remains elusive. Simulations lend themselves conveniently to studying dissolution since they provide the spatio-temporal resolution that can be difficult to obtain experimentally. Nevertheless, the complexity of various inter- and intra-molecular interactions require careful treatment and long time scale simulations, both of which are typically hindered by computational expense. Here, we use advances in machine learning potential methodology to resolve at an ab initio level of theory the dissolution mechanism of NaCl in water. The picture that emerges is that of a steady ion-wise unwrapping of the crystal preceding its rapid disintegration, reminiscent of crumbling. The onset of crumbling can be explained by a strong increase in the ratio of the surface area to volume of the crystal. Overall, dissolution comprises a series of highly dynamical microscopic sub-processes, resulting in an inherently stochastic mechanism. These atomistic level insights contribute to the general understanding of dissolution mechanisms in other crystals, and the methodology is primed for more complex systems of recent interest such as water/salt interfaces under flow and salt crystals under confinement.

Water and Carbon Dioxide Capillary Bridges in Nanoscale Slit Pores: Effects of Temperature, Pressure, and Salt Concentration on the Water Contact Angle

We perform molecular dynamics (MD) simulations of a nanoscale water capillary bridge (WCB) surrounded by carbon dioxide over a wide range of temperatures and pressures (T = 280-400 K and carbon dioxide pressures P CO 2 ≈ 0-80 MPa). The water-carbon dioxide system is confined by two parallel silica-based surfaces (hydroxylated β-cristobalite) separated by h = 5 nm. The aim of this work is to study the WCB contact angle (θc) as a function of T and P CO 2 . Our simulations indicate that θc varies weakly with temperature and pressure: Δθc ≈ 10-20° for P CO 2 increasing from ≈0 to 80 MPa (T = 320 K); Δθc ≈ -10° for T increasing from 320 to 360 K (with a fixed amount of carbon dioxide). Interestingly, at all conditions studied, a thin film of water (1-2 water layers-thick) forms under the carbon dioxide volume. Our MD simulations suggest that this is due to the enhanced ability of water, relative to carbon dioxide, to form hydrogen-bonds with the walls. We also study the effects of adding salt (NaCl) to the WCB and corresponding θc. It is found that at the salt concentrations studied (mole fractions xNa = xCl = 3.50, 9.81%), the NaCl forms a large crystallite within the WCB with the ions avoiding the water-carbon dioxide interface and the walls surface. This results in θc being insensitive to the presence of NaCl.

Electrokinetic properties of NaCl solution via molecular dynamics simulations with scaled-charge electrolytes

Scaling ionic charges has become an alternative to polarizable force fields for representing indirect charge transfer effects in molecular simulations. In our work, we apply molecular dynamics simulations to investigate the properties of NaCl aqueous solutions in homogeneous and confined media. We compare classical integer- and scaled-charge force fields for the ions. In the bulk, we validate the force fields by computing equilibrium and transport properties and comparing them with experimental data. Integer-charge ions overestimate dielectric saturation and ionic association. Both force fields present an excess in ion-ion correlation, which leads to a deviation in the ionic conductivity at higher ionic strengths. Negatively charged quartz is used to simulate the confinement effect. Electrostatic interactions dominate counter-ion adsorption. Full-charge ions have stronger and more defined adsorption planes. We obtain the electroosmotic mobility of the solution by combining the shear plane location from non-equilibrium simulations with the ionic distribution from equilibrium simulations. From the Helmholtz-Smoluchowski equation, the zeta potential and the streaming potential coupling coefficient are computed. From an atomic-scale perspective, our molecular dynamics simulations corroborate the hypothesis of maximum packing of the Stern layer, which results in a stable and non-zero zeta potential at high salinity. The scaled-charge model representation of both properties is in excellent qualitative and quantitative agreement with experimental data. With our work, we demonstrate how useful and precise simple scaled-charge models for electrolytes can be to represent complex systems, such as the electrical double layer.

To Pair or not to Pair? Machine-Learned Explicitly-Correlated Electronic Structure for NaCl in Water

The extent of ion pairing in solution is an important phenomenon to rationalize transport and thermodynamic properties of electrolytes. A fundamental measure of this pairing is the potential of mean force (PMF) between solvated ions. The relative stabilities of the paired and solvent shared states in the PMF and the barrier between them are highly sensitive to the underlying potential energy surface. However, direct application of accurate electronic structure methods is challenging, since long simulations are required. We develop wave function based machine learning potentials with the random phase approximation (RPA) and second order Møller-Plesset (MP2) perturbation theory for the prototypical system of Na and Cl ions in water. We show both methods in agreement, predicting the paired and solvent shared states to have similar energies (within 0.2 kcal/mol). We also provide the same benchmarks for different DFT functionals as well as insight into the PMF based on simple analyses of the interactions in the system.

Solubility evolution of elemental sulfur in natural gas with a varying H2S content

CONTEXT

The solubility variations of elemental sulfur are of great importance in the prevention of sulfur deposition during the development of high-sulfur gas formations. It has been observed that the solubility varies with H2S content, which is the main solvent of elemental sulfur in natural gas. Moreover, the addition of small amounts of CH4 and/or CO2 in H2S leads to a dramatic solubility reduction of which the mechanism remains unclear. Using a modified direct coexistence method, molecular dynamics simulations are conducted to uncover the molecular mechanism of the solubility reduction. The observed solubility variations with H2S content are reproduced, and the solubility reduction is interpreted by the antisolvent effect of CH4 and CO2. While the H2S content varies in a wide range in the known high-sulfur gas formation, our simulations provide useful information for controlling the sulfur deposition in gas development.

METHODS

Molecular dynamics simulations are carried out using the LAMMPS package. The initial models are constructed with the Packmol software. The Ballone and Jones' potential is used for S8, the Galliero's potential for H2S and CO2, and the transferable potentials for phase equilibria-united atom (TraPPE-UA) force field for CH4. The time step is set to 1 fs, and the molecular trajectories of additional 2 ns after equilibrium are collected for analysis.

Modeling oceanic sedimentary methane hydrate growth through molecular dynamics simulation

The crystallization process of methane hydrates in a confined geometry resembling seabed porous silica sedimentary conditions has been studied using molecular dynamics simulations. With this objective in mind, a fully atomistic quartz silica slit pore has been designed, and the temperature stability of a methane hydrate crystalline seed in the presence of water and guest molecule methane has been analyzed. NaCl ion pairs have been added in different concentrations, simulating salinity conditions up to values higher than average oceanic conditions. The structure obtained when the hydrate crystallizes inside the pore is discussed, paying special attention to the presence of ionic doping inside the hydrate and the subsequent induced structural distortion. The shift in the hydrate stability conditions due to the increasing water salinity is discussed and compared with the case of unconfined hydrate, concluding that the influence of the confinement geometry and pore hydrophilicity produces a larger deviation in the confined hydrate phase equilibria.

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.

A streamlined molecular-dynamics workflow for computing solubilities of molecular and ionic crystals

Computing the solubility of crystals in a solvent using atomistic simulations is notoriously challenging due to the complexities and convergence issues associated with free-energy methods, as well as the slow equilibration in direct-coexistence simulations. This paper introduces a molecular-dynamics workflow that simplifies and robustly computes the solubility of molecular or ionic crystals. This method is considerably more straightforward than the state-of-the-art, as we have streamlined and optimised each step of the process. Specifically, we calculate the chemical potential of the crystal using the gas-phase molecule as a reference state, and employ the S0 method to determine the concentration dependence of the chemical potential of the solute. We use this workflow to predict the solubilities of sodium chloride in water, urea polymorphs in water, and paracetamol polymorphs in both water and ethanol. Our findings indicate that the predicted solubility is sensitive to the chosen potential energy surface. Furthermore, we note that the harmonic approximation often fails for both molecular crystals and gas molecules at or above room temperature, and that the assumption of an ideal solution becomes less valid for highly soluble substances.

Scaling Solute-Solvent Distances to Improve Solubility and Ion Paring Predictions in Rigid Ion Models

Force fields based on the rigid ion model (RIM) have been developed to accurately predict the various physical and chemical properties of salts and water. However, the combined use of these models often fails to accurately predict the solubility of salts in water. To address this issue, several approaches, such as charge scaling or reparameterization, have been proposed. Nevertheless, these methods require laborious reparameterization of nonbonded force field parameters. In this article, we propose a scaling solute-solvent distance (SSSD) method to improve force fields in predicting salt solubility without changing the solute-solute and solvent-solvent interactions in the original force fields. This method can also tune the ion pairing of salt in water. One main advantage of the SSSD method is that reparameterization of the crystal and water models is not needed. We use two RIMs for the NaCl-water system (JC-SPC/E and SD-SPC/E) and the CHARMM force field for the KCl-water system to demonstrate the improved accuracy in predicting solubility by the SSSD method. Furthermore, we use the RDG-SPC/Fw force field to show that the SSSD method can also be used to tune the ion pairing of CaCO3 in water. Limitations of this method are also discussed.

Self-Consistent Implementation of a Solvation Free Energy Framework to Predict the Salt Solubilities of Six Alkali Halides

To assess the salt solubilities of six alkali halides in aqueous systems, we proposed a thermodynamic cycle and an efficient molecular modeling methodology. The Gibbs free energy changes for vaporization, dissociation, and dissolution were calculated using the experimental data of ionic thermodynamic properties obtained from the NBS tables. Additionally, the Marcus' and Tissandier's solvation free energy data for Li+, Na+, K+, Cl-, and Br- ions were compared with the conventional solvation free energies by substituting into our self-consistent thermodynamic cycle. Furthermore, Tissandier's absolute solvation free energy data were used as the training set to refit the Lennard-Jones parameters of OPLS-AA force field for ions. To predict salt solubilities, an assumption of a pseudo-solvent was proposed to characterize the coupling work of a solute with its environment from infinitely diluted to saturated solutions, indicating that the Gibbs energy change of solvation process is a function of ionic strength. Following the self-consistency of the cycle, the newly derived formulas were used to determine the salt solubilities by interpolating the intersection of Gibbs free energy of dissolution and the zero free energy line. The refined ion parameters can also predict the structure and thermodynamic properties of aqueous electrolyte solutions, such as densities, pair correlation functions, hydration numbers, mean activity coefficients, vapor pressures, and the radial dependences of the net charge at 298.15 K and 1 bar. Our method can be used to characterize the solid-liquid equilibria of ions or charged particles in aqueous systems. Furthermore, for highly concentrated strong electrolyte systems, it is essential to introduce accurate water models and polarizable force fields.

Decoupling of Interactions between Model-Charged Peptides Reveals Key Factors Responsible for Liquid-Liquid Phase Separation

Liquid-liquid phase separation (LLPS) by disordered proteins has been shown to govern biological processes and cause numerous diseases. Therefore, a deeper understanding of the interactions and their variation with external factors is key to modulating the LLPS behavior of different systems and protecting proteins from pathological aggregation. In this context, we have looked at interactions between similarly charged peptides to understand the molecular features that may drive or prevent condensate formation under various conditions. We have studied dimer formation for model peptides where charged and noncharged amino acids have been placed alternatively. Using arginine and glutamic acid as the charged residues and varying the other residues with glycine, alanine, and proline to alter hydrophobicity, we have obtained the free-energy surface (FES) for the dimer formation for these systems under high salt concentration at two different temperatures using all-atom molecular dynamics simulations and the well-tempered metadynamics method. Our results indicate that a combination of effects such as hydrophobicity, arginine-arginine interactions, or water release from the solvation shell makes dimerization free energy more favorable for the positively charged peptides with lower flexibility. For the negatively charged peptides, the crucial role of water has been found in governing the FES. Systems having charged residues and phenylalanine in the peptide sequence also have been studied at high salt concentrations using unbiased simulations. In this case, only the positively charged peptides were found to aggregate through temperature-dependent hydrophobic and cation-π interactions. Overall, our study indicates that the negatively charged peptides are more likely to remain in the dilute phase under various conditions compared to the positively charged systems. The findings from our study would be helpful in designing and controlling systems to obtain LLPS behavior for therapeutic usage.

A Deep Potential model for liquid-vapor equilibrium and cavitation rates of water

Computational studies of liquid water and its phase transition into vapor have traditionally been performed using classical water models. Here, we utilize the Deep Potential methodology-a machine learning approach-to study this ubiquitous phase transition, starting from the phase diagram in the liquid-vapor coexistence regime. The machine learning model is trained on ab initio energies and forces based on the SCAN density functional, which has been previously shown to reproduce solid phases and other properties of water. Here, we compute the surface tension, saturation pressure, and enthalpy of vaporization for a range of temperatures spanning from 300 to 600 K and evaluate the Deep Potential model performance against experimental results and the semiempirical TIP4P/2005 classical model. Moreover, by employing the seeding technique, we evaluate the free energy barrier and nucleation rate at negative pressures for the isotherm of 296.4 K. We find that the nucleation rates obtained from the Deep Potential model deviate from those computed for the TIP4P/2005 water model due to an underestimation in the surface tension from the Deep Potential model. From analysis of the seeding simulations, we also evaluate the Tolman length for the Deep Potential water model, which is (0.091 ± 0.008) nm at 296.4 K. Finally, we identify that water molecules display a preferential orientation in the liquid-vapor interface, in which H atoms tend to point toward the vapor phase to maximize the enthalpic gain of interfacial molecules. We find that this behavior is more pronounced for planar interfaces than for the curved interfaces in bubbles. This work represents the first application of Deep Potential models to the study of liquid-vapor coexistence and water cavitation.

Direct Calculation of the Interfacial Free Energy between NaCl Crystal and Its Aqueous Solution at the Solubility Limit

Salty water is the most abundant electrolyte aqueous mixture on Earth, however, very little is known about the NaCl-saturated solution interfacial free energy (γ_{s}). Here, we provide the first direct estimation of γ_{s} for several NaCl crystallographic planes by means of the mold integration technique, a highly efficient computational method to evaluate interfacial free energies with anisotropic crystal resolution. Making use of the JC-SPC/E model, one of the most benchmarked force fields for NaCl water solutions, we measure γ_{s} of four different crystal planes, (100), (110), (111), and (112[over ¯]) with the saturated solution at normal conditions. We find high anisotropy between the different crystal orientations with values ranging from 100 to 150  mJ m^{-2}, and the average value of the distinct planes being γ[over ¯]_{s}=137(20)  mJ m^{-2}. This value for the coexistence interfacial free energy is in reasonable agreement with previous extrapolations from nucleation studies. Our Letter represents a milestone in the computational calculation of interfacial free energies between ionic crystals and aqueous solutions.

The Contribution of the Ion-Ion and Ion-Solvent Interactions in a Molecular Thermodynamic Treatment of Electrolyte Solutions

Developing molecular equations of state to treat electrolyte solutions is challenging due to the long-range nature of the Coulombic interactions. Seminal approaches commonly used are the mean spherical approximation (MSA) and the Debye-Hückel (DH) theory to account for ion-ion interactions and, often, the Born theory of solvation for ion-solvent interactions. We investigate the accuracy of the MSA and DH approaches using each to calculate the contribution of the ion-ion interactions to the chemical potential of NaCl in water, comparing these with newly computer-generated simulation data; the ion-ion contribution is isolated by selecting an appropriate primitive model with a Lennard-Jones force field to describe the solvent. A study of mixtures with different concentrations and ionic strengths reveals that the calculations from both MSA and DH theories are of similar accuracy, with the MSA approach resulting in marginally better agreement with the simulation data. We also demonstrate that the Born theory provides a good qualitative description of the contribution of the ion-solvent interactions; we employ an explicitly polar water model in these simulations. Quantitative agreement up to moderate salt concentrations and across the relevant range of temperature is achieved by adjusting the Born radius using simulation data of the free energy of solvation. We compute the radial and orientational distribution functions of the systems, thereby providing further insight on the differences observed between the theory and simulation. We thus provide rigorous benchmarks for use of the MSA, DH, and Born theories as perturbation approaches, which will be of value for improving existing models of electrolyte solutions, especially in the context of equations of state.

Hydrophobicity of arginine leads to reentrant liquid-liquid phase separation behaviors of arginine-rich proteins

Intrinsically disordered proteins rich in cationic amino acid groups can undergo Liquid-Liquid Phase Separation (LLPS) in the presence of charge-balancing anionic counterparts. Arginine and Lysine are the two most prevalent cationic amino acids in proteins that undergo LLPS, with arginine-rich proteins observed to undergo LLPS more readily than lysine-rich proteins, a feature commonly attributed to arginine's ability to form stronger cation-π interactions with aromatic groups. Here, we show that arginine's ability to promote LLPS is independent of the presence of aromatic partners, and that arginine-rich peptides, but not lysine-rich peptides, display re-entrant phase behavior at high salt concentrations. We further demonstrate that the hydrophobicity of arginine is the determining factor giving rise to the reentrant phase behavior and tunable viscoelastic properties of the dense LLPS phase. Controlling arginine-induced reentrant LLPS behavior using temperature and salt concentration opens avenues for the bioengineering of stress-triggered biological phenomena and drug delivery systems.

Obtaining Consistent Free Energies for Ion Binding at Surfaces from Solution: Pathways versus Alchemy for Determining Kink Site Stability

Ion incorporation or removal from a solid at the interface with solution is a fundamental part of crystal growth. Despite this, there have been few quantitative determinations of the thermodynamics for such processes from atomistic molecular dynamics due to the associated technical challenges. In this study, we compute the free energies for ion removal from kink sites at the interface between NaCl and water as an illustrative example. To examine the influence of the free energy technique used, we compare methods that follow an explicit pathway for dissolution with those that focus on the thermodynamics of the initial and final states using metadynamics and free energy perturbation, respectively. While the initial results of the two approaches are found to be completely different, it is demonstrated that the thermodynamics can be reconciled with appropriate corrections for the standard states, thus illustrating the need for caution in interpreting raw free energy curves for ion binding as widely found in the literature. In addition, a new efficient approach is introduced to correct for the system size dependence of kink site energies both due to the periodic interaction of charges in an inhomogeneous dielectric system and due to the dipolar interactions between pairs of kinks along a row. Ultimately, it is shown that with suitable care, both classes of free energy techniques are capable of producing kink site stabilities that are consistent with the solubility of the underlying bulk solid. However, the precise values for individual kink sites exhibit a small systematic offset, which can be ascribed to the contribution of the interfacial potential to the pathway-based results. For the case of NaCl, the free energies of the kink sites relative to a 1 M aqueous solution for Na⁺ and Cl^- are found to be surprisingly different and of opposite sign, despite the ions having very similar hydration free energies.

Computing chemical potentials of solutions from structure factors

The chemical potential of a component in a solution is defined as the free energy change as the amount of the component changes. Computing this fundamental thermodynamic property from atomistic simulations is notoriously difficult, because of the convergence issues in free energy methods and finite size effects. This paper presents the so-called S0 method, which can be used to obtain chemical potentials from static structure factors computed from equilibrium molecular dynamics simulations under the isothermal-isobaric ensemble. This new method is demonstrated on the systems of binary Lennard-Jones particles, urea--water mixtures, a NaCl aqueous solution, and a high-pressure carbon-hydrogen mixture.

Activity Coefficients and Solubilities of NaCl in Water-Methanol Solutions from Molecular Dynamics Simulations

We obtain activity coefficients and solubilities of NaCl in water-methanol solutions at 298.15 K and 1 bar from molecular dynamics (MD) simulations with the Joung-Cheatham, SPC/E, and OPLS-AA force fields for NaCl, water, and methanol, respectively. The Lorentz-Berthelot combining rules were adopted for the unlike-pair interactions of Na⁺, Cl^-, and the oxygen site in SPC/E water, and geometric combining rules were utilized for the remainder of the cross interactions. We found that the selection of appropriate combining rules is important in obtaining physically realistic solubilities. The solvent compositions studied range from pure water to pure methanol. Several salt concentrations were investigated at each solvent composition, from the lowest concentrations permitted by the system size used up to the experimental solubilities. We first calculated individual ion activity coefficients (IIACs) for Na⁺ and Cl^- from the free energy change due to the gradual insertion of a single cation or anion into the solution, accompanied by a neutralizing background. We obtained the salt solubilities by comparing the chemical potentials in solution with solid NaCl chemical potentials calculated previously using the Einstein crystal method. Mean ionic activity coefficients obtained from the IIACs are in reasonable agreement with experimental data, with deviations increasing for solutions of higher methanol content. Predictions for the salt solubility are in surprisingly good agreement with experimental data, despite well-known challenges in the simultaneous calculation of activity coefficients and solubilities with classical MD simulations. The present study demonstrates that good predictions for these two important phase equilibrium properties can be obtained for mixed-solvent electrolyte solutions using existing nonpolarizable models and further suggests that the previously proposed single ion insertion technique can be extended to complex mixed-solvent solutions as well.

An investigation of the kinetics and thermodynamics of NaCl nucleation through composite clusters

Having a good understanding of nucleation is critical for the control of many important processes, such as polymorph selection during crystallization. However, a complete picture of the molecular-level mechanisms of nucleation remains elusive. In this work, we take an in-depth look at the NaCl homogeneous nucleation mechanism through thermodynamics. Distinguished from the classical nucleation theory, we calculate the free energy of nucleation as a function of two nucleus size coordinates: crystalline and amorphous cluster sizes. The free energy surface reveals a thermodynamic preference for a nonclassical mechanism of nucleation through a composite cluster, where the crystalline nucleus is surrounded by an amorphous layer. The thickness of the amorphous layer increases with an increase in supersaturation. The computed free energy landscape agrees well with the composite cluster-free energy model, through which phase specific thermodynamic properties are evaluated. As the supersaturation increases, there is a change in stability of the amorphous phase relative to the solution phase, resulting in a change from one-step to two-step mechanism, seen clearly from the free energy profile along the minimum free energy path crossing the transition curve. By obtaining phase-specific diffusion coefficients, we construct the full mesoscopic model and present a clear roadmap for NaCl nucleation.

Homogeneous nucleation of NaCl in supersaturated solutions

The seeding method is an approximate approach to investigate nucleation that combines molecular dynamics simulations with classical nucleation theory. Recently, this technique has been successfully implemented in a broad range of nucleation studies. However, its accuracy is subject to the arbitrary choice of the order parameter threshold used to distinguish liquid-like from solid-like molecules. We revisit here the crystallization of NaCl from a supersaturated brine solution and show that consistency between seeding and rigorous methods, like Forward Flux Sampling (from previous work) or spontaneous crystallization (from this work), is achieved by following a mislabelling criterion to select such threshold (i.e. equaling the fraction of the mislabelled particles in the bulk parent and nucleating phases). This work supports the use of seeding to obtain fast and reasonably accurate nucleation rate estimates and the mislabelling criterion as one giving the relevant cluster size for classical nucleation theory in crystallization studies.

Physics-driven coarse-grained model for biomolecular phase separation with near-quantitative accuracy

Various physics- and data-driven sequence-dependent protein coarse-grained models have been developed to study biomolecular phase separation and elucidate the dominant physicochemical driving forces. Here, we present Mpipi, a multiscale coarse-grained model that describes almost quantitatively the change in protein critical temperatures as a function of amino-acid sequence. The model is parameterised from both atomistic simulations and bioinformatics data and accounts for the dominant role of π-π and hybrid cation-π/π-π interactions and the much stronger attractive contacts established by arginines than lysines. We provide a comprehensive set of benchmarks for Mpipi and seven other residue-level coarse-grained models against experimental radii of gyration and quantitative in-vitro phase diagrams; Mpipi predictions agree well with experiment on both fronts. Moreover, it can account for protein-RNA interactions, correctly predicts the multiphase behaviour of a charge-matched poly-arginine/poly-lysine/RNA system, and recapitulates experimental LLPS trends for sequence mutations on FUS, DDX4 and LAF-1 proteins.

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

Graphite and related sp² 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.

CμMD reveals multi-layer electrolyte screening in the double layer beyond 0.6 M, which affects ion activities, speciation and mobility; asymmetric charge screening explains concentration dependent changes to electrochemical properties.

New Features of the Open Source Monte Carlo Software Brick-CFCMC: Thermodynamic Integration and Hybrid Trial Moves

We present several new major features added to the Monte Carlo (MC) simulation code Brick-CFCMC for phase- and reaction equilibria calculations (https://gitlab.com/ETh_TU_Delft/Brick-CFCMC). The first one is thermodynamic integration for the computation of excess chemical potentials (μ^ex). For this purpose, we implemented the computation of the ensemble average of the derivative of the potential energy with respect to the scaling factor for intermolecular interactions (). Efficient bookkeeping is implemented so that the quantity is updated after every MC trial move with negligible computational cost. We demonstrate the accuracy and reliability of the calculation of μ^ex for sodium chloride in water. Second, we implemented hybrid MC/MD translation and rotation trial moves to increase the efficiency of sampling of the configuration space. In these trial moves, short Molecular Dynamics (MD) trajectories are performed to collectively displace or rotate all molecules in the system. These trajectories are accepted or rejected based on the total energy drift. The efficiency of these trial moves can be tuned by changing the time step and the trajectory length. The new trial moves are demonstrated using MC simulations of a viscous fluid (deep eutectic solvent).

The seven deadly sins: When computing crystal nucleation rates, the devil is in the details

The formation of crystals has proven to be one of the most challenging phase transformations to quantitatively model-let alone to actually understand-be it by means of the latest experimental technique or the full arsenal of enhanced sampling approaches at our disposal. One of the most crucial quantities involved with the crystallization process is the nucleation rate, a single elusive number that is supposed to quantify the average probability for a nucleus of critical size to occur within a certain volume and time span. A substantial amount of effort has been devoted to attempt a connection between the crystal nucleation rates computed by means of atomistic simulations and their experimentally measured counterparts. Sadly, this endeavor almost invariably fails to some extent, with the venerable classical nucleation theory typically blamed as the main culprit. Here, we review some of the recent advances in the field, focusing on a number of perhaps more subtle details that are sometimes overlooked when computing nucleation rates. We believe it is important for the community to be aware of the full impact of aspects, such as finite size effects and slow dynamics, that often introduce inconspicuous and yet non-negligible sources of uncertainty into our simulations. In fact, it is key to obtain robust and reproducible trends to be leveraged so as to shed new light on the kinetics of a process, that of crystal nucleation, which is involved into countless practical applications, from the formulation of pharmaceutical drugs to the manufacturing of nano-electronic devices.

Free energies of crystals computed using Einstein crystal with fixed center of mass and differing spring constants

Free energies of crystals computed using a center of mass constraint require a finite-size correction, as shown in previous work by Polson et al. [J. Chem. Phys. 112, 5339-5342 (2000)]. Their reference system is an Einstein crystal with equal spring constants. In this paper, we extend the work of Polson et al. [J. Chem. Phys. 112, 5339-5342 (2000)] to the case of differing spring constants. The generalization is convenient for constraining the center of mass in crystals with atoms of differing masses, and it helps to optimize the free energy calculations. To test the theory, we compare the free energies of LiI and NaCl crystals from calculations with differing spring constants to those computed using equal spring constants. Using these center of mass finite size corrections, we compute the true free energies of these crystals for different system sizes to eliminate the intrinsic finite-size effects. These calculations help demonstrate the size of these finite-size corrections relative to other contributions to the absolute free energy of the crystals.

Reentrant liquid condensate phase of proteins is stabilized by hydrophobic and non-ionic interactions

Liquid-liquid phase separation of proteins underpins the formation of membraneless compartments in living cells. Elucidating the molecular driving forces underlying protein phase transitions is therefore a key objective for understanding biological function and malfunction. Here we show that cellular proteins, which form condensates at low salt concentrations, including FUS, TDP-43, Brd4, Sox2, and Annexin A11, can reenter a phase-separated regime at high salt concentrations. By bringing together experiments and simulations, we demonstrate that this reentrant phase transition in the high-salt regime is driven by hydrophobic and non-ionic interactions, and is mechanistically distinct from the low-salt regime, where condensates are additionally stabilized by electrostatic forces. Our work thus sheds light on the cooperation of hydrophobic and non-ionic interactions as general driving forces in the condensation process, with important implications for aberrant function, druggability, and material properties of biomolecular condensates.

Reentrant liquid condensate phase of proteins is stabilized by hydrophobic and non-ionic interactions

Liquid-liquid phase separation of proteins underpins the formation of membraneless compartments in living cells. Elucidating the molecular driving forces underlying protein phase transitions is therefore a key objective for understanding biological function and malfunction. Here we show that cellular proteins, which form condensates at low salt concentrations, including FUS, TDP-43, Brd4, Sox2, and Annexin A11, can reenter a phase-separated regime at high salt concentrations. By bringing together experiments and simulations, we demonstrate that this reentrant phase transition in the high-salt regime is driven by hydrophobic and non-ionic interactions, and is mechanistically distinct from the low-salt regime, where condensates are additionally stabilized by electrostatic forces. Our work thus sheds light on the cooperation of hydrophobic and non-ionic interactions as general driving forces in the condensation process, with important implications for aberrant function, druggability, and material properties of biomolecular condensates.

NaCl aggregation in water at elevated temperatures and pressures: Comparison of classical force fields

The properties of water vary dramatically with temperature and density. This can be exploited to control its effectiveness as a solvent. Thus, supercritical water is of keen interest as solvent in many extraction processes. The low solubility of salts in lower density supercritical water has even been suggested as a means of desalination. The high temperatures and pressures required to reach supercritical conditions can present experimental challenges during collection of required physical property and phase equilibria data, especially in salt-containing systems. Molecular simulations have the potential to be a valuable tool for examining the behavior of solvated ions at these high temperatures and pressures. However, the accuracy of classical force fields under these conditions is unclear. We have, therefore, undertaken a parametric study of NaCl in water, comparing several salt and water models at 200 bar-600 bar and 450 K-750 K for a range of salt concentrations. We report a comparison of structural properties including ion aggregation, hydrogen bonding, density, and static dielectric constants. All of the force fields qualitatively reproduce the trends in the liquid phase density. An increase in ion aggregation with decreasing density holds true for all of the force fields. The propensity to aggregate is primarily determined by the salt force field rather than the water force field. This coincides with a decrease in the water static dielectric constant and reduced charge screening. While a decrease in the static dielectric constant with increasing NaCl concentration is consistent across all model combinations, the salt force fields that exhibit more ionic aggregation yield a slightly smaller dielectric decrement.

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

Absolute chemical potentials for complex molecules in fluid phases: A centroid reference for predicting phase equilibria

Solid-fluid phase equilibria are difficult to predict in simulations because bound degrees of freedom in the crystal phase must be converted to free translations and rotations in the fluid phase. Here, we avoid the solid-to-fluid transformation step by starting with chemical potentials for two reference systems, one for the fluid phase and one for the solid phase. For the solid, we start from the Einstein crystal and transform to the fully interacting molecular crystal. For the fluid phase, we introduce a new reference system, the "centroid," and then transform to gas phase molecules. We illustrate the new calculations by predicting the sublimation vapor pressure of succinic acid in the temperature range of 300 K-350 K.

Dielectric constant and density of aqueous alkali halide solutions by molecular dynamics: A force field assessment

The concentration dependence of the dielectric constant and the density of 11 aqueous alkali halide solutions (LiCl, NaCl, KCl, RbCl, CsCl, LiI, NaI, KI, CsI, KF, and CsF) is investigated by molecular simulation. Predictions using eight non-polarizable ion force fields combined with the TIP4P/ε water model are compared to experimental data. The influence of the water model and the temperature on the results for the NaCl brine are also addressed. The TIP4P/ε water model improves the accuracy of dielectric constant predictions compared to the SPC/E water model. The solution density is predicted well by most ion models. Almost all ion force fields qualitatively capture the decline of the dielectric constant with the increase of concentration for all solutions and with the increase of temperature for NaCl brine. However, the sampled dielectric constant is mostly in poor quantitative agreement with experimental data. These results are related to the microscopic solution structure, ion pairing, and ultimately the force field parameters. Ion force fields with excessive contact ion pairing and precipitation below the experimental solubility limit generally yield higher dielectric constant values. An adequate reproduction of the experimental solubility limit should therefore be a prerequisite for further investigations of the dielectric constant of aqueous electrolyte solutions by molecular simulation.

Lennard-Jones Parameters Determined to Reproduce the Solubility of NaCl and KCl in SPC/E, TIP3P, and TIP4P/2005 Water

Most classical nonpolarizable ion potential models underestimate the solubility values of NaCl and KCl in water significantly. We determine Lennard-Jones parameters of Na⁺, K⁺, and Cl^- that reproduce the solubility as well as the hydration free energy in dilute aqueous solutions for three water potential models, SPC/E, TIP3P, and TIP4P/2005. The ion-oxygen distance in the solution and the cation-anion distance in salt are also considered in the parametrization. In addition to the target properties, the hydration enthalpy, hydration entropy, self-diffusion coefficient, coordination number, lattice energy, enthalpy of solution, density, viscosity, and number of contact ion pairs are calculated for comparison with 17 frequently used or recently developed ion potential models. The overall performance of each ion model is represented by a global score using a scheme that was originally developed for comparison of water potential models. The global score is better for our models than for the other 17 models not only because of the quite good prediction for the solubility but also because of the relatively small deviation from the experimental value for many of the other properties.

Performing solvation free energy calculations in LAMMPS using the decoupling approach

The decoupling approach to solvation free energy calculations requires scaling the interactions between the solute and the solution with all intramolecular interactions preserved. This paper reports a new procedure that makes it possible to these calculations in LAMMPS. The procedure is tested against built-in GROMACS capabilities. The model compounds chosen to test our methodology are ethanol and biphenyl. The LAMMPS and GROMACS results obtained are in good agreement with each other. This work should help perform solvation free energy calculations in LAMMPS and/or other molecular dynamics software having no built-in functions to implement the decoupling approach.

On the transferability of ion parameters to the TIP4P/2005 water model using molecular dynamics simulations

Countless molecular dynamics studies have relied on available ion and water force field parameters to model aqueous electrolyte solutions. The TIP4P/2005 model has proven itself to be among the best rigid water force fields, whereas many of the most successful ion parameters were optimized in combination with SPC/E, TIP3P, or TIP4P/Ew water. Many researchers have combined these ions with TIP4P/2005, hoping to leverage the strengths of both parameter sets. To assess if this widely used approach is justified and to provide a guide in selecting ion parameters, we investigated the transferability of various commonly used monovalent and multivalent ion parameters to the TIP4P/2005 water model. The transferability is evaluated in terms of ion hydration free energy, hydration radius, coordination number, and self-diffusion coefficient at infinite dilution. For selected ion parameters, we also investigated density, ion pairing, chemical potential, and mean ionic activity coefficients at finite concentrations. We found that not all ions are equally transferable to TIP4P/2005 without compromising their performance. In particular, ions optimized for TIP3P water were found to be poorly transferable to TIP4P/2005, whereas ions optimized for TIP4P/Ew water provided nearly perfect transferability. The latter ions also showed good overall agreement with experimental values. The one exception is that no combination of ion parameters and water model considered here was found to accurately reproduce experimental self-diffusion coefficients. Additionally, we found that cations optimized for SPC/E and TIP3P water displayed consistent underpredictions in the hydration free energy, whereas anions consistently overpredicted the hydration free energy.

Brine rejection and hydrate formation upon freezing of NaCl aqueous solutions

Studying the freezing of saltwater on a molecular level is of fundamental importance for improving freeze desalination techniques.

Solubility prediction for a soluble organic molecule via chemical potentials from density of states

While the solubility of a substance is a fundamental property of widespread significance, its prediction from first principles (starting from only the knowledge of the molecular structure of the solute and solvent) remains a challenge. Recently, we proposed a robust and efficient method to predict the solubility from the density of states of a solute-solvent system using classical molecular simulation. The efficiency, and indeed the generality, of the method has now been enhanced by extending it to calculate solution chemical potentials (rather than probability distributions as done previously), from which solubility may be accessed. The method has been employed to predict the chemical potential of Form 1 of urea in both water and methanol for a range of concentrations at ambient conditions and for two charge models. The chemical potential calculations were validated by thermodynamic integration with the two sets of values being in excellent agreement. The solubility determined from the chemical potentials for urea in water ranged from 0.46 to 0.50 mol kg^-1, while that for urea in methanol ranged from 0.62 to 0.85 mol kg^-1, over the temperature range 298-328 K. In common with other recent studies of solubility prediction from molecular simulation, the predicted solubilities differ markedly from experimental values, reflecting limitations of current forcefields.

A force field of Li+, Na+, K+, Mg2+, Ca2+, Cl-, and SO4 2- in aqueous solution based on the TIP4P/2005 water model and scaled charges for the ions

In this work, a force field for several ions in water is proposed. In particular, we consider the cations Li⁺, Na⁺, K⁺, Mg²⁺, and Ca²⁺ and the anions Cl^- and SO₄ ^2-. These ions were selected as they appear in the composition of seawater, and they are also found in biological systems. The force field proposed (denoted as Madrid-2019) is nonpolarizable, and both water molecules and sulfate anions are rigid. For water, we use the TIP4P/2005 model. The main idea behind this work is to further explore the possibility of using scaled charges for describing ionic solutions. Monovalent and divalent ions are modeled using charges of 0.85 and 1.7, respectively (in electron units). The model allows a very accurate description of the densities of the solutions up to high concentrations. It also gives good predictions of viscosities up to 3 m concentrations. Calculated structural properties are also in reasonable agreement with the experiment. We have checked that no crystallization occurred in the simulations at concentrations similar to the solubility limit. A test for ternary mixtures shows that the force field provides excellent performance at an affordable computer cost. In summary, the use of scaled charges, which could be regarded as an effective and simple way of accounting for polarization (at least to a certain extend), improves the overall description of ionic systems in water. However, for purely ionic systems, scaled charges will not adequately describe neither the solid nor the melt.