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

Abstract

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

Figure 1

Scheme of the CμMD method. (a) Solute concentration profile c_i 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). Adapted with permission from ref (1) (copyright 2015, AIP).

Figure 2

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, Δg^E, for ions at the interface vs bulk concentration. 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).

Figure 3

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 (J_z) across a ZIF-8 membrane as a function of the concentration gradient Δc. (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. 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. Adapted with permission from ref (73) (copyright 2021, Elsevier).

Figure 4

(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, z_I^S and z_I^U, 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).

Figure 5

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 (s^O and s^H'), 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.