Ideal conductor/dielectric model (ICDM): A generalized technique to correct for finite-size effects in molecular simulations of hindered ion transport.

Molecular simulations serve as indispensable tools for investigating the kinetics and elucidating the mechanism of hindered ion transport across nanoporous membranes. In particular, recent advancements in advanced sampling techniques have made it possible to access translocation timescales spanning several orders of magnitude. In our prior study [Shoemaker et al., J. Chem. Theory Comput. 18, 7142 (2022)], we identified significant finite size artifacts in simulations of pressure-driven hindered ion transport through nanoporous graphitic membranes. We introduced the ideal conductor model, which effectively corrects for such artifacts by assuming the feed to be an ideal conductor. In the present work, we introduce the ideal conductor dielectric model (Icdm), a generalization of our earlier model, which accounts for the dielectric properties of both the membrane and the filtrate. Using the Icdm model substantially enhances the agreement among corrected free energy profiles obtained from systems of varying sizes, with notable improvements observed in regions proximate to the pore exit. Moreover, the model has the capability to consider secondary ion passage events, including the transport of a co-ion subsequent to the traversal of a counter-ion, a feature that is absent in our original model. We also investigate the sensitivity of the new model to various implementation details. The Icdm model offers a universally applicable framework for addressing finite size artifacts in molecular simulations of ion transport. It stands as a significant advancement in our quest to use molecular simulations to comprehensively understand and manipulate ion transport processes through nanoporous membranes.

How to quantify and avoid finite size effects in computational studies of crystal nucleation: The case of homogeneous crystal nucleation.

Finite size artifacts arise in molecular simulations of nucleation when critical nuclei are too close to their periodic images. A rigorous determination of what constitutes too close is, however, a major challenge. Recently, we devised rigorous heuristics for detecting such artifacts based on our investigation of how system size impacts the rate of heterogeneous ice nucleation [S. Hussain and A. Haji-Akbari, J. Chem. Phys. 154, 014108 (2021)]. We identified the prevalence of critical nuclei spanning across the periodic boundary, and the thermodynamic and structural properties of the liquid occupying the inter-image region as indicators of finite size artifacts. Here, we further probe the performance of such heuristics by examining the dependence of homogeneous crystal nucleation rates in the Lennard-Jones (LJ) liquid on system size. The rates depend non-monotonically on system size and vary by almost six orders of magnitude for the range of system sizes considered here. We confirm that the prevalence of spanning critical nuclei is the primary indicator of finite size artifacts and almost fully explains the observed variations in rate. Proximity, or structuring of the inter-image liquid, however, is not as strong of an indicator due to the fragmented nature of crystalline nuclei. As a result, the dependence of rate on system size is subtle for the systems with a minuscule fraction of spanning critical nuclei. These observations indicate that our heuristics are universally applicable to different modes of nucleation (homogeneous and heterogeneous) in different systems even if they might be overly stringent for homogeneous nucleation, e.g., in the LJ system.

Role of Nanoscale Interfacial Proximity in Contact Freezing in Water.

Contact freezing is a mode of atmospheric ice nucleation in which a collision between a dry ice nucleating particle (INP) and a water droplet results in considerably faster heterogeneous nucleation. The molecular mechanism of such an enhancement is, however, still a mystery. While earlier studies had attributed it to collision-induced transient perturbations, recent experiments point to the pivotal role of nanoscale proximity of the INP and the free interface. By simulating the heterogeneous nucleation of ice within INP-supported nanofilms of two model water-like tetrahedral liquids, we demonstrate that such nanoscale proximity is sufficient for inducing rate increases commensurate with those observed in contact freezing experiments, but only if the free interface has a tendency to enhance homogeneous nucleation. Water is suspected of possessing this latter property, known as surface freezing propensity. Our findings therefore establish a connection between the surface freezing propensity and kinetic enhancement during contact nucleation. We also observe that faster nucleation proceeds through a mechanism markedly distinct from classical heterogeneous nucleation, involving the formation of hourglass-shaped crystalline nuclei that conceive at either interface and that have a lower free energy of formation due to the nanoscale proximity of the interfaces and the modulation of the free interfacial structure by the INP. In addition to providing valuable insights into the physics of contact nucleation, our findings can assist in improving the accuracy of heterogeneous nucleation rate measurements in experiments and in advancing our understanding of ice nucleation on nonuniform surfaces such as organic, polymeric, and biological materials.

How to quantify and avoid finite size effects in computational studies of crystal nucleation: The case of heterogeneous ice nucleation.

Computational studies of crystal nucleation can be impacted by finite size effects, primarily due to unphysical interactions between crystalline nuclei and their periodic images. It is, however, not always feasible to systematically investigate the sensitivity of nucleation kinetics and mechanism to system size due to large computational costs of nucleation studies. Here, we use jumpy forward flux sampling to accurately compute the rates of heterogeneous ice nucleation in the vicinity of square-shaped model structureless ice nucleating particles (INPs) of different sizes and identify three distinct regimes for the dependence of rate on the INP dimension, L. For small INPs, the rate is a strong function of L due to the artificial spanning of critical nuclei across the periodic boundary. Intermediate-sized INPs, however, give rise to the emergence of non-spanning "proximal" nuclei that are close enough to their periodic images to fully structure the intermediary liquid. While such proximity can facilitate nucleation, its effect is offset by the higher density of the intermediary liquid, leading to artificially small nucleation rates overall. The critical nuclei formed at large INPs are neither spanning nor proximal. Yet, the rate is a weak function of L, with its logarithm scaling linearly with 1/L. The key heuristic emerging from these observations is that finite size effects will be minimal if critical nuclei are neither spanning nor proximal and if the intermediary liquid has a region that is structurally indistinguishable from the supercooled liquid under the same conditions.

Single crystal texture by directed molecular self-assembly along dual axes.

Creating well-defined single-crystal textures in materials requires the biaxial alignment of all grains into desired orientations, which is challenging to achieve in soft materials. Here we report the formation of single crystals with rigorously controlled texture over macroscopic areas (>1 cm2) in a soft mesophase of a columnar discotic liquid crystal. We use two modes of directed self-assembly, physical confinement and magnetic fields, to achieve control of the orientations of the columnar axes and the hexagonal lattice along orthogonal directions. Field control of the lattice orientation emerges in a low-temperature phase of tilted discogens that breaks the field degeneracy around the columnar axis present in non-tilted states. Conversely, column orientation is controlled by physical confinement and the resulting imposition of homeotropic anchoring at bounding surfaces. These results extend our understanding of molecular organization in tilted systems and may enable the development of a range of new materials for distinct applications.

Studying rare events using forward-flux sampling: Recent breakthroughs and future outlook.

Rare events are processes that occur upon the emergence of unlikely fluctuations. Unlike what their name suggests, rare events are fairly ubiquitous in nature, as the occurrence of many structural transformations in biology and material sciences is predicated upon crossing large free energy barriers. Probing the kinetics and uncovering the molecular mechanisms of possible barrier crossings in a system is critical to predicting and controlling its structural and functional properties. Due to their activated nature, however, rare events are exceptionally difficult to study using conventional experimental and computational techniques. In recent decades, a wide variety of specialized computational techniques-known as advanced sampling techniques-have been developed to systematically capture improbable fluctuations relevant to rare events. In this perspective, we focus on a technique called forward flux sampling [Allen et al., J. Chem. Phys. 124, 024102 (2006)] and overview its recent methodological variants and extensions. We also provide a detailed overview of its application to study a wide variety of rare events and map out potential avenues for further explorations.

Computational investigation of surface freezing in a molecular model of water.

Water freezes in a wide variety of low-temperature environments, from meteors and atmospheric clouds to soil and biological cells. In nature, ice usually nucleates at or near interfaces, because homogenous nucleation in the bulk can only be observed at deep supercoolings. Although the effect of proximal surfaces on freezing has been extensively studied, major gaps in understanding remain regarding freezing near vapor-liquid interfaces, with earlier experimental studies being mostly inconclusive. The question of how a vapor-liquid interface affects freezing in its vicinity is therefore still a major open question in ice physics. Here, we address this question computationally by using the forward-flux sampling algorithm to compute the nucleation rate in a freestanding nanofilm of supercooled water. We use the TIP4P/ice force field, one of the best existing molecular models of water, and observe that the nucleation rate in the film increases by seven orders of magnitude with respect to bulk at the same temperature. By analyzing the nucleation pathway, we conclude that freezing in the film initiates not at the surface, but within an interior region where the formation of double-diamond cages (DDCs) is favored in comparison with the bulk. This, in turn, facilitates freezing by favoring the formation of nuclei rich in cubic ice, which, as demonstrated by us earlier, are more likely to grow and overcome the nucleation barrier. The films considered here are ultrathin because their interior regions are not truly bulk-like, due to their subtle structural differences with the bulk.

Effect of material flexibility on the thermodynamics and kinetics of hydrophobically induced evaporation of water.

The evaporation of water induced by confinement between hydrophobic surfaces has received much attention due to its suggested functional role in numerous biophysical phenomena and its importance as a general mechanism of hydrophobic self-assembly. Although much progress has been made in understanding the basic physics of hydrophobically induced evaporation, a comprehensive understanding of the substrate material features (e.g., geometry, chemistry, and mechanical properties) that promote or inhibit such transitions remains lacking. In particular, comparatively little research has explored the relationship between water's phase behavior in hydrophobic confinement and the mechanical properties of the confining material. Here, we report the results of extensive molecular simulations characterizing the rates, free energy barriers, and mechanism of water evaporation when confined between model hydrophobic materials with tunable flexibility. A single-order-of-magnitude reduction in the material's modulus results in up to a nine-orders-of-magnitude increase in the evaporation rate, with the corresponding characteristic time decreasing from tens of seconds to tens of nanoseconds. Such a modulus reduction results in a 24-orders-of-magnitude decrease in the reverse rate of condensation, with time scales increasing from nanoseconds to tens of millions of years. Free energy calculations provide the barriers to evaporation and confirm our previous theoretical predictions that making the material more flexible stabilizes the confined vapor with respect to liquid. The mechanism of evaporation involves surface bubbles growing/coalescing to form a subcritical gap-spanning tube, which then must grow to cross the barrier.

Thermodynamically driven assemblies and liquid-liquid phase separations in biology.

The sustenance of life depends on the high degree of organization that prevails through different levels of living organisms, from subcellular structures such as biomolecular complexes and organelles to tissues and organs. The physical origin of such organization is not fully understood, and even though it is clear that cells and organisms cannot maintain their integrity without consuming energy, there is growing evidence that individual assembly processes can be thermodynamically driven and occur spontaneously due to changes in thermodynamic variables such as intermolecular interactions and concentration. Understanding the phase separation in vivo requires a multidisciplinary approach, integrating the theory and physics of phase separation with experimental and computational techniques. This paper aims at providing a brief overview of the physics of phase separation and its biological implications, with a particular focus on the assembly of membraneless organelles. We discuss the underlying physical principles of phase separation from its thermodynamics to its kinetics. We also overview the wide range of methods utilized for experimental verification and characterization of phase separation of membraneless organelles, as well as the utility of molecular simulations rooted in thermodynamics and statistical physics in understanding the governing principles of thermodynamically driven biological self-assembly processes.

Forward-flux sampling with jumpy order parameters.

Forward-flux sampling (FFS) is a path sampling technique that has gained increased popularity in recent years and has been used to compute rates of rare event phenomena such as crystallization, condensation, hydrophobic evaporation, DNA hybridization, and protein folding. The popularity of FFS is not only due to its ease of implementation but also because it is not very sensitive to the particular choice of an order parameter. The order parameter utilized in conventional FFS, however, still needs to satisfy a stringent smoothness criterion in order to assure sequential crossing of FFS milestones. This condition is usually violated for order parameters utilized for describing aggregation phenomena such as crystallization. Here, we present a generalized FFS algorithm for which this smoothness criterion is no longer necessary and apply it to compute homogeneous crystal nucleation rates in several systems. Our numerical tests reveal that conventional FFS can sometimes underestimate the nucleation rate by several orders of magnitude.

Forward flux sampling calculation of homogeneous nucleation rates from aqueous NaCl solutions.

We used molecular dynamics simulations and the path sampling technique known as forward flux sampling to study homogeneous nucleation of NaCl crystals from supersaturated aqueous solutions at 298 K and 1 bar. Nucleation rates were obtained for a range of salt concentrations for the Joung-Cheatham NaCl force field combined with the Extended Simple Point Charge (SPC/E) water model. The calculated nucleation rates are significantly lower than the available experimental measurements. The estimates for the nucleation rates in this work do not rely on classical nucleation theory, but the pathways observed in the simulations suggest that the nucleation process is better described by classical nucleation theory than an alternative interpretation based on Ostwald's step rule, in contrast to some prior simulations of related models. In addition to the size of NaCl nucleus, we find that the crystallinity of a nascent cluster plays an important role in the nucleation process. Nuclei with high crystallinity were found to have higher growth probability and longer lifetimes, possibly because they are less exposed to hydration water.

Direct calculation of ice homogeneous nucleation rate for a molecular model of water.

Ice formation is ubiquitous in nature, with important consequences in a variety of environments, including biological cells, soil, aircraft, transportation infrastructure, and atmospheric clouds. However, its intrinsic kinetics and microscopic mechanism are difficult to discern with current experiments. Molecular simulations of ice nucleation are also challenging, and direct rate calculations have only been performed for coarse-grained models of water. For molecular models, only indirect estimates have been obtained, e.g., by assuming the validity of classical nucleation theory. We use a path sampling approach to perform, to our knowledge, the first direct rate calculation of homogeneous nucleation of ice in a molecular model of water. We use TIP4P/Ice, the most accurate among existing molecular models for studying ice polymorphs. By using a novel topological approach to distinguish different polymorphs, we are able to identify a freezing mechanism that involves a competition between cubic and hexagonal ice in the early stages of nucleation. In this competition, the cubic polymorph takes over because the addition of new topological structural motifs consistent with cubic ice leads to the formation of more compact crystallites. This is not true for topological hexagonal motifs, which give rise to elongated crystallites that are not able to grow. This leads to transition states that are rich in cubic ice, and not the thermodynamically stable hexagonal polymorph. This mechanism provides a molecular explanation for the earlier experimental and computational observations of the preference for cubic ice in the literature.

Perspective: Surface freezing in water: A nexus of experiments and simulations.

Surface freezing is a phenomenon in which crystallization is enhanced at a vapor-liquid interface. In some systems, such as n-alkanes, this enhancement is dramatic and results in the formation of a crystalline layer at the free interface even at temperatures slightly above the equilibrium bulk freezing temperature. There are, however, systems in which the enhancement is purely kinetic and only involves faster nucleation at or near the interface. The first, thermodynamic, type of surface freezing is easier to confirm in experiments, requiring only the verification of the existence of crystalline order at the interface. The second, kinetic, type of surface freezing is far more difficult to prove experimentally. One material that is suspected of undergoing the second type of surface freezing is liquid water. Despite strong indications that the freezing of liquid water is kinetically enhanced at vapor-liquid interfaces, the findings are far from conclusive, and the topic remains controversial. In this perspective, we present a simple thermodynamic framework to understand conceptually and distinguish these two types of surface freezing. We then briefly survey fifteen years of experimental and computational work aimed at elucidating the surface freezing conundrum in water.

Computational investigation of structure, dynamics and nucleation kinetics of a family of modified Stillinger-Weber model fluids in bulk and free-standing thin films.

In recent years, computer simulations have found increasingly widespread use as powerful tools for studying phase transitions in wide variety of systems. In the particular and very important case of aqueous systems, the commonly used force-fields tend to offer quite different predictions with respect to a wide range of thermodynamic and kinetic properties, including the ease of ice nucleation, the propensity to freeze at a vapor-liquid interface, and the existence of a liquid-liquid phase transition. It is thus of fundamental and practical interest to understand how different features of a given water model affect its thermodynamic and kinetic properties. In this work, we use the forward-flux sampling technique to study the crystallization kinetics of a family of modified Stillinger-Weber (SW) potentials with energy (ε) and length (σ) scales taken from the monoatomic water (mW) model, but with different tetrahedrality parameters (λ). By increasing λ from 21 to 24, we observe the nucleation rate increases by 48 orders of magnitude at a supercooling of ζ = T/Tm = 0.845. Using classical nucleation theory, we are able to demonstrate that this change can largely be accounted for by the increase in |Δµ|, the thermodynamic driving force. We also perform rate calculations in freestanding thin films of the supercooled liquid, and observe a crossover from surface-enhanced crystallization at λ = 21 to bulk-dominated crystallization for λ ≥ 22.

Disordered, quasicrystalline and crystalline phases of densely packed tetrahedra.

All hard, convex shapes are conjectured by Ulam to pack more densely than spheres, which have a maximum packing fraction of phi = pi/ radical18 approximately 0.7405. Simple lattice packings of many shapes easily surpass this packing fraction. For regular tetrahedra, this conjecture was shown to be true only very recently; an ordered arrangement was obtained via geometric construction with phi = 0.7786 (ref. 4), which was subsequently compressed numerically to phi = 0.7820 (ref. 5), while compressing with different initial conditions led to phi = 0.8230 (ref. 6). Here we show that tetrahedra pack even more densely, and in a completely unexpected way. Following a conceptually different approach, using thermodynamic computer simulations that allow the system to evolve naturally towards high-density states, we observe that a fluid of hard tetrahedra undergoes a first-order phase transition to a dodecagonal quasicrystal, which can be compressed to a packing fraction of phi = 0.8324. By compressing a crystalline approximant of the quasicrystal, the highest packing fraction we obtain is phi = 0.8503. If quasicrystal formation is suppressed, the system remains disordered, jams and compresses to phi = 0.7858. Jamming and crystallization are both preceded by an entropy-driven transition from a simple fluid of independent tetrahedra to a complex fluid characterized by tetrahedra arranged in densely packed local motifs of pentagonal dipyramids that form a percolating network at the transition. The quasicrystal that we report represents the first example of a quasicrystal formed from hard or non-spherical particles. Our results demonstrate that particle shape and entropy can produce highly complex, ordered structures.

Thermodynamic and kinetic anisotropies in octane thin films.

Confinement breaks the translational symmetry of materials, making all thermodynamic and kinetic quantities functions of position. Such symmetry breaking can be used to obtain configurations that are not otherwise accessible in the bulk. Here, we use computer simulations to explore the effect of substrate-liquid interactions on thermodynamic and kinetic anisotropies induced by a solid substrate. We consider n-octane nano-films that are in contact with substrates with varying degrees of attraction, parameterized by an interaction parameter ϵS. Complete freezing of octane nano-films is observed at low temperatures, irrespective of ϵS, while at intermediate temperatures, a frozen monolayer emerges at solid-liquid and vapor-liquid interfaces. By carefully inspecting the profiles of translational and orientational relaxation times, we confirm that the translational and orientational degrees of freedom are decoupled at these frozen monolayers. At sufficiently high temperatures, however, free interfaces and solid-liquid interfaces close to loose (low-ϵS) substrates undergo "pre-freezing," characterized by mild peaks in several thermodynamic quantities. Two distinct dynamic regimes are observed at solid-liquid interfaces. The dynamics is accelerated in the vicinity of loose substrates, while sticky (high-ϵS) substrates decelerate dynamics, sometimes by as much as two orders of magnitude. These two distinct dynamical regimes have been previously reported by Haji-Akbari and Debenedetti [J. Chem. Phys. 141, 024506 (2014)] for a model atomic glass-forming liquid. We also confirm the existence of two correlations-proposed in the above-mentioned work-in solid-liquid subsurface regions of octane thin films, i.e., a correlation between atomic density and normal stress, and between atomic translational relaxation time and lateral stress. Finally, we inspect the ability of different regions of an octane film to explore the potential energy landscape by performing inherent structure calculations, and observe no noticeable difference between the free surface and the bulk in efficiently exploring the potential energy landscape. This is unlike the films of model atomic glass formers that tend to sample their respective landscape more efficiently at free surfaces. We discuss the implications of this finding to the ability of octane-and other n-alkanes-to form ultrastable glasses.

Suppression of sub-surface freezing in free-standing thin films of a coarse-grained model of water.

Freezing in the vicinity of water-vapor interfaces is of considerable interest to a wide range of disciplines, most notably the atmospheric sciences. In this work, we use molecular dynamics and two advanced sampling techniques, forward flux sampling and umbrella sampling, to study homogeneous nucleation of ice in free-standing thin films of supercooled water. We use a coarse-grained monoatomic model of water, known as mW, and we find that in this model a vapor-liquid interface suppresses crystallization in its vicinity. This suppression occurs in the vicinity of flat interfaces where no net Laplace pressure in induced. Our free energy calculations reveal that the pre-critical crystalline nuclei that emerge near the interface are thermodynamically less stable than those that emerge in the bulk. We investigate the origin of this instability by computing the average asphericity of nuclei that form in different regions of the film, and observe that average asphericity increases closer to the interface, which is consistent with an increase in the free energy due to increased surface-to-volume ratios.

The effect of substrate on thermodynamic and kinetic anisotropies in atomic thin films.

Glasses have a wide range of technological applications. The recent discovery of ultrastable glasses that are obtained by depositing the vapor of a glass-forming liquid onto the surface of a cold substrate has sparked renewed interest in the effects of confinements on physicochemical properties of liquids and glasses. Here, we use molecular dynamics simulations to study the effect of substrate on thin films of a model glass-forming liquid, the Kob-Andersen binary Lennard-Jones system, and compute profiles of several thermodynamic and kinetic properties across the film. We observe that the substrate can induce large oscillations in profiles of thermodynamic properties such as density, composition, and stress, and we establish a correlation between the oscillations in total density and the oscillations in normal stress. We also demonstrate that the kinetic properties of an atomic film can be readily tuned by changing the strength of interactions between the substrate and the liquid. Most notably, we show that a weakly attractive substrate can induce the emergence of a highly mobile region in its vicinity. In this highly mobile region, structural relaxation is several times faster than in the bulk, and the exploration of the potential energy landscape is also more efficient. In the subsurface region near a strongly attractive substrate, however, the dynamics is decelerated and the sampling of the potential energy landscape becomes less efficient than the bulk. We explain these two distinct behaviors by establishing a correlation between the oscillations in kinetic properties and the oscillations in lateral stress. Our findings offer interesting opportunities for designing better substrates for the vapor deposition process or developing alternative procedures for situations where vapor deposition is not feasible.

Estimating Position-Dependent and Anisotropic Diffusivity Tensors from Molecular Dynamics Trajectories: Existing Methods and Future Outlook.

Confinement can substantially alter the physicochemical properties of materials by breaking translational isotropy and rendering all physical properties position-dependent. Molecular dynamics (MD) simulations have proven instrumental in characterizing such spatial heterogeneities and probing the impact of confinement on materials' properties. For static properties, this is a straightforward task and can be achieved via simple spatial binning. Such an approach, however, cannot be readily applied to transport coefficients due to lack of natural extensions of autocorrelations used for their calculation in the bulk. The prime example of this challenge is diffusivity, which, in the bulk, can be readily estimated from the particles' mobility statistics, which satisfy the Fokker-Planck equation. Under confinement, however, such statistics will follow the Smoluchowski equation, which lacks a closed-form analytical solution. This brief review explores the rich history of estimating profiles of the diffusivity tensor from MD simulations and discusses various approximate methods and algorithms developed for this purpose. Besides discussing heuristic extensions of bulk methods, we overview more rigorous algorithms, including kernel-based methods, Bayesian approaches, and operator discretization techniques. Additionally, we outline methods based on applying biasing potentials or imposing constraints on tracer particles. Finally, we discuss approaches that estimate diffusivity from mean first passage time or committor probability profiles, a conceptual framework originally developed in the context of collective variable spaces describing rare events in computational chemistry and biology. In summary, this paper offers a concise survey of diverse approaches for estimating diffusivity from MD trajectories, highlighting challenges and opportunities in this area.