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

Nanocrystal Assemblies: Current Advances and Open Problems

We explore the potential of nanocrystals (a term used equivalently to nanoparticles) as building blocks for nanomaterials, and the current advances and open challenges for fundamental science developments and applications. Nanocrystal assemblies are inherently multiscale, and the generation of revolutionary material properties requires a precise understanding of the relationship between structure and function, the former being determined by classical effects and the latter often by quantum effects. With an emphasis on theory and computation, we discuss challenges that hamper current assembly strategies and to what extent nanocrystal assemblies represent thermodynamic equilibrium or kinetically trapped metastable states. We also examine dynamic effects and optimization of assembly protocols. Finally, we discuss promising material functions and examples of their realization with nanocrystal assemblies.

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.

ANNaMo: Coarse-grained modeling for folding and assembly of RNA and DNA systems

The folding of RNA and DNA strands plays crucial roles in biological systems and bionanotechnology. However, studying these processes with high-resolution numerical models is beyond current computational capabilities due to the timescales and system sizes involved. In this article, we present a new coarse-grained model for investigating the folding dynamics of nucleic acids. Our model represents three nucleotides with a patchy particle and is parameterized using well-established nearest-neighbor models. Thanks to the reduction of degrees of freedom and to a bond-swapping mechanism, our model allows for simulations at timescales and length scales that are currently inaccessible to more detailed models. To validate the performance of our model, we conducted extensive simulations of various systems: We examined the thermodynamics of DNA hairpins, capturing their stability and structural transitions, the folding of an MMTV pseudoknot, which is a complex RNA structure involved in viral replication, and also explored the folding of an RNA tile containing a k-type pseudoknot. Finally, we evaluated the performance of the new model in reproducing the melting temperatures of oligomers and the dependence on the toehold length of the displacement rate in toehold-mediated displacement processes, a key reaction used in molecular computing. All in all, the successful reproduction of experimental data and favorable comparisons with existing coarse-grained models validate the effectiveness of the new model.

Effect of Pressure on the Conformational Landscape of Human γD-Crystallin from Replica Exchange Molecular Dynamics Simulations

Human γD-crystallin belongs to a crucial family of proteins known as crystallins located in the fiber cells of the human lens. Since crystallins do not undergo any turnover after birth, they need to possess remarkable thermodynamic stability. However, their sporadic misfolding and aggregation, triggered by environmental perturbations or genetic mutations, constitute the molecular basis of cataracts, which is the primary cause of blindness in the globe according to the World Health Organization. Here, we investigate the impact of high pressure on the conformational landscape of wild-type HγD-crystallin using replica exchange molecular dynamics simulations augmented with principal component analysis. We find pressure to have a modest impact on global measures of protein stability, such as root-mean-square displacement and radius of gyration. Upon projecting our trajectories along the first two principal components from principal component analysis, however, we observe the emergence of distinct free energy basins at high pressures. By screening local order parameters previously shown or hypothesized as markers of HγD-crystallin stability, we establish correlations between a tyrosine-tyrosine aromatic contact within the N-terminal domain and the protein's end-to-end distance with projections along the first and second principal components, respectively. Furthermore, we observe the simultaneous contraction of the hydrophobic core and its intrusion by water molecules. This exploration sheds light on the intricate responses of HγD-crystallin to elevated pressures, offering insights into potential mechanisms underlying its stability and susceptibility to environmental perturbations, crucial for understanding cataract formation.

Comparing brute force to transition path sampling for gas hydrate nucleation with a flat interface: comments on time reversal symmetry

Fluid to solid nucleation is often investigated with the rare event method transition path sampling (TPS). I claim that the inherent irreversibility of solid nucleation, even at stationary conditions, calls into question TPS's applicability for determining solid nucleation mechanisms, especially for pre-critical behavior. Even when applied to a phenomenon which displays time reversal asymmetry like solid nucleation, TPS is a good means of exploring phase space and giving trends in post-critical structure, and its ability to facilitate nucleation rate and free energy calculations remains outstanding. Forward-only splitting and ratcheting methods such as forward flux sampling are more attractive for understanding nucleation mechanisms as they do not require time reversal symmetry, but at low driving forces may suffer from the same limitations as brute force: they may never make it to the first ratchet. Here I briefly summarize the TPS method and gas hydrate nucleation simulation literature, focusing on topics within both to facilitate a comparison of brute force hydrate nucleation to transition path sampling of hydrate nucleation. Perhaps anecdotally, the brute force technique results in more crystalline trajectories despite having higher driving forces than TPS. I maintain this difference is because of the inherent irreversibility of hydrate nucleation, meaning its pre-critical behavior cannot accurately be determined by the melting trajectories that comprise approximately half of the configurations in TPS's path ensemble.

Divergence among Local Structure, Dynamics, and Nucleation Outcome in Heterogeneous Nucleation of Close-Packed Crystals

Heterogeneous crystal nucleation is the dominant mechanism of crystallization in most systems, yet its underlying physics remains an enigma. While emergent interfacial crystalline order precedes heterogeneous nucleation, its importance in the nucleation mechanism is unclear. Here, we use path sampling simulations of two model systems to demonstrate that crystalline order in its traditional sense is not predictive of the outcome of the heterogeneous nucleation of close-packed crystals. Consequently, structure-based collective variables (CVs) that reliably describe homogeneous nucleation can be poor descriptors of heterogeneous nucleation. This divergence between structure and nucleation outcome is accompanied by an intriguing dynamical anomaly, wherein low-coordinated crystalline particles outpace their liquid-like counterparts. We use committor analysis, high-throughput screening, and machine learning to devise CV optimization strategies and present suitable structural heuristics within the metastable fluid for CV prescreening. Employing such optimized CVs is pivotal for properly characterizing the mechanism of heterogeneous nucleation in metallic and colloidal systems.

Correlations in Charged Multipore Systems: Implications for Enhancing Selectivity and Permeability in Nanoporous Membranes

Nanoporous membranes have emerged as powerful tools for diverse applications, including gas separation and water desalination. Achieving high permeability for desired molecules alongside exceptional rejection of other species presents a significant design challenge. One potential strategy involves optimizing the chemistry and geometry of isolated nanopores to enhance permeability and selectivity while maximizing their density within a membrane. However, the impact of the pore proximity on membrane performance remains an open question. Through path sampling simulations of model graphitic membranes with multiple subnanometer pores, we reveal that nanoscale proximity between pores detrimentally affects water permeability and salt rejection. Specifically, counterion transport is decelerated, while co-ion transport is accelerated, due to direct interactions among water molecules, salt ions, and the dipoles within neighboring pores. Notably, the observed ionic transport time scales significantly deviate from established theories such as the access resistance model but are well explained using the simple phenomenological model that we develop in this work. We use this model to prescreen and optimize pore arrangements that elicit minimal correlations at a target pore density. These findings deepen our understanding of multipore systems, informing the rational design of nanoporous membranes for enhanced separation processes such as water desalination. They also shed light on the physiology of biological cells that employ ion channel proteins to modulate ion transport and reversal potentials.

Kinetics of Oriented Attachment of Mica Crystals

Oriented attachment (OA), that is, the coalescence of crystals through attachment on coaligned crystal faces, is a nonclassical crystal growth process. Before attachment, a mesocrystal consisting of coaligned parallel crystals but with liquid separating them was observed. Fundamental questions such as why OA is kinetically favored and whether a mesocrystal stage is a prerequisite for OA are raised. Through combining brute-force molecular dynamics simulations and path samplings based on extensive umbrella simulations, we address these questions with a case study on the OA of a mica nanocrystal onto a mica crystal substrate in water. Brute-force simulations show that if two mica crystals are attached but largely misaligned, coalignment hardly appears. Thus, if OA is possible, then coalignment must appear before the attachment between crystals. Electrophoresis of the nanocrystal toward the substrate surface is spontaneous, but mesocrystal formation is occasional, also shown by brute-force simulations. Free energies along different pathways show that OA is spontaneous and kinetically favored over non-OA, and a mesocrystal formation is just a bifurcation in the pathway. OA is through a pathway in which the nanocrystal is tilted with respect to the substrate. Part of the nanocrystal is attached to the substrate first, and then, OA is gradually completed. Once a mesocrystal is occasionally formed, then a jump event is needed for the nanocrystal to get back to the OA pathway. The sampling technique here can hopefully guide the design of nanostructured materials facilitated by OA.

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.

Predicting Gaseous Solute Diffusion in Viscous Multivalent Ionic Liquid Solvents

Calculating solute diffusion in dense, viscous solvents can be particularly challenging in molecular dynamics simulations due to the long time scales involved. Here, a new scaling approach is developed for predicting solute diffusion based on analyses of CO2 and SO2 diffusion in two different multivalent ionic liquid solvents. Various scaling approaches are initially evaluated, including single and separate thermostats for the solute and solvent, as well as the application of the Arrhenius relationship and the Speedy-Angell power law. A very strong logarithmic correlation is established between the solvent-accessible surface area and solute diffusion. This relationship, reflecting Danckwerts' surface renewal theory and the Vrentas-Duda free volume model, presents a valuable method for estimating diffusion behavior from short simulation trajectories at elevated temperatures. The approach may be beneficial for enhancing predictive modeling in similar challenging systems and should be more broadly evaluated.

Does supercooled water retain its universal nucleation behavior under shear at high pressure?

Understanding the nucleation of homogeneous flow systems at high pressures is vital in protein crystallization and cryopreservation, where high pressure prevents the freezing of biological samples. This study examines the behavior of ice nucleation under shear at various pressures and explores the universal nucleation behavior of the sheared systems applied to supercooled water at higher pressures. In this study, the nucleation rates for the TIP4P/Ice model via a seeding method based on extended classical nucleation theory (CNT) are computed at pressures of 1, 100, 500, 700, and 1000 bar and a constant temperature of 240 K. Using extended CNT with explicitly embodying the shear rate, we analyzed the dependence of pressure on the transport and thermodynamic properties. In line with previous studies, we observed that Δμliq-ice and viscosity decrease while diffusivity increases with an increase in pressure. Furthermore, we showed that the dependence of the nucleation rate on shear at higher pressure is non-monotonic, with the maximum at optimal shear rates between 107 and 108 s-1. Our results demonstrate a non-monotonic pressure dependence of the optimal shear rates, which could originate from a violation of the Stokes-Einstein relation.

Interplay of multiple clusters and initial interface positioning for forward flux sampling simulations of crystal nucleation

Forward flux sampling (FFS) is a path sampling technique widely used in computer simulations of crystal nucleation from the melt. In such studies, the order parameter underpinning the progress of the FFS algorithm is often the size of the largest crystalline nucleus. In this work, we investigate the effects of two computational aspects of FFS simulations, using the prototypical Lennard-Jones liquid as our computational test bed. First, we quantify the impact of the positioning of the liquid basin and first interface in the space of the order parameter. In particular, we demonstrate that these choices are key to ensuring the consistency of the FFS results. Second, we focus on the frequently encountered scenario where the population of crystalline nuclei is such that there are multiple clusters of size comparable to the largest one. We demonstrate the contribution of clusters other than the largest cluster to the initial flux; however, we also show that they can be safely ignored for the purposes of converging a full FFS calculation. We also investigate the impact of different clusters merging, a process that appears to be facilitated by substantial spatial correlations-at least at the supercooling considered here. Importantly, all of our results have been obtained as a function of system size, thus contributing to the ongoing discussion on the impact of finite size effects on simulations of crystal nucleation. Overall, this work either provides or justifies several practical guidelines for performing FFS simulations that can also be applied to more complex and/or computationally expensive models.

Nonadiabatic Forward Flux Sampling for Excited-State Rare Events

We present a rare event sampling scheme applicable to coupled electronic excited states. In particular, we extend the forward flux sampling (FFS) method for rare event sampling to a nonadiabatic version (NAFFS) that uses the trajectory surface hopping (TSH) method for nonadiabatic dynamics. NAFFS is applied to two dynamically relevant excited-state models that feature an avoided crossing and a conical intersection with tunable parameters. We investigate how nonadiabatic couplings, temperature, and reaction barriers affect transition rate constants in regimes that cannot be otherwise obtained with plain, traditional TSH. The comparison with reference brute-force TSH simulations for limiting cases of rareness shows that NAFFS can be several orders of magnitude cheaper than conventional TSH and thus represents a conceptually novel tool to extend excited-state dynamics to time scales that are able to capture rare nonadiabatic events.

A study on prenucleation and heterogeneous nucleation in liquid Pb on solid Al using molecular dynamics simulations

In this paper, we investigate prenucleation and heterogeneous nucleation in the liquid Pb/solid Al system as an example of systems with large lattice misfit using molecular dynamics simulation. Solid Pb and Al have a large positive lattice misfit (f) of 18.2% along the densely packed [110] direction. This study reveals that prenucleation occurs at 600 K (an undercooling of 15 K), and a 2-dimensional (2D) ordered structure forms at the interface with a coincidence site lattice (CSL) between the first Pb and first Al layers. The CSL accommodates the major part of the f, and only a small residual lattice misfit (f_r) of 1.9% remains. The formation of the CSL transforms the original substrate into a considerably potent nucleant, where the first Pb layer becomes the new surface layer of the substrate. At an undercooling of about 22 K, nucleation proceeds by merging 2D ordered structure through structural templating: the second Pb layer is epitaxial to the CSL Pb layer, the third Pb layer largely accommodates the f_r, and the fourth Pb layer is a nearly perfect crystalline plane. Further analysis indicates that the interface with the CSL has a lower interfacial energy than with a cube-to-cube orientation relationship. For the first time, we established that the CSL was an effective mechanism to accommodate the f for systems with a large positive misfits. Heterogeneous nucleation is governed not by a single mechanism (misfit dislocations in Turnbull's model), but instead by various mechanisms depending on f. This study sheds new light on the atomistic mechanism of heterogeneous nucleation.

Ideal Conductor Model: An Analytical Finite-Size Correction for Nonequilibrium Molecular Dynamics Simulations of Ion Transport through Nanoporous Membranes

Modulating ion transport through nanoporous membranes is critical to many important chemical and biological separation processes. The corresponding transport timescales, however, are often too long to capture accurately using conventional molecular dynamics (MD). Recently, path sampling techniques, such as forward-flux sampling (FFS), have emerged as attractive alternatives for efficiently and accurately estimating arbitrarily long ionic passage times. Here, we use non-equilibrium MD and FFS to explore how the kinetics and mechanisms of pressure-driven chloride transport through a nanoporous graphitic membrane are affected by its lateral dimensions. We not only find ionic passage times and free energy barriers to decrease dramatically upon increasing the membrane surface area but also observe an abrupt and discontinuous change in the locus of the transition state. These strong finite size effects arise due to the cumulative effect of the periodic images of the leading ion entering the pore on the distribution of the induced excess charge at the membrane surface in the feed. By assuming that the feed is an ideal conductor, we analytically derive a finite size correction term that can be computed from the information obtained from a single simulation and successfully use it to obtain corrected free energy profiles with no dependence on the system size. We then estimate ionic passage times in the thermodynamic limit by assuming an Eyring-type dependence of rates on barriers with a size-independent prefactor. This approach constitutes a universal framework for removing finite size artifacts in molecular simulations of ion transport through nanoporous membranes and biological channel proteins.

Extending and validating bubble nucleation rate predictions in a Lennard-Jones fluid with enhanced sampling methods and transition state theory

We calculate bubble nucleation rates in a Lennard-Jones fluid through explicit molecular dynamics simulations. Our approach—based on a recent free energy method (dubbed reweighted Jarzynski sampling), transition state theory, and a simple recrossing correction—allows us to probe a fairly wide range of rates in several superheated and cavitation regimes in a consistent manner. Rate predictions from this approach bridge disparate independent literature studies on the same model system. As such, we find that rate predictions based on classical nucleation theory, direct brute force molecular dynamics simulations, and seeding are consistent with our approach and one another. Published rates derived from forward flux sampling simulations are, however, found to be outliers. This study serves two purposes: First, we validate the reliability of common modeling techniques and extrapolation approaches on a paradigmatic problem in materials science and chemical physics. Second, we further test our highly generic recipe for rate calculations, and establish its applicability to nucleation processes.

Identifying a Feasible Transition Pathway between Two Conformational States for a Protein

Proteins usually need to transit between different conformational states to fulfill their biological functions. In the mechanistic study of such transition processes by molecular dynamics simulations, identification of the minimum free energy path (MFEP) can substantially reduce the sampling space, thus enabling rigorous thermodynamic evaluation of the process. Conventionally, the MFEP is derived by iterative local optimization from an initial path, which is typically generated by simple brute force techniques like the targeted molecular dynamics (tMD). Therefore, the quality of the initial path determines the successfulness of MFEP estimation. In this work, we propose a method to improve derivation of the initial path. Through iterative relaxation-biasing simulations in a bidirectional manner, this method can construct a feasible transition pathway connecting two known states for a protein. Evaluation on small, fast-folding proteins against long equilibrium trajectories supports the good sampling efficiency of our method. When applied to larger proteins including the catalytic domain of human c-Src kinase as well as the converter domain of myosin VI, the paths generated by our method deviate significantly from those computed with the generic tMD approach. More importantly, free energy profiles and intermediate states obtained from our paths exhibit remarkable improvements over those from tMD paths with respect to both physical rationality and consistency with a priori knowledge.

Practical guide to replica exchange transition interface sampling and forward flux sampling

Path sampling approaches have become invaluable tools to explore the mechanisms and dynamics of the so-called rare events that are characterized by transitions between metastable states separated by sizable free energy barriers. Their practical application, in particular to ever more complex molecular systems, is, however, not entirely trivial. Focusing on replica exchange transition interface sampling (RETIS) and forward flux sampling (FFS), we discuss a range of analysis tools that can be used to assess the quality and convergence of such simulations, which is crucial to obtain reliable results. The basic ideas of a step-wise evaluation are exemplified for the study of nucleation in several systems with different complexities, providing a general guide for the critical assessment of RETIS and FFS simulations.

Modeling of biomolecular machines in non-equilibrium steady states

Numerical computations have become a pillar of all modern quantitative sciences. Any computation involves modeling-even if often this step is not made explicit-and any model has to neglect details while still being physically accurate. Equilibrium statistical mechanics guides both the development of models and numerical methods for dynamics obeying detailed balance. For systems driven away from thermal equilibrium, such a universal theoretical framework is missing. For a restricted class of driven systems governed by Markov dynamics and local detailed balance, stochastic thermodynamics has evolved to fill this gap and to provide fundamental constraints and guiding principles. The next step is to advance stochastic thermodynamics from simple model systems to complex systems with tens of thousands or even millions of degrees of freedom. Biomolecules operating in the presence of chemical gradients and mechanical forces are a prime example for this challenge. In this Perspective, we give an introduction to isothermal stochastic thermodynamics geared toward the systematic multiscale modeling of the conformational dynamics of biomolecular and synthetic machines, and we outline some of the open challenges.

Gaussian-Accelerated Molecular Dynamics with the Weighted Ensemble Method: A Hybrid Method Improves Thermodynamic and Kinetic Sampling

Gaussian-accelerated molecular dynamics (GaMD) is a well-established enhanced sampling method for molecular dynamics simulations that effectively samples the potential energy landscape of the system by adding a boost potential, which smoothens the surface and lowers the energy barriers between states. GaMD is unable to give time-dependent properties such as kinetics directly. On the other hand, the weighted ensemble (WE) method can efficiently sample transitions between states with its many weighted trajectories, which directly yield rates and pathways. However, convergence to equilibrium conditions remains a challenge for the WE method. Hence, we have developed a hybrid method that combines the two methods, wherein GaMD is first used to sample the potential energy landscape of the system and WE is subsequently used to further sample the potential energy landscape and kinetic properties of interest. We show that the hybrid method can sample both thermodynamic and kinetic properties more accurately and quickly compared to using either method alone.

Concept of the Energy-Dependent Temperature for Direct Langevin Dynamics Simulations of Rare Events in Systems with Arbitrary High Energy Barriers

We suggest an algorithm which allows single-stage direct Langevin dynamics simulations of transitions over arbitrarily high energy barriers. For this purpose, we propose a concept of the energy-dependent temperature (EDT): near the energy minima this temperature is high, but it tends toward room temperature for energies approaching the barrier value. In the resulting algorithm simulation time required for the computation of the escape rate over the barrier does not increase with barrier height. Switching times computed via our EDT algorithm agree very well with those obtained with the forward flux sampling (FFS). As the simulation time required by our method does not increase with the energy barrier, we achieve a very large speed-up compared even to our strongly optimized version of FFS (and all other multistage algorithms). In addition, our approach is free from the instability occurring in all multistage "climbing" methods where a product of a large number of transition probabilities between the interfaces must be computed.

Deciphering the molecular mechanism of water boiling at heterogeneous interfaces

Water boiling control evolution of natural geothermal systems is widely exploited in industrial processes due to the unique non-linear thermophysical behavior. Even though the properties of water both in the liquid and gas state have been extensively studied experimentally and by numerical simulations, there is still a fundamental knowledge gap in understanding the mechanism of the heterogeneous nucleate boiling controlling evaporation and condensation. In this study, the molecular mechanism of bubble nucleation at the hydrophilic and hydrophobic solid-water interface was determined by performing unbiased molecular dynamics simulations using the transition path sampling scheme. Analyzing the liquid to vapor transition path, the initiation of small void cavities (vapor bubbles nuclei) and their subsequent merging mechanism, leading to successively growing vacuum domains (vapor phase), has been elucidated. The molecular mechanism and the boiling nucleation sites' location are strongly dependent on the solid surface hydrophobicity and hydrophilicity. Then simulations reveal the impact of the surface functionality on the adsorbed thin water molecules film structuring and the location of high probability nucleation sites. Our findings provide molecular-scale insights into the computational aided design of new novel materials for more efficient heat removal and rationalizing the damage mechanisms.

Nonparametric Analysis of Nonequilibrium Simulations

We extend the nonparametric framework of reaction coordinate optimization to nonequilibrium ensembles of (short) trajectories. For example, we show how, starting from such an ensemble, one can obtain an equilibrium free-energy profile along the committor, which can be used to determine important properties of the dynamics exactly. A new adaptive sampling approach, the transition-state ensemble enrichment, is suggested, which samples the configuration space by "growing" committor segments toward each other starting from the boundary states. This framework is suggested as a general tool, alternative to the Markov state models, for a rigorous and accurate analysis of simulations of large biomolecular systems, as it has the following attractive properties. It is immune to the curse of dimensionality, does not require system-specific information, can approximate arbitrary reaction coordinates with high accuracy, and has sensitive and rigorous criteria to test optimality and convergence. The approaches are illustrated on a 50-dimensional model system and a realistic protein folding trajectory.

Multiscale aspects of molecular motions: From molecular vibrations, conformational changes of biomolecules to cellular dynamics

Molecular aspects of living systems are important because it is the most basic aspects of life as exemplified in biochemistry and structural biology. Since molecules move due to interactive forces between atoms, physics plays an important role to understand the dynamic phenomena of living systems. Here we review our multiscale approaches for computationally treating different levels of molecular motions: vibrational dynamics of molecules, conformational change of biomolecules, and cellular dynamics using statistical-mechanics-based models.

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.

Parasitic crystallization of colloidal electrolytes: growing a metastable crystal from the nucleus of a stable phase

Colloidal particles have been extensively used to comprehend the main principles governing liquid-crystal nucleation. Multiple mechanisms and frameworks have been proposed, through either experiments or computational approaches, to rationalise the ubiquitous formation of colloidal crystals. In this work, we elucidate the nucleation scenario behind the crystallization of oppositely charged colloids. By performing molecular dynamics simulations of colloidal electrolytes in combination with the Seeding technique, we evaluate the fundamental factors, such as the nucleation rate, free energy barrier, surface tension and kinetic pre-factor, that determine the liquid-to-solid transition of several crystalline polymorphs. Our results show that at a high packing fraction, there is a cross-over between the nucleation of the CsCl structure and that of a substitutionally disordered fcc phase, despite the CuAu crystal being the most stable phase. We demonstrate that the crucial factor in determining which phase nucleates the fastest is the free energy cost of the cluster formation rather than their kinetic ability to grow from the liquid. While at a low packing fraction, the stable phase, CsCl, is the one that nucleates and subsequently grows, we show how at moderate and high packing fractions, a disordered fcc phase subsequently grows regardless of the nature of the nucleating phase, termed parasitic crystallization. Taken together, our results provide a panoramic perspective of the complex nucleation scenario of oppositely charged colloids at moderate temperature and rationalise the different thermodynamic and kinetic aspects behind it.

Target Search of Active Agents Crossing High Energy Barriers

Target search by active agents in rugged energy landscapes has remained a challenge because standard enhanced sampling methods do not apply to irreversible dynamics. We overcome this nonequilibrium rare-event problem by developing an algorithm generalizing transition-path sampling to active Brownian dynamics. This method is exemplified and benchmarked for a paradigmatic two-dimensional potential with a high barrier. We find that even in such a simple landscape the structure and kinetics of the ensemble of transition paths changes drastically in the presence of activity. Indeed, active Brownian particles reach the target more frequently than passive Brownian particles, following longer and counterintuitive search patterns.

Homogeneous nucleation of sheared liquids: advances and insights from simulations and theory

One of the most ubiquitous and technologically important phenomena in nature is the nucleation of homogeneous flowing systems. The microscopic effects of shear on a nucleating system are still imperfectly understood, although in recent years a consistent picture has emerged. The opposing effects of shear can be split into two major contributions for simple atomic and molecular liquids: increase of the energetic cost of nucleation, and enhancement of the kinetics. In this perspective, we describe the latest computational and theoretical techniques which have been developed over the past two decades. We collate and unify the overarching influences of shear, temperature, and supersaturation on the process of homogeneous nucleation. Experimental techniques and capabilities are discussed, against the backdrop of results from simulations and theory. Although we primarily focus on simple systems, we also touch upon the sheared nucleation of more complex systems, including glasses and polymer melts. We speculate on the promising directions and possible advances that could come to fruition in the future.

Validation of DBFOLD: An efficient algorithm for computing folding pathways of complex proteins

Atomistic simulations can provide valuable, experimentally-verifiable insights into protein folding mechanisms, but existing ab initio simulation methods are restricted to only the smallest proteins due to severe computational speed limits. The folding of larger proteins has been studied using native-centric potential functions, but such models omit the potentially crucial role of non-native interactions. Here, we present an algorithm, entitled DBFOLD, which can predict folding pathways for a wide range of proteins while accounting for the effects of non-native contacts. In addition, DBFOLD can predict the relative rates of different transitions within a protein’s folding pathway. To accomplish this, rather than directly simulating folding, our method combines equilibrium Monte-Carlo simulations, which deploy enhanced sampling, with unfolding simulations at high temperatures. We show that under certain conditions, trajectories from these two types of simulations can be jointly analyzed to compute unknown folding rates from detailed balance. This requires inferring free energies from the equilibrium simulations, and extrapolating transition rates from the unfolding simulations to lower, physiologically-reasonable temperatures at which the native state is marginally stable. As a proof of principle, we show that our method can accurately predict folding pathways and Monte-Carlo rates for the well-characterized Streptococcal protein G. We then show that our method significantly reduces the amount of computation time required to compute the folding pathways of large, misfolding-prone proteins that lie beyond the reach of existing direct simulation. Our algorithm, which is available online, can generate detailed atomistic models of protein folding mechanisms while shedding light on the role of non-native intermediates which may crucially affect organismal fitness and are frequently implicated in disease.

Many proteins must adopt a specific structure in order to function. Computational simulations have been used to shed light on the mechanisms of protein folding, but unfortunately, realistic simulations can typically only be run for small proteins, due to severe limits in computational speed. Here, we present a method to solve this problem, whereby instead of directly simulating folding from an unfolded state, we run simulations that allow for computation of equilibrium folding free energies, alongside high temperature simulations to compute unfolding rates. From these quantities, folding rates can be computed using detailed balance. Importantly, our method can account for the effects of nonnative contacts which transiently form during folding and must be broken prior to adoption of the native state. Such contacts, which are often excluded from simple models of folding, may crucially affect real protein folding pathways and are often observed in folding intermediates implicated in disease.

Enantiomeric Transitions in the Chiral Cluster Au15 Studied by a Reaction Coordinate Deduced from Molecular Dynamics Simulations

A recently developed modified basin hopping (MBH) optimization algorithm, combined with an energy function calculated by the semiempirical density functional tight-binding (DFTB) theory, was applied to determine the lowest-energy structures of Au_n clusters with size n = 3-20. It was predicted from the DFTB/MBH optimization algorithm calculations that clusters Au₁₀, Au₁₅, and Au₁₈ exhibit chiral properties; i.e., each of these three clusters possesses the same energy value and associated with it are two nonsuperposable mirror-image clusters. In the potential energy landscape, there thus exist multidimensional barriers separating the two enantiomers, and this lowest-energy double-well morphology is surrounded by potential-energy minima of higher energies. In this paper, we have chosen to study the chiral cluster Au₁₅ by employing an isothermal Brownian-type molecular dynamics simulation to discern in greater detail its conformational transition from one enantiomer, say left, to its right counterpart. To facilitate our analysis of the simulation data, we transpose the multidimensional configurational space description to a lower dimensional collective variable (CV) space spanned by two geometry-relevant CVs. The thermally driven progression and mechanism of enantiomeric transitions between the left and right enantiomers will be our main focus, and the strategy is to dissect the time development of the CVs collected from different sets of independent simulation runs. From simulation data, we found that an understanding of the dynamics of enantiomeric transitions needs first to seek out the left and right enantiomers through a molecular modeling and visualizing program, then to ferret out and identify between the left and right enantiomers a symmetrical structure, and finally to define from the latter a reaction coordinate. We showed in this work that this single reaction coordinate is predictive in unraveling the left ⇌ right enantiomeric transition events, providing a specific inkling of the transition time span and its associated distribution which can be checked further for its reasonableness by the autocorrelation function and a vibrational analysis, all of which shed light on the mechanisms of transition.

Forward Flux Sampling of Polymer Desorption Paths from a Solid Surface into Dilute Solution

We compute desorption rates for isolated polymers adsorbed to a solid wall with a rare event sampling technique called multilevel splitting, also known as forward flux sampling. We interpret computed rates with theories based on the conjecture that the product tdesDRg2 of the desorption time tdes and diffusivity D divided by squared radius of gyration Rg scales with exp(h/Rg) where h is the equilibrium ratio of adsorbed surface concentration of polymer Γ to bulk concentration of polymer c. As the polymer-wall interaction energy is increased, the slope of lntdesDRg2 vs. NVMFkBT nearly approaches unity, as expected for strongly-adsorbing chains, where N is the degree of polymerization and VMF is the height-averaged monomer-wall interaction energy for a strongly adsorbed chain. However, we also find that this scaling law is only accurate when adsorption strength per monomer exceeds a threshold value on the order of 0.3-0.5 kBT for a freely jointed chain without or with excluded volume effects. Below the critical value, we observe that tdesDRg2 becomes nearly constant with N, so that tdes∝Nα, with α≈2. This suggests a crossover from "strong" detachment-controlled to a "weak" diffusion-controlled desorption rate as VMF/kBT drops below some threshold. These results may partially explain experimental data, that in some cases show "strong" exponential dependence of desorption time on chain length, while in others a "weak" power-law dependence is found. However, in the "strong" adsorption case, our results suggest much longer desorption times than those measured, while the reverse is true in the weak adsorption limit. We discuss possible reasons for these discrepancies.