Grand-canonical molecular dynamics simulations powered by a hybrid 4D nonequilibrium MD/MC method: Implementation in LAMMPS and applications to electrolyte solutions.

Molecular simulations in an open environment, involving ion exchange, are necessary to study various systems, from biosystems to confined electrolytes. However, grand-canonical simulations are often computationally demanding in condensed phases. A promising method [L. Belloni, J. Chem. Phys. 151, 021101 (2019)], one of the hybrid nonequilibrium molecular dynamics/Monte Carlo algorithms, was recently developed, which enables efficient computation of fluctuating number or charge density in dense fluids or ionic solutions. This method facilitates the exchange through an auxiliary dimension, orthogonal to all physical dimensions, by reducing initial steric and electrostatic clashes in three-dimensional systems. Here, we report the implementation of the method in LAMMPS with a Python interface, allowing facile access to grand-canonical molecular dynamics simulations with massively parallelized computation. We validate our implementation with two electrolytes, including a model Lennard-Jones electrolyte similar to a restricted primitive model and aqueous solutions. We find that electrostatic interactions play a crucial role in the overall efficiency due to their long-range nature, particularly for water or ion-pair exchange in aqueous solutions. With properly screened electrostatic interactions and bias-based methods, our approach enhances the efficiency of salt-pair exchange in Lennard-Jones electrolytes by approximately four orders of magnitude, compared to conventional grand-canonical Monte Carlo. Furthermore, the acceptance rate of NaCl-pair exchange in aqueous solutions at moderate concentrations reaches about 3% at the maximum efficiency.

Dielectric response of confined water films from a classical density functional theory perspective.

We re-examine the problem of the dielectric response of highly polar liquids such as water in confinement between two walls using simple two-variable density functional theory involving number and polarisation densities. In the longitudinal polarisation case where a perturbing field is applied perpendicularly to the walls, we show that the notion of the local dielectric constant, although ill-defined at a microscopic level, makes sense when coarse-graining over the typical size of a particle is introduced. The approach makes it possible to study the effective dielectric response of thin liquid films of various thicknesses in connection with the recent experiments of Fumagalli et al., [Science, 2018, 360, 1339-1342], and to discuss the notion of the interfacial dielectric constant. We argue that the observed properties as a function of slab dimensions, in particular the very low dielectric constants of the order of 2-3 measured for thin slabs of ∼1 nm thickness do not highlight any special properties of water but can be recovered for a generic polar solvent having similar particle size and the same high dielectric constant. Regarding the transverse polarisation case where the perturbing field is parallel to the walls, the associated effective dielectric constant as a function of slab dimensions reaches bulk-like values at much shorter widths than in the longitudinal case. In both cases, we find an oscillatory behaviour for slab thicknesses in the one nanometer range due to packing effects.

Coexistence of Transient Liquid Droplets and Amorphous Solid Particles in Nonclassical Crystallization of Cerium Oxalate.

Crystallization from solution often occurs via "nonclassical" routes; that is, it involves transient, non-crystalline states like reactant-rich liquid droplets and amorphous particles. However, in mineral crystals, the well-defined thermodynamic character of liquid droplets and whether they convert─or not─into amorphous phases have remained unassessed. Here, by combining cryo-transmission electron microscopy and X-ray scattering down to a 250 ms reaction time, we unveil that crystallization of cerium oxalate involves a metastable chemical equilibrium between transient liquid droplets and solid amorphous particles: contrary to the usual expectation, reactant-rich droplets do not evolve into amorphous solids. Instead, at concentrations above 2.5 to 10 mmol L-1, both amorphous and reactant-rich liquid phases coexist for several tens of seconds and their molar fractions remain constant and follow the lever rule in a multicomponent phase diagram. Such a metastable chemical equilibrium between solid and liquid precursors has been so far overlooked in multistep nucleation theories and highlights the interest of rationalizing phase transformations using multicomponent phase diagrams not only when designing and recycling rare earths materials but also more generally when describing nonclassical crystallization.

Accurate prediction of hydration free energies and solvation structures using molecular density functional theory with a simple bridge functional.

This paper assesses the ability of molecular density functional theory to predict efficiently and accurately the hydration free energies of molecular solutes and the surrounding microscopic water structure. A wide range of solutes were investigated, including hydrophobes, water as a solute, and the FreeSolv database containing 642 drug-like molecules having a variety of shapes and sizes. The usual second-order approximation of the theory is corrected by a third-order, angular-independent bridge functional. The overall functional is parameter-free in the sense that the only inputs are bulk water properties, independent of the solutes considered. These inputs are the direct correlation function, compressibility, liquid-gas surface tension, and excess chemical potential of the solvent. Compared to molecular simulations with the same force field and the same fixed solute geometries, the present theory is shown to describe accurately the solvation free energy and structure of both hydrophobic and hydrophilic solutes. Overall, the method yields a precision of order 0.5 kBT for the hydration free energies of the FreeSolv database, with a computer speedup of 3 orders of magnitude. The theory remains to be improved for a better description of the H-bonding structure and the hydration free energy of charged solutes.

Simple Parameter-Free Bridge Functionals for Molecular Density Functional Theory. Application to Hydrophobic Solvation.

Computer simulations have been fundamental in understanding the fine details of hydrophobic solvation and hydrophobic interactions. Alternative approaches based on liquid-state theories have been proposed, but are not yet at the same degree of completeness and accuracy. In this vein, a classical, molecular density functional theory approach to hydrophobic solvation is introduced. The lowest, second-order approximation of the theory, equivalent to the hypernetted chain approximation in integral equations, fails in describing correctly cavitation free-energies. It is corrected here by two simple, angular-independent, so-called bridge functionals; they are parameter-free in the sense that all variables can be fixed unambiguously from the water bulk properties, including pressure, isothermal compressibility, and liquid-gas surface tension. A hard-sphere bridge functional, based on the known functional of a reference hard fluid system, turns out to face strong limitations for water. A simpler weighted density approximation is shown to properly reproduce the solvation free energy of hydrophobes of various sizes, from microscopic ones to the nanoscale, and predicting the solvation free energy of a data set of more than 600 model hydrophobic molecules having a variety of shapes and sizes with an accuracy of a quarter of kBT compared to Monte Carlo simulations values. It constitutes an excellent starting point for a general functional describing accurately both hydrophobic and hydrophilic solvation, and making it possible to study nonidealized hydrophobic interactions.

Predicting Hydration Free Energies of the FreeSolv Database of Drug-like Molecules with Molecular Density Functional Theory.

We assess the performance of molecular density functional theory (MDFT) to predict hydration free energies of the small drug-like molecules benchmark, FreeSolv. The MDFT in the hypernetted chain approximation (HNC) coupled with a pressure correction predicts experimental hydration free energies of the FreeSolv database within 1 kcal/mol with an average computation time of 2 cpu·min per molecule. This is the same accuracy as for simulation-based free energy calculations that typically require hundreds of cpu·h or tens of gpu·h per molecule.

Hydration free energies and solvation structures with molecular density functional theory in the hypernetted chain approximation.

The capability of molecular density functional theory in its lowest, second-order approximation, equivalent to the hypernetted chain approximation in integral equations, to predict accurately the hydration free-energies and microscopic structure of molecular solutes is explored for a variety of systems: spherical hydrophobic solutes, ions, water as a solute, and the Mobley's dataset of organic molecules. The successes and the caveats of the approach are carefully pinpointed. Compared to molecular simulations with the same force field and the same fixed solute geometries, the theory describes accurately the solvation of cations, less so that of anions or generally H-bond acceptors. Overall, the electrostatic contribution to solvation free-energies of neutral molecules is correctly reproduced. On the other hand, the cavity contribution is poorly described but can be corrected using scaled-particle theory ideas. Addition of a physically motivated, one-parameter cavity correction accounting for both pressure and surface effects in the nonpolar solvation contribution yields a precision of 0.8 kcal/mol for the overall hydration free energies of the whole Mobley's dataset. Inclusion of another one-parameter cavity correction for the electrostatics brings it to 0.6 kcal/mol, that is, kBT. This is accomplished with a three-orders of magnitude numerical speed-up with respect to molecular simulations.

Structure Factor of EuCl3 Aqueous Solutions via Coupled Molecular Dynamics Simulations and Integral Equations.

Identifying the structure of an aqueous solution is essential to rationalize various phenomena such as crystallization in solution, chemical reactivity, extraction of rare earth elements, and so forth. Despite this, the efforts to describe the structure of an aqueous solution have been hindered by the difficulty to retrieve structural data both from experiments and simulations. To overcome this, first, undersaturated EuCl3 aqueous solutions of concentrations varying from 0.15 to 1.8 mol/kg were studied using X-ray scattering. Second, for the first time, the theoretical X-ray signal of 1.8 mol/kg EuCl3 aqueous solution was simulated, with precise details for the complete range of scattering vectors using coupled molecular dynamics and hypernetted chain integral equations, and satisfactorily compared with the 1.8 mol/kg experimental X-ray scattering signal. The theoretical calculations demonstrate that the experimental structure factor is dominated by Eu3+-Eu3+ correlations.

Non-equilibrium hybrid insertion/extraction through the 4th dimension in grand-canonical simulation.

The process of inserting/deleting a particle during grand-canonical Monte-Carlo (MC) simulations is investigated using a novel, original technique: the trial event is made of a short nonequilibrium molecular dynamics (MD) trajectory during which a coordinate w along a 4th dimension is added to the particle in the course of insertion/deletion and is forced to decrease from large values down to zero (for insertion) or increased from 0 up to large values (for extraction) at imposed vw velocity. The probability of acceptation of the whole MC move is controlled by the chemical potential and the external work applied during the trajectory. Contrary to the standard procedures which create/delete suddenly a particle, the proposed technique gives time to the fluid environment to relax during the gradual insertion/extraction before the acceptation decision. The reward for this expensive trial move is a gain of many orders of magnitude in the success rate. The power and wide domain of interest of this hybrid "H4D" algorithm which marries stochastic MC and nonequilibrium deterministic MD flavors are briefly illustrated with hard sphere, water, and electrolyte systems. The same approach can be easily adapted in order to measure the chemical potential of a solute particle immersed in a fluid during canonical or isobaric simulations. It then becomes an efficient application of the Jarzynski theorem for the determination of solvation free energy.

Finite-size corrections in numerical simulation of liquid water.

Monte Carlo (MC) simulations of the SPC/E liquid water model are performed at two numbers of molecules N = 100 and 512 and in canonical NVT, isobaric NPT, and grand canonical µVT ensembles. The molecular non-spherically symmetric pair distribution function g(r, Ω) (pdf) is accumulated in terms of projections gµνmnl(r) onto rotational invariants. The precisely measured differences between N values and between ensembles are in very good agreement with the theoretical predictions for the expected finite-size corrections of different origins: (1) the canonical simulation is affected by explicit corrections due to the absence of density fluctuations. Beyond the well-known shift in the long-range asymptote, all projections exhibit rich short-range contributions. These corrections vanish exactly in the isobaric ensemble provided that the pdf is constructed using the rigorous definition. (2) In the presence of dielectric discontinuity with the external medium surrounding the central box and its replica within the Ewald treatment of the Coulombic interactions, the disagreement with the prediction of the literature, already noticed recently for dipolar fluids, is confirmed in the present site-site model. (3) The implicit corrections originate from the coupling between the environment around a given particle and that around its periodic images in the neighboring cells. The recent, powerful MC/HNC mixed integral equation, which offers a complete and exact description of the molecular correlations in the whole real and Fourier spaces, enables us to quantify the observed N-dependence in the pdf projections down to the sub 10-3 levels.

What Does Second-Harmonic Scattering Measure in Diluted Electrolytes?

We derive a theoretical expression of the second harmonic scattering signal in diluted electrolytes compared with bulk water. We show that the enhancement of the signal with respect to pure water observed recently for electrolytes at very low dilution in the micromolar range is a mere manifestation of the Debye screening that makes the infinite-range dipole-dipole solvent correlations in 1/ r3 disappear as soon as the ionic concentration becomes finite. In q space, this translates into a correlation function having a well known singular behavior around q = 0, which drives the observed ionic effects. We find that the signal is independent of the ion-induced long-range behavior of the function ⟨cos ϕ( r)⟩ that has been recently discussed. We find also that the enhancement depends on the experimental geometry and occurs only for in-plane polarization detection, as observed experimentally. On the contrary, the measured isotope effect between light and heavy water cannot be fully explained.

Screened Coulombic Orientational Correlations in Dilute Aqueous Electrolytes.

The ion-induced long-range orientational order between water molecules recently observed in second harmonic scattering experiments and illustrated with large scale molecular dynamics simulations is quantitatively explained using the Ornstein-Zernike integral equation approach of liquid physics. This general effect, not specific to hydrogen-bonding solvents, is controlled by electroneutrality conditions, dipolar interactions, and dielectric+ionic screening. As expected, all numerical theories recover the well-known analytical expressions established 40 years ago.

Finite-size corrections in simulation of dipolar fluids.

Monte Carlo simulations of dipolar fluids are performed at different numbers of particles N = 100-4000. For each size of the cubic cell, the non-spherically symmetric pair distribution function g(r,Ω) is accumulated in terms of projections gmnl(r) onto rotational invariants. The observed N dependence is in very good agreement with the theoretical predictions for the finite-size corrections of different origins: the explicit corrections due to the absence of fluctuations in the number of particles within the canonical simulation and the implicit corrections due to the coupling between the environment around a given particle and that around its images in the neighboring cells. The latter dominates in fluids of strong dipolar coupling characterized by low compressibility and high dielectric constant. The ability to clean with great precision the simulation data from these corrections combined with the use of very powerful anisotropic integral equation techniques means that exact correlation functions both in real and Fourier spaces, Kirkwood-Buff integrals, and bridge functions can be derived from box sizes as small as N ≈ 100, even with existing long-range tails. In the presence of dielectric discontinuity with the external medium surrounding the central box and its replica within the Ewald treatment of the Coulombic interactions, the 1/N dependence of the gmnl(r) is shown to disagree with the, yet well-accepted, prediction of the literature.

Exact molecular direct, cavity, and bridge functions in water system.

The exact molecular bridge function of the extended simple point charge model of liquid water at room temperature is extracted from Monte Carlo (MC) simulation data. The projections gµνmnl(r) onto rotational invariants of the highly directional pair distribution function g(r,Ω) are accumulated during simulation performed with N = 512 molecules (cubic box size L ≈ 25 Å). Making intensive use of anisotropic integral equation techniques, the molecular Ornstein-Zernike equation fed with the MC data available at short distances and completed beyond L/2 with the hypernetted chain closure valid at long distances is then inverted in order to derive on the whole r range the direct correlation function cµνmnl(r), the cavity function yµνmnl(r), the negative excess potential of mean force lnyµνmnl(r), and, finally, the holy grail in such liquid state theory, the bridge function bµνmnl(r) projections. For completeness, the short distance domain inside the soft core can be reached, thanks to the use of a specially designed anisotropic finite potential which replaces the true one between a single pair of molecules in the simulation. The final bridge function b(r,Ω) of bulk water presents strong, non-universal directional features and can now serve as a reference for approximated bridge functions or functionals in liquid physics of aqueous solvents and solutions.

Efficient molecular density functional theory using generalized spherical harmonics expansions.

We show that generalized spherical harmonics are well suited for representing the space and orientation molecular density in the resolution of the molecular density functional theory. We consider the common system made of a rigid solute of arbitrary complexity immersed in a molecular solvent, both represented by molecules with interacting atomic sites and classical force fields. The molecular solvent density ρ(r,Ω) around the solute is a function of the position r≡(x,y,z) and of the three Euler angles Ω≡(θ,ϕ,ψ) describing the solvent orientation. The standard density functional, equivalent to the hypernetted-chain closure for the solute-solvent correlations in the liquid theory, is minimized with respect to ρ(r,Ω). The up-to-now very expensive angular convolution products are advantageously replaced by simple products between projections onto generalized spherical harmonics. The dramatic gain in speed of resolution enables to explore in a systematic way molecular solutes of up to nanometric sizes in arbitrary solvents and to calculate their solvation free energy and associated microscopic solvent structure in at most a few minutes. We finally illustrate the formalism by tackling the solvation of molecules of various complexities in water.

Electrolytes at interfaces: accessing the first nanometers using X-ray standing waves.

Ion-surface interactions are of high practical importance in a wide range of technological, environmental and biological problems. In particular, they ultimately control the electric double layer structure, hence the interaction between particles in aqueous solutions. Despite numerous achievements, progress in their understanding is still limited by the lack of experimental determination of the surface composition with appropriate resolution. Tackling this challenge, we have developed a method based on X-ray standing waves coupled to nano-confinement which allows the determination of ion concentrations at a solid-solution interface with a sub-nm resolution. We have investigated mixtures of KCl/CsCl and KCl/KI in 0.1 mM to 10 mM concentrations on silica surfaces and obtained quantitative information on the partition of ions between bulk and Stern layer as well as their distribution in the Stern layer. Regarding partition of potassium ions, our results are in agreement with a recent AFM study. We show that in a mixture of KCl and KI, chloride ions exhibit a higher surface propensity than iodide ions, having a higher concentration within the Stern layer and being on average closer to the surface by ≈1-2 Å, in contrast to the solution water interface. Confronting such data with molecular simulations will lead to a precise understanding of ionic distributions at aqueous interfaces.

Multipod-like silica/polystyrene clusters.

Multipod-like clusters composed of a silica core and PS satellites are prepared according to a seeded-growth emulsion polymerization of styrene in the presence of size-monodisperse silica particles previously surface-modified with methacryloxymethyltriethoxysilane. Tuning the diameter and concentration of the silica seeds affords homogeneous batches of tetrapods, hexapods, octopods, nonapods and dodecapods with morphology yields as high as 80%. Three-dimensional reconstructions by cryo-electron tomography are presented on large fields for the first time to show the high symmetry and regularity of the clusters demonstrating the good control of the synthesis process. These synthesis experiments are visited again digitally, in order to successfully refine an original simulation model and better understand the correlation between the history of the cluster growth and the final composition of the cluster mixture. Finally, using the model as a predictive tool and varying the extra experimental conditions, e.g. the composition of the surfactant mixture and the styrene concentration, result in trapping other cluster morphologies, such as tripods.

Single-step formation of micron long (OH)3Al2O3Ge(OH) imogolite-like nanotubes.

Micron-long germanium-based double-walled imogolite nanotubes were synthesized at high concentrations, as evidenced by cryo-TEM, AFM, SAXS and IR characterization methods. In addition, the spontaneous formation of a liquid-crystalline phase was observed. The novel synthesis route made it possible for the first time to obtain both long and concentrated germanium-based imogolite-like nanotubes in a single step.

X-ray reflectivity at polarized liquid-Hg-aqueous-electrolyte interface: challenging macroscopic approaches for ion-specificity issues.

We report Angstrom-resolved x-ray reflectivity analysis of externally polarized liquid-Hg surface in contact with molar LiCl, LiBr, and MgSO4 aqueous electrolytes. Interpretation of reflectivity curves demonstrates a dependence of Hg-surface layering on both applied potential and ion nature. It further highlights how interfacial polarization degree impacts electron density profiles at a molecular scale. These profiles indicate accumulation of anions and cations at the Hg surface. Upon decrease of the potential from the point of zero charge, anions are gradually expelled from the Hg surface. The study challenges traditional thermodynamic approaches for deriving countercharge composition at the Hg-electrolyte-solution interface from macroscopic Hg-surface tension data. It further dismisses the long-standing approximation that assimilates the Hg surface to a smooth, perfect chemically inert conductor with a uniformly smeared-out surface charge density.

Spheres growing on a sphere: a model to predict the morphology yields of colloidal molecules obtained through a heterogeneous nucleation route.

Through the heterogeneous nucleation of polymer nodules on a surface-modified silica particle, the high-yield achievement of hybrid colloidal molecules with a well-controlled multipod-like morphology was recently demonstrated. However, as the formation mechanism of these colloidal molecules has not been completely understood yet, some opportunities remain to reduce the tedious empirical process needed to optimize the chemical recipes. In this work, we propose a model to help understand the formation mechanism of almost pure suspensions of well-defined colloidal molecules. The outcomes of the model allow proposing probable nucleation growth scenario able to explain the experimental results. Such a model should make easier the determination of the optimal recipe parameters for a targeted morphology. The reasonably good agreements between the model and the experimental results show that the most important processes have been captured. It is thus a first step toward the rational design of large quantities of chemically prepared colloidal molecules.