Skip to main content
NCBI home page
ACS AuthorChoice logoLink to ACS AuthorChoice
. 2023 Apr 25;56(10):1156–1167. doi: 10.1021/acs.accounts.2c00811

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

Tarak Karmakar †,*, Aaron R Finney ‡, Matteo Salvalaglio ‡,*, A Ozgur Yazaydin ‡, Claudio Perego ¶,*
PMCID: PMC10193523  PMID: 37120847

Conspectus

graphic file with name ar2c00811_0006.jpg

Concentration-driven processes in solution, i.e., phenomena that are sustained by persistent concentration gradients, such as crystallization and surface adsorption, are fundamental chemical processes. Understanding such phenomena is crucial for countless applications, from pharmaceuticals to biotechnology. Molecular dynamics (MD), both in- and out-of-equilibrium, plays an essential role in the current understanding of concentration-driven processes. Computational costs, however, impose drastic limitations on the accessible scale of simulated systems, hampering the effective study of such phenomena. In particular, due to these size limitations, closed system MD of concentration-driven processes is affected by solution depletion/enrichment that unavoidably impacts the dynamics of the chemical phenomena under study. As a notable example, in simulations of crystallization from solution, the transfer of monomers between the liquid and crystal phases results in a gradual depletion/enrichment of solution concentration, altering the driving force for phase transition. In contrast, this effect is negligible in experiments, given the macroscopic size of the solution volume. Because of these limitations, accurate MD characterization of concentration-driven phenomena has proven to be a long-standing simulation challenge. While disparate equilibrium and nonequilibrium simulation strategies have been proposed to address the study of such processes, the methodologies are in continuous development.

In this context, a novel simulation technique named constant chemical potential molecular dynamics (CμMD) was recently proposed. CμMD employs properly designed, concentration-dependent external forces that regulate the flux of solute species between selected subregions of the simulation volume. This enables simulations of systems under a constant chemical drive in an efficient and straightforward way. The CμMD scheme was originally applied to the case of crystal growth from solution and then extended to the simulation of various physicochemical processes, resulting in new variants of the method. This Account illustrates the CμMD method and the key advances enabled by it in the framework of in silico chemistry. We review results obtained in crystallization studies, where CμMD allows growth rate calculations and equilibrium shape predictions, and in adsorption studies, where adsorption thermodynamics on porous or solid surfaces was correctly characterized via CμMD. Furthermore, we will discuss the application of CμMD variants to simulate permeation through porous materials, solution separation, and nucleation upon fixed concentration gradients. While presenting the numerous applications of the method, we provide an original and comprehensive assessment of concentration-driven simulations using CμMD. To this end, we also shed light on the theoretical and technical foundations of CμMD, underlining the novelty and specificity of the method with respect to existing techniques while stressing its current limitations. Overall, the application of CμMD to a diverse range of fields provides new insight into many physicochemical processes, the in silico study of which has been hitherto limited by finite-size effects. In this context, CμMD stands out as a general-purpose method that promises to be an invaluable simulation tool for studying molecular-scale concentration-driven phenomena.

Key References

  • Perego C.; Salvalaglio M.; Parrinello M.. Molecular dynamics simulations of solutions at constant chemical potential. J. Chem. Phys. 2015, 142, 144113. 10.1063/1.4917200.1 This article introduces the CμMD method, reporting its application to the test case of urea crystallization from an aqueous solution.

  • Ozcan A.; Perego C.; Salvalaglio M.; Parrinello M.; Yazaydin O.. Concentration gradient driven molecular dynamics: a new method for simulations of membrane permeation and separation. Chemical Science 2017, 8, 3858–3865 10.1039/C6SC04978H.2 This article reports the first application of asymmetric CμMD to study the permeation of organic compounds through a zeolitic imidazolate framework-8 membrane.

  • Karmakar T.; Piaggi P. M.; Parrinello M.. Molecular dynamics simulations of crystal nucleation from solution at constant chemical potential. J. Chem. Theory Comput. 2019, 15, 6923–6930 10.1021/acs.jctc.9b00795(3) This article presents the CμMD method in its spherical implementation, applied to the homogeneous nucleation of sodium chloride from an aqueous solution.

  • Finney A. R.; McPherson I. J.; Unwin P. R.; Salvalaglio M.. Electrochemistry, ion adsorption and dynamics in the double layer: a study of NaCl(aq) on graphite. Chemical Science 2021, 12, 11166–11180 10.1039/D1SC02289J(4) This article presents symmetric CμMD to investigate the adsorption of ions at carbon surfaces with insight into the concentration dependence of solution electrochemistry and activity.

1. Introduction

Molecular dynamics (MD) simulation is a useful technique for investigating in microscopic detail physicochemical and biological processes where imposing the thermodynamic conditions that control a process of interest, including temperature and pressure, is essential. In this Account, we address the key challenge of simulating out-of-equilibrium concentration-driven processes, i.e., physicochemical processes that are sustained by solution concentration gradients, subject to persistent differences in chemical potential and to the steady transfer of molecules. Whereas the simulation of these processes is of considerable interest for the study of, e.g., crystallization, transport through porous media, biomolecular condensation, it is less than straightforward. This is mainly because out-of-equilibrium simulations cannot exploit equilibrium statistical mechanics and must rely on the complexity and lower predictive capabilities of nonequilibrium theories.5,6 In addition, these simulations are most often affected by finite-size artifacts; typical MD simulations contain up to ∼107 atoms/particles, with periodic boundary conditions (PBCs) that remove surface effects. This limited size is often inadequate to study concentration-driven processes, as a constant thermodynamic force that “feeds” the process of interest is required. For example, during crystallization, solute molecules are drawn from the solution to the crystallite. In standard closed-system simulations with a fixed number of molecules, this results in solute depletion and a reduction in crystal growth rates due to a decreasing supersaturation. Therefore, comparisons with experiments, where the bulk composition is almost constant over μs time scales, are challenging at best. Whereas theoretical corrections to such finite-size effects have been proposed,7−15ad hoc simulation techniques are crucial for this purpose.

Overcoming finite-size limitations ideally requires open-system simulation, a crucial challenge in molecular modeling.16 The core feature of open-system molecular simulation is the free exchange of molecules between an explicitly modeled finite-sized region of the system, which hosts the molecular process of interest, and an external reservoir. This strategy is typically employed in equilibrium molecular simulations to sample the grand canonical (GC) ensemble, where the total number of molecules varies and the chemical potential μ is constant. This is the case in simulation methods such as GC-MD,17−20 that introduced variable-particle Hamiltonians, GC-Monte Carlo (GC-MC),21−23 or hybrid GC schemes,20,24 coupling MD or MC engines with stochastic particle insertion/deletion algorithms. While variable-particle MD schemes are still not broadly exploited, particle insertion methods have been applied more extensively, also for in silico chemical potential calculation25 and in nonequilibrium simulations (see below). Nonetheless, when simulating condensed matter systems, the inevitable very low probabilities for particle insertion can be problematic.26 This issue can be addressed by enhancing the cavity location,27 using nonequilibrium MC steps,28 or introducing extra dimensionalities.29 Recently, the effect of open boundaries was also implemented by employing a machine-learning strategy based on neural networks.30

Given the difficulties in open-system simulation, alternatives were also proposed, for example, by combining carefully prepared closed simulations, representing subsequent steps along a proper reaction coordinate, to maintain the solution concentration as the chemical process, e.g., nucleation, progresses.31

A widespread closed-system strategy to reduce finite-size effects is the multiple resolution approach in which a region of interest is described with high resolution, typically atomistic, and the remaining system is described using a lower-resolution (coarse-grained) model, which acts as a molecule reservoir for the system. While this allows the simulation of relatively large volumes,32−34 special care is required to couple these two regions. To this aim, different strategies were developed, including AdResS35 and H-AdResS,36−38 that have been successfully applied to explore a broad class of systems and phenomena, including the dynamics of polymers and biomolecules in water. Multiple resolution methods were recently coupled with insertion/deletion algorithms as well, exploiting the simplified reservoir description to overcome the low probability issues.39,40

Whereas the aforementioned schemes were typically proposed for equilibrium GC simulation, the possibility to supply/remove molecules to/from a region of interest makes these methods suitable for studying concentration-driven, out-of-equilibrium processes.16 For example, MC moves are adopted to insert/remove particles in the dual-control volume GCMD41,42 and nonequilibrium chemical potential gradient-based methods43,44 in which a chemical potential gradient is generated between the “source” (high-concentration) and “sink” (low-concentration) regions of the simulation cell. These methods have been successfully applied to study gas diffusion and separation,45−48 but in the presence of condensed matter, they still suffer from low particle insertion probabilities. In closed systems, concentration-driven simulations can also be achieved by applying external forces to establish molecular flows. This field-driven nonequilibrium approach was employed to study osmosis,49,50 diffusiophoresis,51 and transport in porous media.52,53 This extends to studies of crystallization that exploit external forces or local temperature control to establish a molecular flow facilitated by the simultaneous dissolution and growth of a crystal slab.54,55

In this Account, we focus on the constant chemical potential molecular dynamics (CμMD) technique,1 which allows us to model processes evolving out of equilibrium under the effect of concentration gradients. CμMD establishes steady-state concentration gradients in a closed system by locally controlling the concentration of the species that populate it. This is attained by a proper subdivision of the simulation volume and by the application of concentration-dependent forces, which regulate the molecular diffusion throughout the system without directly perturbing the dynamics of the process of interest. Here we discuss the theoretical foundations of CμMD and present results from its application to a wide range of systems, describing current limitations and potential developments of the method.

2. Constant Chemical Potential Molecular Dynamics (CμMD)

When a fluid phase is in contact with a phase boundary acting either as a “sink” or “source” for component i, mass transfer to/from the interphase boundary is established. By regulating the composition of the fluid phase, CμMD controls the chemical potential difference of i across the phase boundary, i.e., its mass-transfer driving force. For instance, when the sink/source is a crystal at constant temperature T and pressure P, the chemical potential difference is56

2. 1

In a system with N components, the equilibrium activity ai(xeq, T, P) is constant, so controlling N – 1 concentrations keeps the molar fraction vector x = {xi}i=1N (and with it ai) constant, consequently controlling Inline graphic(further details in the Supporting Information).57

2.1. Symmetric Scheme

CμMD was originally developed to simulate urea crystal growth in solutions with constant chemical potential (Figure 1).1 The urea concentration, ci, was assumed to be homogeneous throughout the liquid volume, except in the region close to the crystallite, where the phase transition perturbs ci from the bulk value cb, creating a concentration gradient perpendicular to the crystal–solution interface (Figure 1a). As crystallization proceeds, the phase boundary shifts and urea is depleted from the solution, lowering cb. The interface, therefore, acts as a “sink” for solute species. To maintain the driving force for crystallization, namely, the solution supersaturation, extra solute molecules must be supplied so that Inline graphic is constant (eq 1). Alternatively, if the urea crystal dissolves in an undersaturated solution, then it acts as a “source” of molecules. In this case, the driving force for dissolution can be kept constant by removing the additional solute.

Figure 1.

Figure 1

Open in a new tab

Scheme of the CμMD method. (a) Solute concentration profile ci of a planar (urea) crystal in solution. (b) CμMD action and definition of different regions of the solution volume. (See the text.) (c) Symmetric CμMD scheme applied to a urea crystal surrounded by aqueous urea solution. As in (b), the TR is highlighted in gray and the FR is in blue. (d) Urea crystal growth simulated via CμMD. The crystal size variation (molecule number, top) and molar fraction in the CR (bottom) are shown vs MD time. Constant molar fraction is enforced for several ns, enabling the growth rate estimate. This demonstrates the linear growth rate dependence on the solution composition (inset).1 Adapted with permission from ref (1) (copyright 2015, AIP).

CμMD applies external forces to control the composition of the solution phase in selected regions surrounding the phenomenon of interest. This is schematized in Figure 1b, which shows the solute concentration profile ci(z) enforced by CμMD in the presence of an interface region. As in Figure 1a, a transition region (TR) delimits the solution layer adjacent to the interface, where ci is perturbed by the sink/source effect of the phase boundary. Outside the TR, ci(z) converges to the bulk value cb. CμMD partitions the solution bulk into (i) the control region (CR), where ci is kept constant as crystallization/dissolution proceeds, (ii) the force region (FR), where the external controlling forces act, and (iii) the remaining molecule reservoir (Figure 1b). In planar crystal growth, two interfaces are simulated, as shown in Figure 1c. The instantaneous concentration of species i inside the CR is

2.1. 2

where VCR is the CR volume, N is the total number of solute particles (j), and θ(z) indicates if the position of the molecule z is inside (θ = 1) or outside (θ = 0) the CR. For practical reasons, e.g., reducing fluctuations in c, θ is defined as a continuous, differentiable switching function.1 To regulate the molecular flow toward the CR, enforcing a target ci (ci0), CμMD applies a harmonic-like force

2.1. 3

where ki is a properly chosen force constant and G(z) is a bell-shaped function that localizes Fi(z) in the FR. In planar geometry, G is chosen as

2.1. 4

where σ determines the FR width. When the CR concentration ci deviates from the desired value c0, as a result of crystallization/dissolution, the solute molecules in the FR are accelerated to balance this deviation. The user selects the location of the CRs and FRs, which species are subject to Fi, the force constants ki, and the target concentrations ci0. This scheme constitutes the core of the CμMD method, which is then conveniently adapted using system-dependent precautions for the concentration-driven phenomenon of interest. We underline that the CμMD external forces interfere with the equilibrium ensemble sampled in the entire simulation box. However, we assume that this interference does not affect the dynamics (and the desired sampling) in the active region where the process of interest takes place. To support this assumption, in the study of urea crystallization at constant chemical potential,1 geometry adaptivity was implemented: the crystal interface is located on-the-fly by monitoring the solvent distribution, and the positions of CR and FR are regulated accordingly, maintaining a fixed distance from the phase boundary. Moreover, the CR and FR geometries were appropriately set after thoroughly studying the solute and solvent concentration profiles with and without Fi, aiming at decoupling the active region from the FR and reservoir volume. Since the simulated system is finite, this decoupling holds temporarily, but this is often sufficient to estimate crystallization/dissolution rates at constant composition (Figure 1d). In general, to minimize any direct influence of Fi on the “active” region, preliminary studies to determine the optimal control (such as in ref (1)) are required. This procedure can be more or less demanding, according to the requirements of the simulation (for example, if rigorous sampling is needed for a long time) and depending on the properties of the solution (such as diffusivity). Overall, this means that the CμMD method cannot be applied in a “black box” fashion before having attained sufficient knowledge of the simulated system. However, this preparatory stage typically does not exceed the standard preliminary system assessment that MD simulations normally require (for example, at the equilibration stage), involving relatively short simulation steps.

2.1.1. Crystallization and Solubility

As mentioned above, CμMD was first applied to urea crystallization,1 investigating both the rapidly (001) and slowly growing (110) urea crystal faces.7,8 CμMD allows estimating growth rates at constant supersaturation, thus revealing crystallization trends (Figure 1d) and enabling comparisons to experiments. In another CμMD application, crystal growth in binary Lennard-Jones mixtures was simulated,58 highlighting how CμMD effectively controls the local composition of specific species near a growing crystal surface. The method was further utilized to investigate solvent effects on the crystal morphologies of an antituberculosis drug, isoniazid (INH), in alcohol solvents.59 The solvent’s role in crystal growth was investigated here, fixing the solution supersaturation and estimating the relative growth rates of the two INH surfaces. This allowed the prediction of the equilibrium, solvent-specific crystal geometry. Recently, CμMD was used to calculate relative growth rates for the (011) and (010) faces of glycine crystals in aqueous solutions and investigate the impact of ions on the crystal morphology.60

CμMD was also combined with well-tempered metadynamics61 (WTMetaD) to predict supersaturation-dependent naphthalene crystal shapes in ethanol.62 WTMetaD enhances the sampling of crystallization/dissolution events over simulation time scales by driving the system along carefully selected reaction variables, the collective variables (CV). Here, for example, a CV indicating the surface-layer crystallinity of solutes was used so that WTMetaD could efficiently sample the growth and dissolution of a new layer on the crystal. This allowed the characterization of the growth rates of the 201̅, 11̅0, and 001̅ faces of naphthalene. During enhanced growth/dissolution, concentration fluctuations were controlled by CμMD, entailing the prediction of the steady-state crystal shape. The combined WTMetaD-CμMD approach has been further applied to investigate kink-site growth at different solution concentrations, from undersaturated to supersaturated regimes.63,64 Herein, CμMD controlled the solution concentration during the calculation of free-energy differences between the crystalline and dissolved solute states.

CμMD was also successfully applied to solubility calculations. Coupling the method with WTMetaD enabled the prediction of the equilibrium solubility of molecular solids and salts.63,64 A standard symmetric CμMD scheme was instead sufficient for computing gaseous CO2 solubility in liquid water. Ansari et al. observed that standard MD of gas solubility resulted in an ∼10–20% vapor pressure depletion. CμMD alleviated this issue, providing a calculated CO2 solubility (mCO2(calc) = 0.0135) in excellent agreement with the experimental solubility (mCO2(exp) = 0.0140).65

2.1.2. Adsorption

The adsorption of fluids in porous materials is often simulated using GCMC. Particle insertion/removal in a pore region from/to a reservoir at fixed chemical potential is used here to establish the equilibrium distribution of fluids in the pores.23,66 While, GCMC does not provide information on the diffusion of molecules from the bulk phase into the pores, CμMD can provide this mechanistic insight.

Loganathan et al.67 applied CμMD to investigate CO2 and CH4 adsorption in clay mineral slit pores. They created constant composition reservoirs of equimolar CH4/CO2 gas mixtures in control volumes external to Na-montmorillonite slit pores with variable thicknesses. Simulations indicate that both the basal pore surfaces and broken edge surfaces of montmorillonite favor CO2 over CH4 adsorption, and near the aforementioned surfaces, CO2 readily displaces CH4. The preference for CO2 was due to its favorable interaction with the surfaces out to ∼20 Å. In larger pore thicknesses, however, the fluid composition far from the surfaces approached that in the external control volumes. Simulations were later extended to consider Illite slit pores with varying thicknesses, investigating the role of K+ and Na+ cations on the adsorption of CO2 and CH4 from an equimolar mixture.68 Here, basal surfaces interact more favorably with CO2, especially when K+ is the exchangeable cation, and the near-surface CO2 structures were more stable than in the Na+ case. In contrast, the protonated edge surfaces were shown to favor CH4 adsorption.

CμMD was recently adopted by Finney et al. to investigate the adsorption/depletion of Na+ and Cl– ions at graphite basal plane surfaces over a wide range of bulk aqueous solution concentrations while maintaining a constant thermodynamic driving force for adsorption.4 Here, favorable interactions between Na+ and carbon π-orbital electrons (implicitly captured by the classical force field fitted to electronic structure calculations) led to an accumulation of cations at the graphite/solution interface (Figure 2a). The asymmetric ion adsorption modified the electric potential drop on the solution side of the so-called double layer (DL) as a function of concentration, consistent with experimental measurements. At the highest concentrations, the adsorption of cations also led to the accumulation of anions in a partially saturated, multilayered interface region that emerges due to ion crowding (Figure 2a,b). The coordinate x (perpendicular to the basal surface) where ion concentrations converge marks the outermost DL edge, the analysis of which showed an inflection point in DL size as a function of bulk concentration, in contrast to mean-field predictions. This is due to overscreening and ion crowding in the DL at high concentrations, as highlighted by the electrical screening factor (f(x′) = ∫0x′cCl(x) dx/∫0cNa(x) dx) in Figure 2c, which shows deviations from the predicted monotonically increasing f as the bulk concentration increases.

Figure 2.

Figure 2

Open in a new tab

Sodium (a) and chloride (b) solution concentration profiles perpendicular to a graphite basal plane (gray bar at x ≈ 1.2 nm). Arrows and dashed lines indicate changes to the profiles as the bulk solution concentration changes, highlighted by the color scale. The rightmost gray region indicates the CμMD CR. (c) Screening factor, f, identifying how the asymmetric charge screening is affected by the bulk concentration. (d) Integral capacitance was measured in experiments (black points) and evaluated in CμMD simulations (blue points) at 1 M bulk solution concentration. (e) Excess free energy, ΔgE, for ions at the interface vs bulk concentration.57 The squares and circles represent Cl– and Na+, respectively. (a–d) Adapted with permission from ref (4) (copyright 2021, the authors). Published by the Royal Society of Chemistry under a Creative Commons Attribution 3.0 Unported License. (e) Reproduced from ref (57) (copyright 2022, Elsevier).

CμMD was essential to avoid ion depletion in the bulk and maintain electroneutral bulk solutions, which would otherwise result in unphysical electrochemical artifacts when computing the capacitance of the graphite basal surface. Finney et al. varied the surface charge density (σ) of graphite in ∼40 simulations from ±0.04 to ±1 e nm–3.4 The resulting total integral capacitance was evaluated according to C = σ/Δψ0, where Δψ0 is the potential of zero charge shifted potential drop. Excellent agreement between C evaluated in simulations and measured in experiments is apparent in Figure 2d. From these results and the analysis of the differential capacitance, Finney et al. were able to re-evaluate the long-held assumption that ion accumulation in the DL represents a minor contribution to the capacitance of graphite and demonstrated that changes in electrochemical properties are related to changes in the screening factor across the interface.

The equilibrium between solutions in the DL and the bulk (kept at controlled composition) facilitated the evaluation of the excess ion adsorption free energies (ΔgE), obtained using Maxwell–Stefan diffusive mass transfer to describe ion fluxes.57 This description of ΔgE accounts explicitly for the electric fields induced by the mobile charges in the double layer. The asymmetric driving force for ion surface adsorption is shown in Figure 2e by differences in ΔgE for Na+ and Cl–, which indicate a clear concentration dependence. Indeed, a change in the most negative ΔgE occurs at 3 to 4 mol/L, with anion adsorption around 4 kJ/mol more favorable than cation adsorption at the highest sampled concentrations due to changes in DL structure and charge screening. Small fluctuations in the CR concentrations, independent of the target, ensured accurate identification of the concentration range where this transition is observed.

2.2. Asymmetric Scheme

The symmetric CμMD scheme has been successfully applied to study crystallization, solubility, and surface adsorption in solution.1,4,62,65,69 However, symmetric CμMD is limited by the finite capacity of the reservoir in supplying or consuming species whose concentration is controlled in the CR. As a result, the constant gradient regime lasts as long as the reservoir population is sufficient to maintain the CR concentration (Figure 1d).1 Afterward, finite-size effects will interfere with the dynamics of the CRs, affecting the quantitative output of CμMD. This finite-reservoir limitation is avoided in the asymmetric variant of the CμMD, first applied to the study of fluid permeation through porous membranes2 and, subsequently, to the simultaneous investigation of the growth and dissolution of crystals in solution.70

2.2.1. Permeation Studies

The first asymmetric implementation of CμMD, dubbed concentration-gradient-driven MD (CGD-MD),2 was applied to study the permeation of fluids and separation of mixtures through porous membranes. Herein, external forces given by eq 3 are applied asymmetrically at the inlet and at the outlet sides of a porous membrane (Figure 3a) to induce a constant concentration gradient across the membrane. As a result, the solute molecules flow through the membrane toward lower concentrations while PBCs allow solute recirculation and prevent depletion at the inlet, analogous to the dual control volume GCMD method.41 Herein, CμMD generates both a source and a sink of solute molecules to investigate the permeability and permselectivity of porous materials. By allowing accurate composition control of the gaseous/liquid phase at the inlet and outlet of a membrane (Figure 3a), via modulation of molecule recirculation across PBCs, the CGD-MD scheme enables the simulation of fluid mixture permeation and separation, without the necessity of particle insertion/deletion algorithms. This provides a valid alternative to the existing methodologies.41,53,71

Figure 3.

Figure 3

Open in a new tab

CμMD for permeation studies. (a) CGD-MD scheme: a concentration gradient is generated between the inlet (left) and the outlet (right) sides of a porous membrane. PBCs allow recirculating solute from the outlet to the inlet. (b) Variation of methane flux (Jz) across a ZIF-8 membrane as a function of the concentration gradient Δc.2 (a) and (b) Adapted with permission from ref (2) (copyright 2017, the authors). Published by the Royal Society of Chemistry under a Creative Commons Attribution 3.0 Unported License. (c) Concentration profiles of methane (top) and hydrogen (bottom) across a composite PIM-1/ZIF-8 membrane in single-component and mixture simulations. Dashed lines indicate the interfaces of the PIM-1/ZIF-8 framework.72 Adapted with permission from ref (72) (copyright 2020 American Chemical Society). (d) Concentration profiles of methanol, ethanol, and water across an MFI membrane (in transparency), obtained via CGD-MD simulations.73 Adapted with permission from ref (73) (copyright 2021, Elsevier).

Ozcan et al. applied CGD-MD to evaluate the permeabilities of methane, ethane, and ethylene across a ZIF-8 membrane. The approach demonstrated that the methane flux Jz across ZIF-8 strongly increases with the Δc, reaching a plateau beyond a threshold (Figure 3b), which separates transport control by concentration difference and mass transfer resistance. Remarkably, the CμMD-computed permselectivities of the ZIF-8 membrane for methane/ethane/ethylene separation follow the same trend of the experiments. While CGD-MD performs well when the target concentrations correspond to medium/high pressure, at low gas concentrations (∼1 atm) the instantaneous concentration in the CRs undergoes unavoidably large fluctuations, reducing the effectiveness of the method. This finite-size effect restricts the applicability range of CGD-MD as reported here; very large simulation volumes or reformulating the external forces could alleviate these issues.

Namsani et al.74 recently demonstrated that CGD-MD could be reliably used to simulate ternary mixture separation in porous membranes. Here, the separation of H2 in a ZIF-8 membrane from a H2, N2, and CO2 effluent syngas mixture, at the typical composition obtained from commercial production (air-blown autothermal reforming and water–gas shift reactions), was studied. The CGD-MD simulation at the feed temperature and pressure of a typical H2 purification unit (300 K and 35 atm) showed that the ZIF-8 membrane could select H2 over N2 and CO2 species, in agreement with experiments. These results contrasted with the widely used ideal permselectivity predictions, computed using self-diffusion coefficients and infinite dilution solubility factors. The study concluded that such approximate methods are inaccurate when mixtures of strongly associating gases or high-pressure ranges are considered. Methods that allow direct calculation of permeabilities, such as CGD-MD, are therefore crucial in these cases.

Ozcan et al. simulated the H2/CH4 transport and separation across the composite PIM-1/ZIF-8 membrane. The flux of H2 and CH4 gases and the H2/CH4 permselectivity in the composite PIM-1/ZIF-8 membrane (Figure 3c) were compared to results obtained for the individual PIM-1 and ZIF-8 membranes.72 The reported H2/CH4 permselectivity is reduced by ∼20% compared to the ideal permselectivity estimated via a macroscopic model75 using the data obtained for ZIF-8 and PIM-1 membranes separately. This demonstrated that macroscopic models, which compute the permeability of composite membranes by combining the permeabilities of individual components, can be inaccurate. These models do not account for interfacial effects and defects, such as nonselective void spaces present between the polymer and crystal membranes, which affect the permeability of the composite system. An approach such as CGD-MD, combined with an accurate description of the polymer/crystal interface, can directly quantify these interface-driven deviations from ideal permselectivty.

CGD-MD can also be used to study liquid-phase membrane separations. In ref (73), the permeation of alcohol/water mixtures through a silicalite (MFI) membrane was addressed via CμMD. Pure ethanol and methanol can diffuse through the membrane, while pure water diffusion is prevented (Figure 3d) due to hydrogen bonding with the interfacial silanol groups of MFI. When instead water/ethanol or methanol mixtures were tested, water permeation was observed. This is attributed to the competition of alcohol species against MFI silanol groups in forming hydrogen bonds with water, as the bonding of water with alcohol can enable the diffusion of both species across the silicalite membrane.

2.2.2. “Cannibalistic” Crystallization

The asymmetric CμMD scheme was applied to simulate simultaneous growth and dissolution on the two opposite sides of a crystal slab, in an approach dubbed cannibalistic, referring to the feeding of the crystallization by means of the molecule of the very same crystallite. As done in the symmetric case (Figure 1c), also here a solid slab is exposed to two fluid phases; however, as in the CGD-MD scheme (Figure 3), the target concentrations of the two CR regions are different: one side of a crystal slab is exposed to a supersaturated solution leading to growth, while the other side is exposed to an undersaturated solution which favors the dissolution of the slab (Figure 4a).70 These two CR regions thus remove and supply solutes to the same reservoir region connected by PBCs. A proper balance between the growth and dissolution rates ensures that the reservoir can maintain the user-defined target concentrations irrespective of the total simulation time.

Figure 4.

Figure 4

Open in a new tab

(a) Cannibalistic CμMD setup. The solution in the left CR is supersaturated (S), and in the right CR, it is undersaturated (U). The two crystal–solution interface positions, zIS and zIU, are adaptively identified at every time step. (b) Growth profiles from three independent CμMD simulations at concentrations of 1.5, 2.0, and 2.4 nm–3. (c) Representative concentration profile: running averages demonstrate the effective control of the solution concentrations compared with target values. Adapted with permission from ref (70) (copyright 2018, American Chemical Society).

To check the effectiveness of this approach, we applied it to study the paradigmatic case of urea growth and dissolution in an aqueous solution.70 A series of independent CμMD simulations at different combinations of solution super- and undersaturation in the CRs were carried out.70 In each simulation, we observed multiple events of crystal layer growth and dissolution. The CR concentrations on both supersaturated and undersaturated sides remained close to their corresponding target values (Figure 4b,c) throughout the simulations, indicating that our method can maintain a constant chemical potential environment for an extended period of simulation. As a result, we could accurately calculate the growth and dissolution rates. It is worth noting that our cannibalistic method differs from the approach of Kusalik et al.54,55 in that CμMD implements solution concentration control, whereas temperature changes in direct forces are applied in that scheme.

2.3. Spherical Scheme: Nucleation

Even more than in crystal growth, finite-size confinement effects have the potential to suppress nucleation by radically changing the thermodynamics of the process compared to a macroscopic open system.9−14 The idea of the CμMD simulation of crystallization has therefore been extended to the study of crystal nucleation. Unlike in the case of growth and adsorption simulations, in the case of nucleation, a more convenient spherical setup would allow the concentration to be controlled as a function of the radial distance (ri) from the center of nascent clusters of an emerging phase, in line with predictions (following classical nucleation theory) that nuclei are pseudospherical objects. In this setup, a cubic simulation box is considered. From its center, a spherical region is defined, which acts as the nucleus growth region (GR) (schematic in Figure 5a). Surrounding the GR, a set of concentric shells, which can be compared with the rectangular slabs in the planar CμMD setup (Figure 1), are defined to calculate the concentration of solute species (using eqs 2 and 3, the z coordinates should be replaced with ri). Similar to the planar model, spherical volume elements representing the transition region (TR), control region (CR), and force region (FR), respectively, are defined. The region beyond the FR is the molecular reservoir, which encapsulates the remaining region of the cubic cell with periodic boundary conditions applied in all directions.

Figure 5.

Figure 5

Open in a new tab

Spherical CμMD of NaCl nucleation from aqueous solution: (a) Spherical CμMD scheme. The GR is shown in pink, followed by a thin TR, and CR is shown in green. The FR surrounds the CR. The rest of the box acts as a molecular reservoir. A crystal nucleus is depicted in the GR (red spheres). (b) Solution concentration in the CR shell vs simulation time, obtained from standard NVT metadynamics, and (c) in the CμMD setup with two different target concentrations (upper and lower panels). The color indicates the nucleus size. The corresponding free-energy surfarces, as a function of the selected CVs (sO and sH'), are depicted in (d)–(f); a representative critical nucleus and a rock-salt structure are shown in panel (f). Adapted with permission from ref (3). Copyright 2019 American Chemical Society.

This CμMD spherical variant was applied to investigate the nucleation of sodium chloride from an aqueous supersaturated solution.3 For comparison, one simulation in the standard NVT setup was carried out, followed by two CμMD simulations at two target solution concentrations, 4.2 nm–3 (molality 9.6 m) and 5.0 nm–3 (11.6 m). Note that the equilibrium solubility calculated using the Joung–Cheatham potential model76 is 3.7 m,77,78 and thus the two concentrations correspond to supersaturations, (S = m/meq) 2.6 and 3.1, respectively.

Since nucleation involves high free-energy barriers at moderate supersaturation levels, WTMetaD was used in both the NVT and CμMD simulations to enhance the sampling of NaCl crystal states on the simulation time scale. Two CVs were used in the WTMetaD simulations. The first one (sO) was based on the environment-similarity CV developed by Piaggi and Parrinello,79 which describes the local ordering of atoms in a given crystal structure. The second CV (sH') was the average number of water molecules surrounding Na+ ions within a given cutoff radius. Furthermore, these two CVs were modified to ensure that a nucleus forms inside the GR. Thanks to WTMetaD, multiple cycles of nucleation and dissolution were observed in the simulations. The concentration profiles obtained from a standard NVT (Figure 5a) and the two WTMetaD-CμMD simulations (Figure 5b) are shown. In the NVT simulation, each nucleation cycle is accompanied by a significant drop in the solution concentration in the CR. In contrast, the CR solution concentration in WTMetaD-CμMD simulations fluctuates around the target value during all crystal nucleation and dissolution events. Good control of the solution concentration during each nucleation cycle allowed us to accurately calculate nucleation free energies (Figure 5c–e) and to explicitly estimate supersaturation-dependent nucleation rates.3

3. Future Perspectives

We have herein summarized the main features and applications of the CμMD method. In the landscape of different techniques that have been proposed to implement concentration gradient simulations (Introduction section), CμMD has distinguished itself for its versatility and relatively facile application. The method takes advantage of the adaptive control of solution concentration, which provides a crucial vantage point in simulating the out-of-equilibrium evolution of many concentration-driven processes such as permeation, crystallization, adsorption, dissolution, etc. Originally developed to study crystal growth, CμMD has evolved into a multipurpose method. It is now regularly applied to study various physicochemical processes driven by concentration gradients beyond its original purpose, including nucleation and assembly, adsorption, membrane permeation, and separation. Following the applications discussed herein, CμMD could be applied to studying liquid–liquid transitions or the effect of solution environments (e.g., pH, reactant/product concentrations) on chemically reactive systems. Asymmetric CμMD generates a concentration gradient, enforcing the proper conditions to study permeation and crystal growth. This strategy can be adapted to probe concentration gradients typical of biological or electrochemical systems at nano/microscales and can be coupled with adaptive methods for the calculation of interfacial charges.80 While the spherical variant of CμMD is well-suited to investigate homogeneous nucleation in systems containing fast-diffusing solutes, this scheme could be easily adapted to probe the self-assembly of soft materials. As already demonstrated in crystallization studies, CμMD is easily coupled to enhanced sampling techniques such as Metadynamics, indicating a pathway for extending the validity range and potential of CμMD in future applications.

The future development of CμMD could also involve multiresolution schemes, such as the adaptive resolution technique of Kremer et al., allowing the simulation of very large systems with controlled compositions.38 Overall, CμMD applications are already wide-ranging, and exciting future developments will benefit simulation communities across multiple disciplines.

Acknowledgments

We acknowledge Michele Parrinello for his input during the development of the method and in the follow-up work. We thank Pablo M. Piaggi, Zoran Bjelobrk, and Dandan Han for contributing to various applications. T.K. acknowledges IIT Delhi for the institute seed grant and the Science and Engineering Research Board (SERB), New Delhi, India for the start-up research grant (file no. SRG/2022/000969). M.S. and A.R.F. acknowledge funding from the Crystallization in the Real World EPSRC Programme Grant (grant EP/R018820/1). M.S. acknowledges support from the EPSRC via the UKRI Frontier Research Guaranteee Grant (EP/X033139/1).

Biographies

Tarak Karmakar is currently an assistant professor in the Department of Chemistry at the Indian Institute of Technology, Delhi, India. He was a postdoctoral fellow in the group of Prof. Michele Parrinello at ETH Zurich, USI Lugano, Switzerland, and IIT Genova, Italy from 2017 to 2021. He received his Ph.D. degree from JNCASR, India. His current research focuses on the development and applications of advanced MD simulations and machine learning methods to study soft condensed matter and biological systems.

Aaron R. Finney is a senior research fellow in the Department of Chemical Engineering, University College London and a member of the Thomas Young Centre. He completed his Ph.D. in the Centre for Scientific Computing and Department of Chemistry, University of Warwick. Following this, he was a Doctoral Prize Fellow and Research Fellow in the Materials Science Department, University of Sheffield. His current research focuses on interfaces and crystallization.

Matteo Salvalaglio is an associate professor of chemical engineering at University College London and a member of the Thomas Young Centre and the Sargent Centre for Process Systems Engineering. He received his Ph.D. from Politecnico di Milano (Italy). Before joining UCL, he was a postdoctoral researcher at ETH Zurich (CH) in a joint appointment between chemistry and mechanical engineering. His current research focuses on using computational methods to understand nucleation, polymorphism, and crystal growth.

A. Ozgur Yazaydin is an associate professor of chemical engineering at University College London and a member of the Thomas Young Centre. He received his Ph.D. in chemical engineering from Worcester Polytechnic Institute in Massachusetts, United States. His research focuses on the computational design and multiscale modeling of nanostructured materials for applications in energy and environment.

Claudio Perego is a tenured researcher in the Department of Innovative Technologies at SUPSI, Lugano, Switzerland. Previously, he was a postdoctoral researcher in the polymer theory group at Max-Planck Institute, Mainz, Germany and in the Parrinello research group at ETH Zurich/USI Lugano, Switzerland. He received his Ph.D. from the University of Milano-Bicocca, Milan, Italy. His current research concerns the modeling of self-assembled supramolecular materials, with a focus on their in- and out-of-equilibrium dynamics.

Data Availability Statement

The original CμMD method and its variants have been implemented in the open-source software PLUMED development versions. The CμMD code, input files, and related tutorials are available from the authors’ GitHub pages (https://github.com/mme-ucl/CmuMD,https://github.com/Tarakk/plumed-cumd), materials clouds archive (10.24435/materialscloud:2020.0015/v1 and 10.24435/materialscloud:2020.0013/v1,10.24435/materialscloud:k5-t2), and PLUMED-NEST81 portal (IDs). A tutorial on applying symmetric and asymmetric planar symmetry CμMD is available on the PLUMED website as a part of the 2022 PLUMED Masterclass (https://www.plumed.org/doc-master/user-doc/html/masterclass-22-8.html).

Supporting Information Available

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.accounts.2c00811.

  • Theoretical discussion on concentration-driven processes and solution composition control (PDF)

Author Contributions

CRediT: Tarak Karmakar conceptualization (equal), writing-original draft (equal), writing-review & editing (equal); Aaron R. Finney conceptualization (equal), writing-original draft (equal), writing-review & editing (equal); Matteo Salvalaglio conceptualization (equal), writing-original draft (equal), writing-review & editing (equal); A. Ozgur Yazaydin conceptualization (equal), writing-original draft (equal), writing-review & editing (equal); Claudio Perego conceptualization (equal), writing-original draft (equal), writing-review & editing (equal).

The authors declare no competing financial interest.

Supplementary Material

ar2c00811_si_001.pdf (130KB, pdf)

References

  1. Perego C.; Salvalaglio M.; Parrinello M. Molecular dynamics simulations of solutions at constant chemical potential. J. Chem. Phys. 2015, 142, 144113. 10.1063/1.4917200. [DOI] [PubMed] [Google Scholar]
  2. Ozcan A.; Perego C.; Salvalaglio M.; Parrinello M.; Yazaydin O. Concentration gradient driven molecular dynamics: a new method for simulations of membrane permeation and separation. Chemical Science 2017, 8, 3858–3865. 10.1039/C6SC04978H. [DOI] [PMC free article] [PubMed] [Google Scholar]
  3. Karmakar T.; Piaggi P. M.; Parrinello M. Molecular dynamics simulations of crystal nucleation from solution at constant chemical potential. J. Chem. Theory Comput. 2019, 15, 6923–6930. 10.1021/acs.jctc.9b00795. [DOI] [PubMed] [Google Scholar]
  4. Finney A. R.; McPherson I. J.; Unwin P. R.; Salvalaglio M. Electrochemistry, ion adsorption and dynamics in the double layer: a study of NaCl(aq) on graphite. Chemical Science 2021, 12, 11166–11180. 10.1039/D1SC02289J. [DOI] [PMC free article] [PubMed] [Google Scholar]
  5. Jarzynski C. Nonequilibrium Equality for Free Energy Differences. Phys. Rev. Lett. 1997, 78, 2690–2693. 10.1103/PhysRevLett.78.2690. [DOI] [Google Scholar]
  6. Husic B. E.; Pande V. S. Markov State Models: From an Art to a Science. J. Am. Chem. Soc. 2018, 140, 2386–2396. 10.1021/jacs.7b12191. [DOI] [PubMed] [Google Scholar]
  7. Salvalaglio M.; Vetter T.; Giberti F.; Mazzotti M.; Parrinello M. Uncovering molecular details of urea crystal growth in the presence of additives. J. Am. Chem. Soc. 2012, 134, 17221–17233. 10.1021/ja307408x. [DOI] [PubMed] [Google Scholar]
  8. Salvalaglio M.; Vetter T.; Mazzotti M.; Parrinello M. Controlling and predicting crystal shapes: The case of urea. Angew. Chem., Int. Ed. 2013, 52, 13369–13372. 10.1002/anie.201304562. [DOI] [PubMed] [Google Scholar]
  9. Salvalaglio M.; Perego C.; Giberti F.; Mazzotti M.; Parrinello M. Molecular-dynamics simulations of urea nucleation from aqueous solution. Proc. Natl. Acad. Sci. U.S.A. 2015, 112, E6–E14. 10.1073/pnas.1421192111. [DOI] [PMC free article] [PubMed] [Google Scholar]
  10. Salvalaglio M.; Tiwary P.; Maggioni G. M.; Mazzotti M.; Parrinello M. Overcoming time scale and finite size limitations to compute nucleation rates from small scale well tempered metadynamics simulations. J. Chem. Phys. 2016, 145, 211925. 10.1063/1.4966265. [DOI] [PubMed] [Google Scholar]
  11. Agarwal V.; Peters B. Nucleation near the eutectic point in a Potts-lattice gas model. J. Chem. Phys. 2014, 140, 084111. 10.1063/1.4865338. [DOI] [PubMed] [Google Scholar]
  12. Reguera D.; Bowles R.; Djikaev Y.; Reiss H. Phase transitions in systems small enough to be clusters. J. Chem. Phys. 2003, 118, 340–353. 10.1063/1.1524192. [DOI] [Google Scholar]
  13. Wedekind J.; Reguera D.; Strey R. Finite-size effects in simulations of nucleation. J. Chem. Phys. 2006, 125, 214505. 10.1063/1.2402167. [DOI] [PubMed] [Google Scholar]
  14. Grossier R.; Veesler S. Reaching one single and stable critical cluster through finite-sized systems. Cryst. Growth Des. 2009, 9, 1917–1922. 10.1021/cg801165b. [DOI] [Google Scholar]
  15. Sevilla M.; Cortes-Huerto R. Connecting Density Fluctuations and irkwood–Buff Integrals for Finite-Size Systems. J. Chem. Phys. 2022, 156, 044502. 10.1063/5.0076744. [DOI] [PubMed] [Google Scholar]
  16. Ciccotti G.; Delle Site L. The physics of open systems for the simulation of complex molecular environments in soft matter. Soft Matter 2019, 15, 2114–2124. 10.1039/C8SM02523A. [DOI] [PubMed] [Google Scholar]
  17. Çağin T.; Pettitt B. M. Molecular Dynamics with a Variable Number of Molecules. Mol. Phys. 1991, 72, 169–175. 10.1080/00268979100100111. [DOI] [Google Scholar]
  18. Lo C.; Palmer B. Alternative Hamiltonian for molecular dynamics simulations in the grand canonical ensemble. J. Chem. Phys. 1995, 102, 925–931. 10.1063/1.469159. [DOI] [Google Scholar]
  19. Eslami H.; Müller-Plathe F. Molecular dynamics simulation in the grand canonical ensemble. Journal of computational chemistry 2007, 28, 1763–1773. 10.1002/jcc.20689. [DOI] [PubMed] [Google Scholar]
  20. Yang X.-D.; Chen W.; Ren Y.; Chu L.-Y. Exploration of the Adsorption Kinetics of Surfactants at the Water–Oil Interface via Grand-Canonical Molecular Dynamics Simulations. Langmuir 2022, 38, 1277–1286. 10.1021/acs.langmuir.1c03205. [DOI] [PubMed] [Google Scholar]
  21. Adams D. Grand Canonical Ensemble Monte Carlo for a Lennard-Jones Fluid. Mol. Phys. 1975, 29, 307–311. 10.1080/00268977500100221. [DOI] [Google Scholar]
  22. Papadopoulou A.; Becker E. D.; Lupkowski M.; van Swol F. Molecular dynamics and Monte Carlo simulations in the grand canonical ensemble: Local versus global control. J. Chem. Phys. 1993, 98, 4897–4908. 10.1063/1.464945. [DOI] [Google Scholar]
  23. Wang S.; Feng Q.; Javadpour F.; Hu Q.; Wu K. Competitive adsorption of methane and ethane in montmorillonite nanopores of shale at supercritical conditions: A grand canonical Monte Carlo simulation study. Chem. Eng. J. 2019, 355, 76–90. 10.1016/j.cej.2018.08.067. [DOI] [Google Scholar]
  24. Boinepalli S.; Attard P. Grand Canonical Molecular Dynamics. J. Chem. Phys. 2003, 119, 12769–12775. 10.1063/1.1629079. [DOI] [Google Scholar]
  25. Daly K. B.; Benziger J. B.; Debenedetti P. G.; Panagiotopoulos A. Z. Massively Parallel Chemical Potential Calculation on Graphics Processing Units. Comput. Phys. Commun. 2012, 183, 2054–2062. 10.1016/j.cpc.2012.05.006. [DOI] [Google Scholar]
  26. Delgado-Buscalioni R.; Coveney P. V. USHER: An Algorithm for Particle Insertion in Dense Fluids. J. Chem. Phys. 2003, 119, 978–987. 10.1063/1.1579475. [DOI] [Google Scholar]
  27. Delgado-Buscalioni R.; De Fabritiis G.; Coveney P. V. Determination of the Chemical Potential Using Energy-Biased Sampling. J. Chem. Phys. 2005, 123, 054105. 10.1063/1.2000244. [DOI] [PubMed] [Google Scholar]
  28. Ross G. A.; Rustenburg A. S.; Grinaway P. B.; Fass J.; Chodera J. D. Biomolecular Simulations under Realistic Macroscopic Salt Conditions. J. Phys. Chem. B 2018, 122, 5466–5486. 10.1021/acs.jpcb.7b11734. [DOI] [PMC free article] [PubMed] [Google Scholar]
  29. Belloni L. Non-Equilibrium Hybrid Insertion/Extraction through the 4th Dimension in Grand-Canonical Simulation. J. Chem. Phys. 2019, 151, 021101. 10.1063/1.5110478. [DOI] [PubMed] [Google Scholar]
  30. Floyd J. E.; Lukes J. R. A Neural Network-Assisted Open Boundary Molecular Dynamics Simulation Method. J. Chem. Phys. 2022, 156, 184114. 10.1063/5.0083198. [DOI] [PubMed] [Google Scholar]
  31. Liu C.; Wood G. P.; Santiso E. E. Modelling nucleation from solution with the string method in the osmotic ensemble. Mol. Phys. 2018, 116, 2998–3007. 10.1080/00268976.2018.1482016. [DOI] [Google Scholar]
  32. Praprotnik M.; Site L. D.; Kremer K. Multiscale simulation of soft matter: From scale bridging to adaptive resolution. Annu. Rev. Phys. Chem. 2008, 59, 545–571. 10.1146/annurev.physchem.59.032607.093707. [DOI] [PubMed] [Google Scholar]
  33. Praprotnik M.; Cortes-Huerto R.; Potestio R.; Delle Site L. Adaptive resolution molecular dynamics technique. Handbook of materials modeling: methods: theory and modeling 2018, 1–15. 10.1007/978-3-319-42913-7_89-1. [DOI] [Google Scholar]
  34. Cortes-Huerto R.; Praprotnik M.; Kremer K.; Delle Site L. From Adaptive Resolution to Molecular Dynamics of Open Systems. European Physical Journal B 2021, 94, 189. 10.1140/epjb/s10051-021-00193-w. [DOI] [PMC free article] [PubMed] [Google Scholar]
  35. Wang H.; Hartmann C.; Schütte C.; Delle Site L. Grand-canonical-like molecular-dynamics simulations by using an adaptive-resolution technique. Physical Review X 2013, 3, 011018. 10.1103/PhysRevX.3.011018. [DOI] [Google Scholar]
  36. Potestio R.; Fritsch S.; Espanol P.; Delgado-Buscalioni R.; Kremer K.; Everaers R.; Donadio D. Hamiltonian adaptive resolution simulation for molecular liquids. Physical review letters 2013, 110, 108301. 10.1103/PhysRevLett.110.108301. [DOI] [PubMed] [Google Scholar]
  37. Potestio R.; Español P.; Delgado-Buscalioni R.; Everaers R.; Kremer K.; Donadio D. Monte Carlo Adaptive Resolution Simulation of Multicomponent Molecular Liquids. Phys. Rev. Lett. 2013, 111, 060601. 10.1103/PhysRevLett.111.060601. [DOI] [PubMed] [Google Scholar]
  38. Tarenzi T.; Calandrini V.; Potestio R.; Carloni P. Open-boundary molecular mechanics/coarse-grained framework for simulations of low-resolution G-protein-coupled receptor–ligand complexes. J. Chem. Theory Comput. 2019, 15, 2101–2109. 10.1021/acs.jctc.9b00040. [DOI] [PMC free article] [PubMed] [Google Scholar]
  39. Heidari M.; Kremer K.; Golestanian R.; Potestio R.; Cortes-Huerto R. Open-Boundary Hamiltonian Adaptive Resolution. From Grand Canonical to Non-Equilibrium Molecular Dynamics Simulations. J. Chem. Phys. 2020, 152, 194104. 10.1063/1.5143268. [DOI] [PubMed] [Google Scholar]
  40. Ebrahimi Viand R.; Höfling F.; Klein R.; Delle Site L. Theory and Simulation of Open Systems out of Equilibrium. J. Chem. Phys. 2020, 153, 101102. 10.1063/5.0014065. [DOI] [PubMed] [Google Scholar]
  41. Heffelfinger G. S.; Swol F. v. Diffusion in Lennard-Jones fluids using dual control volume grand canonical molecular dynamics simulation (DCV-GCMD). J. Chem. Phys. 1994, 100, 7548–7552. 10.1063/1.466849. [DOI] [Google Scholar]
  42. Thompson A. P.; Ford D. M.; Heffelfinger G. S. Direct molecular simulation of gradient-driven diffusion. J. Chem. Phys. 1998, 109, 6406–6414. 10.1063/1.477284. [DOI] [Google Scholar]
  43. Maginn E. J.; Bell A. T.; Theodorou D. N. Transport diffusivity of methane in silicalite from equilibrium and nonequilibrium simulations. J. Phys. Chem. 1993, 97, 4173–4181. 10.1021/j100118a038. [DOI] [Google Scholar]
  44. Cracknell R. F.; Nicholson D.; Quirke N. Direct molecular dynamics simulation of flow down a chemical potential gradient in a slit-shaped micropore. Physical review letters 1995, 74, 2463. 10.1103/PhysRevLett.74.2463. [DOI] [PubMed] [Google Scholar]
  45. Pohl P. I.; Heffelfinger G. S.; Smith D. M. Molecular dynamics computer simulation of gas permeation in thin silicalite membranes. Mol. Phys. 1996, 89, 1725–1731. 10.1080/00268979609482570. [DOI] [Google Scholar]
  46. Xu L.; Sedigh M. G.; Sahimi M.; Tsotsis T. T. Nonequilibrium molecular dynamics simulation of transport of gas mixtures in nanopores. Physical review letters 1998, 80, 3511. 10.1103/PhysRevLett.80.3511. [DOI] [PubMed] [Google Scholar]
  47. Wu T.; Firoozabadi A. Molecular simulations of binary gas mixture transport and separation in slit nanopores. J. Phys. Chem. C 2018, 122, 20727–20735. 10.1021/acs.jpcc.8b04976. [DOI] [Google Scholar]
  48. Fayaz-Torshizi M.; Xu W.; Vella J. R.; Marshall B. D.; Ravikovitch P. I.; Mülller E. A. Use of Boundary-Driven Nonequilibrium Molecular Dynamics for Determining Transport Diffusivities of Multicomponent Mixtures in Nanoporous Materials. J. Phys. Chem. B 2022, 126, 1085–1100. 10.1021/acs.jpcb.1c09159. [DOI] [PMC free article] [PubMed] [Google Scholar]
  49. Yoshida H.; Marbach S.; Bocquet L. Osmotic and Diffusio-Osmotic Flow Generation at High Solute Concentration. II. Molecular Dynamics Simulations. J. Chem. Phys. 2017, 146, 194702. 10.1063/1.4981794. [DOI] [PubMed] [Google Scholar]
  50. Ganti R.; Liu Y.; Frenkel D. Hamiltonian Transformation to Compute Thermo-osmotic Forces. Phys. Rev. Lett. 2018, 121, 068002. 10.1103/PhysRevLett.121.068002. [DOI] [PubMed] [Google Scholar]
  51. Ramírez-Hinestrosa S.; Yoshida H.; Bocquet L.; Frenkel D. Studying polymer diffusiophoresis with non-equilibrium molecular dynamics. J. Chem. Phys. 2020, 152, 164901. 10.1063/5.0007235. [DOI] [PubMed] [Google Scholar]
  52. Arya G.; Chang H.-C.; Maginn E. J. A Critical Comparison of Equilibrium, Non-Equilibrium and Boundary-Driven Molecular Dynamics Techniques for Studying Transport in Microporous Materials. J. Chem. Phys. 2001, 115, 8112–8124. 10.1063/1.1407002. [DOI] [Google Scholar]
  53. Frentrup H.; Avendaño C.; Horsch M.; Salih A.; Müller E. A. Transport diffusivities of fluids in nanopores by non-equilibrium molecular dynamics simulation. Mol. Simul. 2012, 38, 540–553. 10.1080/08927022.2011.636813. [DOI] [Google Scholar]
  54. Vatamanu J.; Kusalik P. Molecular dynamics methodology to investigate steady-state heterogeneous crystal growth. J. Chem. Phys. 2007, 126, 124703. 10.1063/1.2710263. [DOI] [PubMed] [Google Scholar]
  55. Liang S.; Hall K. W.; Laaksonen A.; Zhang Z.; Kusalik P. G. Characterizing key features in the formation of ice and gas hydrate systems. Philosophical Transactions of the Royal Society A 2019, 377, 20180167. 10.1098/rsta.2018.0167. [DOI] [PMC free article] [PubMed] [Google Scholar]
  56. Prausnitz J. M.; Lichtenthaler R. N.; De Azevedo E. G.. Molecular Thermodynamics of Fluid-Phase Equilibria; Pearson Education, 1998. [Google Scholar]
  57. Finney A. R.; Salvalaglio M. Bridging the gap between mesoscopic and molecular models of solid/liquid interfaces out-of-equilibrium. Chem. Eng. Res. Des. 2022, 180, 285–295. 10.1016/j.cherd.2022.02.018. [DOI] [Google Scholar]
  58. Radu M.; Kremer K. Enhanced crystal growth in binary Lennard-Jones mixtures. Phys. Rev. Lett. 2017, 118, 055702. 10.1103/PhysRevLett.118.055702. [DOI] [PubMed] [Google Scholar]
  59. Han D.; Karmakar T.; Bjelobrk Z.; Gong J.; Parrinello M. Solvent-mediated morphology selection of the active pharmaceutical ingredient isoniazid: Experimental and simulation studies. Chem. Eng. Sci. 2019, 204, 320. 10.1016/j.ces.2018.10.022. [DOI] [Google Scholar]
  60. Elts E.; Luxenburger F.; Briesen H. Influence of Monovalent Salts on α-Glycine Crystal Growth from Aqueous Solution: Molecular Dynamics Simulations at Constant Supersaturation Conditions. J. Phys. Chem. B 2021, 125, 11732–11741. 10.1021/acs.jpcb.1c07168. [DOI] [PubMed] [Google Scholar]
  61. Barducci A.; Bussi G.; Parrinello M. Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100, 020603. 10.1103/PhysRevLett.100.020603. [DOI] [PubMed] [Google Scholar]
  62. Bjelobrk Z.; Piaggi P. M.; Weber T.; Karmakar T.; Mazzotti M.; Parrinello M. Naphthalene crystal shape prediction from molecular dynamics simulations. CrystEngComm 2019, 21, 3280. 10.1039/C9CE00380K. [DOI] [Google Scholar]
  63. Bjelobrk Z.; Mendels D.; Karmakar T.; Parrinello M.; Mazzotti M. Solubility prediction of organic molecules with molecular dynamics simulations. Cryst. Growth Des. 2021, 21, 5198–5205. 10.1021/acs.cgd.1c00546. [DOI] [Google Scholar]
  64. Bjelobrk Z.; Rajagopalan A. K.; Mendels D.; Karmakar T.; Parrinello M.; Mazzotti M. Solubility of Organic Salts in Solvent–Antisolvent Mixtures: A Combined Experimental and Molecular Dynamics Simulations Approach. J. Chem. Theory Comput. 2022, 18, 4952–4959. 10.1021/acs.jctc.2c00304. [DOI] [PMC free article] [PubMed] [Google Scholar]
  65. Ansari N.; Karmakar T.; Parrinello M. Molecular mechanism of gas solubility in liquid: constant chemical potential molecular dynamics simulations. J. Chem. Theory Comput. 2020, 16, 5279–5286. 10.1021/acs.jctc.0c00450. [DOI] [PubMed] [Google Scholar]
  66. Zhang D.; Jing X.; Sholl D. S.; Sinnott S. B. Molecular simulation of capture of sulfur-containing gases by porous aromatic frameworks. J. Phys. Chem. C 2018, 122, 18456–18467. 10.1021/acs.jpcc.8b03767. [DOI] [Google Scholar]
  67. Loganathan N.; Bowers G. M.; Wakou B. F. N.; Kalinichev A. G.; Kirkpatrick R. J.; Yazaydin A. O. Understanding methane/carbon dioxide partitioning in clay nano-and meso-pores with constant reservoir composition molecular dynamics modeling. Phys. Chem. Chem. Phys. 2019, 21, 6917–6924. 10.1039/C9CP00851A. [DOI] [PubMed] [Google Scholar]
  68. Loganathan N.; Yazaydin A. O.; Bowers G. M.; Ngouana-Wakou B. F.; Kalinichev A. G.; Kirkpatrick R. J. Role of Cations in the Methane/Carbon Dioxide Partitioning in Nano-and Mesopores of Illite Using Constant Reservoir Composition Molecular Dynamics Simulation. J. Phys. Chem. C 2020, 124, 2490–2500. 10.1021/acs.jpcc.9b10051. [DOI] [Google Scholar]
  69. Han D.; Karmakar T.; Bjelobrk Z.; Gong J.; Parrinello M. Solvent-mediated morphology selection of the active pharmaceutical ingredient isoniazid: Experimental and simulation studies. Chem. Eng. Sci. 2019, 204, 320–328. 10.1016/j.ces.2018.10.022. [DOI] [Google Scholar]
  70. Karmakar T.; Piaggi P. M.; Perego C.; Parrinello M. A Cannibalistic Approach to Grand Canonical Crystal Growth. J. Chem. Theory Comput. 2018, 14, 2678–2683. 10.1021/acs.jctc.8b00191. [DOI] [PubMed] [Google Scholar]
  71. Gupta K. M.; Qiao Z.; Zhang K.; Jiang J. Seawater Pervaporation through Zeolitic Imidazolate Framework Membranes: Atomistic Simulation Study. ACS Appl. Mater. Interfaces 2016, 8, 13392–13399. 10.1021/acsami.6b01626. [DOI] [PubMed] [Google Scholar]
  72. Ozcan A.; Semino R.; Maurin G.; Yazaydin A. O. Modeling of gas transport through polymer/MOF interfaces: a microsecond-scale concentration gradient-driven molecular dynamics study. Chem. Mater. 2020, 32, 1288–1296. 10.1021/acs.chemmater.9b04907. [DOI] [PMC free article] [PubMed] [Google Scholar]
  73. Madero-Castro R. M.; Calero S.; Yazaydin A. O. The role of hydrogen bonding in the dehydration of bioalcohols in hydrophobic pervaporation membranes. J. Mol. Liq. 2021, 340, 117297. 10.1016/j.molliq.2021.117297. [DOI] [Google Scholar]
  74. Namsani S.; Ozcan A.; Yazaydın A. O. Direct Simulation of Ternary Mixture Separation in a ZIF-8 Membrane at Molecular Scale. Advanced Theory and Simulations 2019, 2, 1900120. 10.1002/adts.201900120. [DOI] [Google Scholar]
  75. Vinh-Thang H.; Kaliaguine S. Predictive Models for Mixed-Matrix Membrane Performance: A Review. Chem. Rev. 2013, 113, 4980–5028. 10.1021/cr3003888. [DOI] [PubMed] [Google Scholar]
  76. Joung I. S.; Cheatham T. E. III Molecular dynamics simulations of the dynamic and energetic properties of alkali and halide ions using water-model-specific ion parameters. J. Phys. Chem. B 2009, 113, 13279–13290. 10.1021/jp902584c. [DOI] [PMC free article] [PubMed] [Google Scholar]
  77. Mester Z.; Panagiotopoulos A. Z. Mean ionic activity coefficients in aqueous NaCl solutions from molecular dynamics simulations. J. Chem. Phys. 2015, 142, 044507. 10.1063/1.4906320. [DOI] [PubMed] [Google Scholar]
  78. Espinosa J.; Young J.; Jiang H.; Gupta D.; Vega C.; Sanz E.; Debenedetti P. G.; Panagiotopoulos A. Z. On the calculation of solubilities via direct coexistence simulations: Investigation of NaCl aqueous solutions and Lennard-Jones binary mixtures. J. Chem. Phys. 2016, 145, 154111. 10.1063/1.4964725. [DOI] [PubMed] [Google Scholar]
  79. Piaggi P. M.; Parrinello M. Calculation of phase diagrams in the multithermal-multibaric ensemble. J. Chem. Phys. 2019, 150, 244119. 10.1063/1.5102104. [DOI] [PubMed] [Google Scholar]
  80. Di Pasquale N.; Finney A.; Elliott J.; Carbone P.; Salvalaglio M. Constant chemical potential-quantum mechanical-molecular dynamics simulations of the graphene-electrolyte double layer. J. Chem. Phys. 2023, 158, 134714. 10.1063/5.0138267. [DOI] [PubMed] [Google Scholar]
  81. Bonomi M.; Bussi G.; Camilloni C.; Tribello G. A.; Banáš P.; Barducci A.; Bernetti M.; Bolhuis P. G.; Bottaro S.; Branduardi D.; Capelli R.; Carloni P.; Ceriotti M.; Cesari A.; Chen H.; Chen W.; Colizzi F.; De S.; De La Pierre M.; Donadio D.; Drobot V.; Ensing B.; Ferguson A. L.; Filizola M.; Fraser J. S.; Fu H.; Gasparotto P.; Gervasio F. L.; Giberti F.; Gil-Ley A.; Giorgino T.; Heller G. T.; Hocky G. M.; Iannuzzi M.; Invernizzi M.; Jelfs K. E.; Jussupow A.; Kirilin E.; Laio A.; Limongelli V.; Lindorff-Larsen K.; Löhr T.; Marinelli F.; Martin-Samos L.; Masetti M.; Meyer R.; Michaelides A.; Molteni C.; Morishita T.; Nava M.; Paissoni C.; Papaleo E.; Parrinello M.; Pfaendtner J.; Piaggi P.; Piccini G.; Pietropaolo A.; Pietrucci F.; Pipolo S.; Provasi D.; Quigley D.; Raiteri P.; Raniolo S.; Rydzewski J.; Salvalaglio M.; Sosso G. C.; Spiwok V.; Šponer J.; Swenson D. W. H.; Tiwary P.; Valsson O.; Vendruscolo M.; Voth G. A.; White A. Promoting transparency and reproducibility in enhanced molecular simulations. Nat. Methods 2019, 16, 670–673. 10.1038/s41592-019-0506-8. [DOI] [PubMed] [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

ar2c00811_si_001.pdf (130KB, pdf)

Data Availability Statement

The original CμMD method and its variants have been implemented in the open-source software PLUMED development versions. The CμMD code, input files, and related tutorials are available from the authors’ GitHub pages (https://github.com/mme-ucl/CmuMD,https://github.com/Tarakk/plumed-cumd), materials clouds archive (10.24435/materialscloud:2020.0015/v1 and 10.24435/materialscloud:2020.0013/v1,10.24435/materialscloud:k5-t2), and PLUMED-NEST81 portal (IDs). A tutorial on applying symmetric and asymmetric planar symmetry CμMD is available on the PLUMED website as a part of the 2022 PLUMED Masterclass (https://www.plumed.org/doc-master/user-doc/html/masterclass-22-8.html).

Articles from Accounts of Chemical Research are provided here courtesy of American Chemical Society

RESOURCES