Comparison of kinetic Monte Carlo and molecular dynamics simulations of diffusion in a model glass former

Glassy dynamics in a liquid of anisotropic molecules: Bifurcation of relaxation spectrum

In experimental and theoretical studies of glass transition phenomena, one often finds a sharp crossover in dynamical properties at a temperature Tcr. A bifurcation of a relaxation spectrum is also observed at a temperature TB≈Tcr; both lie significantly above the glass transition temperature. In order to better understand these phenomena, we introduce a new model of glass-forming liquids, a binary mixture of prolate and oblate ellipsoids. This model system exhibits sharp thermodynamic and dynamic anomalies, such as the specific heat jump during heating and a sharp variation in the thermal expansion coefficient around a temperature identified as the glass transition temperature, Tg. The same temperature is obtained from the fit of the calculated relaxation times to the Vogel-Fulcher-Tammann (VFT) form. As the temperature is lowered, the calculated single peak rotational relaxation spectrum splits into two peaks at TB above the estimated Tg. Similar bifurcation is also observed in the distribution of short-to-intermediate time translational diffusion. Interrogation of the two peaks reveals a lower extent of dynamic heterogeneity in the population of the faster mode. We observe an unexpected appearance of a sharp peak in the product of rotational relaxation time τ2 and diffusion constant D at a temperature Tcr, close to TB, but above the glass transition temperature. Additionally, we coarse-grain the system into cubic boxes, each containing, on average, ∼62 particles, to study the average dynamical properties. Clear evidence of large-scale sudden changes in the diffusion coefficient and rotational correlation time signals first-order transitions between low and high-mobility domains.

A lattice kinetic Monte-Carlo method for simulating chromosomal dynamics and other (non-)equilibrium bio-assemblies

Biological assemblies in living cells such as chromosomes constitute large many-body systems that operate in a fluctuating, out-of-equilibrium environment. Since a brute-force simulation of that many degrees of freedom is currently computationally unfeasible, it is necessary to perform coarse-grained stochastic simulations. Here, we develop all tools necessary to write a lattice kinetic Monte-Carlo (LKMC) algorithm capable of performing such simulations. We discuss the validity and limits of this approach by testing the results of the simulation method in simple settings. Importantly, we illustrate how at large external forces Metropolis-Hastings kinetics violate the fluctuation-dissipation and steady-state fluctuation theorems and discuss better alternatives. Although this simulation framework is rather general, we demonstrate our approach using a DNA polymer with interacting SMC condensin loop-extruding enzymes. Specifically, we show that the scaling behavior of the loop-size distributions that we obtain in our LKMC simulations of this SMC-DNA system is consistent with that reported in other studies using Brownian dynamics simulations and analytic approaches. Moreover, we find that the irreversible dynamics of these enzymes under certain conditions result in frozen, sterically jammed polymer configurations, highlighting a potential pitfall of this approach.

Vertical drying of a suspension of sticks: Monte Carlo simulation for continuous two-dimensional problem

The vertical drying of a two-dimensional colloidal film containing zero-thickness sticks (lines) was studied by means of kinetic Monte Carlo (MC) simulations. The continuous two-dimensional problem for both the positions and orientations was considered. The initial state before drying was produced using a model of random sequential adsorption with isotropic orientations of the sticks. During the evaporation, an upper interface falls with a linear velocity in the vertical direction, and the sticks undergo translational and rotational Brownian motions. The MC simulations were run at different initial number concentrations (the numbers of sticks per unit area), p_{i}, and solvent evaporation rates, u. For completely dried films, the spatial distributions of the sticks, the order parameters, and the electrical conductivities of the films in both the horizontal, x, and vertical, y, directions were examined. Significant evaporation-driven self-assembly and stratification of the sticks in the vertical direction was observed. The extent of stratification increased with increasing values of u. The anisotropy of the electrical conductivity of the film can be finely regulated by changes in the values of p_{i} and u.

Pathways for diffusion in the potential energy landscape of the network glass former SiO2

We study the dynamical behaviour of a computer model for viscous silica, the archetypal strong glass former, and compare its diffusion mechanism with earlier studies of a fragile binary Lennard-Jones liquid. Three different methods of analysis are employed. First, the temperature and time scale dependence of the diffusion constant is analysed. Negative correlation of particle displacements influences transport properties in silica as well as in fragile liquids. We suggest that the difference between Arrhenius and super-Arrhenius diffusive behaviour results from competition between the correlation time scale and the caging time scale. Second, we analyse the dynamics using a geometrical definition of cage-breaking transitions that was proposed previously for fragile glass formers. We find that this definition accurately captures the bond rearrangement mechanisms that control transport in open network liquids, and reproduces the diffusion constants accurately at low temperatures. As the same method is applicable to both strong and fragile glass formers, we can compare correlation time scales in these two types of systems. We compare the time spent in chains of correlated cage breaks with the characteristic caging time and find that correlations in the fragile binary Lennard-Jones system persist for an order of magnitude longer than those in the strong silica system. We investigate the origin of the correlation behaviour by sampling the potential energy landscape for silica and comparing it with the binary Lennard-Jones model. We find no qualitative difference between the landscapes, but several metrics suggest that the landscape of the fragile liquid is rougher and more frustrated. Metabasins in silica are smaller than those in binary Lennard-Jones and contain fewer high-barrier processes. This difference probably leads to the observed separation of correlation and caging time scales.

Reaction sampling and reactivity prediction using the stochastic surface walking method

The new theoretical method demonstrates the ability of automated reaction sampling and activity prediction for complex organic reactions.

Temporal disconnectivity of the energy landscape in glassy systems

An alternative graphical representation of the potential energy landscape (PEL) has been developed and applied to a binary Lennard-Jones glassy system, providing insight into the unique topology of the system's potential energy hypersurface. With the help of this representation one is able to monitor the different explored basins of the PEL, as well as how--and mainly when--subsets of basins communicate with each other via transitions in such a way that details of the prior temporal history have been erased, i.e., local equilibration between the basins in each subset has been achieved. In this way, apart from detailed information about the structure of the PEL, the system's temporal evolution on the PEL is described. In order to gather all necessary information about the identities of two or more basins that are connected with each other, we consider two different approaches. The first one is based on consideration of the time needed for two basins to mutually equilibrate their populations according to the transition rate between them, in the absence of any effect induced by the rest of the landscape. The second approach is based on an analytical solution of the master equation that explicitly takes into account the entire explored landscape. It is shown that both approaches lead to the same result concerning the topology of the PEL and dynamical evolution on it. Moreover, a "temporal disconnectivity graph" is introduced to represent a lumped system stemming from the initial one. The lumped system is obtained via a specially designed algorithm [N. Lempesis, D. G. Tsalikis, G. C. Boulougouris, and D. N. Theodorou, J. Chem. Phys. 135, 204507 (2011)]. The temporal disconnectivity graph provides useful information about both the lumped and the initial systems, including the definition of "metabasins" as collections of basins that communicate with each other via transitions that are fast relative to the observation time. Finally, the two examined approaches are compared to an "on the fly" molecular dynamics-based algorithm [D. G. Tsalikis, N. Lempesis, G. C. Boulougouris, and D. N. Theodorou, J. Chem. Theory Comput. 6, 1307 (2010)].

The Activation-Relaxation Technique: ART Nouveau and Kinetic ART

The evolution of many systems is dominated by rare activated events that occur on timescale ranging from nanoseconds to the hour or more. For such systems, simulations must leave aside the full thermal description to focus specifically on mechanisms that generate a configurational change. We present here the activation relaxation technique (ART), an open-ended saddle point search algorithm, and a series of recent improvements to ART nouveau and kinetic ART, an ART-based on-the-fly off-lattice self-learning kinetic Monte Carlo method.

Kinetic activation-relaxation technique

We present a detailed description of the kinetic activation-relaxation technique (k-ART), an off-lattice, self-learning kinetic Monte Carlo (KMC) algorithm with on-the-fly event search. Combining a topological classification for local environments and event generation with ART nouveau, an efficient unbiased sampling method for finding transition states, k-ART can be applied to complex materials with atoms in off-lattice positions or with elastic deformations that cannot be handled with standard KMC approaches. In addition to presenting the various elements of the algorithm, we demonstrate the general character of k-ART by applying the algorithm to three challenging systems: self-defect annihilation in c-Si (crystalline silicon), self-interstitial diffusion in Fe, and structural relaxation in a-Si (amorphous silicon).

Kinetic Monte Carlo method for simulating reactions in solutions

We present an off-lattice kinetic Monte Carlo method, which is useful to simulate reactions in solutions. We derive the method from first-principles. We assume that diffusion leads to a Gaussian distribution for the position of the particles. This allows us to deal with the diffusion analytically, and we only need to simulate the reactive processes. The rate constants of these reactions can be computed before a simulation is started, and need not be computed on-the-fly as in other off-lattice kinetic Monte Carlo methods. We show how solvent molecules can be removed from the simulations, which minimizes the number of particles that have to be simulated explicitly. We present the relation with the customary macroscopic rate equations, and compare the results of these equations and our method on a variation of the Lotka model.

Phase space theory of evaporation in neon clusters: the role of quantum effects

Unimolecular evaporation of neon clusters containing between 14 and 148 atoms is theoretically investigated in the framework of phase space theory. Quantum effects are incorporated in the vibrational densities of states, which include both zero-point and anharmonic contributions, and in the possible tunneling through the centrifugal barrier. The evaporation rates, kinetic energy released, and product angular momentum are calculated as a function of excess energy or temperature in the parent cluster and compared to the classical results. Quantum fluctuations are found to generally increase both the kinetic energy released and the angular momentum of the product, but the effects on the rate constants depend nontrivially on the excess energy. These results are interpreted as due to the very few vibrational states available in the product cluster when described quantum mechanically. Because delocalization also leads to much narrower thermal energy distributions, the variations of evaporation observables as a function of canonical temperature appear much less marked than in the microcanonical ensemble. While quantum effects tend to smooth the caloric curve in the product cluster, the melting phase change clearly keeps a signature on these observables. The microcanonical temperature extracted from fitting the kinetic energy released distribution using an improved Arrhenius form further suggests a backbending in the quantum Ne(13) cluster that is absent in the classical system. Finally, in contrast to delocalization effects, quantum tunneling through the centrifugal barrier does not play any appreciable role on the evaporation kinetics of these rather heavy clusters.

Equilibrium density of states and thermodynamic properties of a model glass former

This paper presents an analysis of the thermodynamics of a model glass former. We have performed equilibrium sampling of a popular binary Lennard-Jones model, employing parallel tempering Monte Carlo to cover the crystalline, amorphous, and liquid regions of configuration space. Disconnectivity graphs are used to visualize the potential energy landscape in the vicinity of a crystalline geometry and in an amorphous region of configuration space. The crystalline global minimum is separated from the bulk of the minima by a large potential energy gap, leading to broken ergodicity in conventional simulations. Our sampling reveals crystalline global minima that are lower in potential energy than some of the previous candidates. We present equilibrium thermodynamic properties based on parallel tempering simulations, including heat capacities and free energy profiles, which depend explicitly on the crystal structure. We also report equilibrium melting temperatures.

Super-Arrhenius diffusion in an undercooled binary Lennard-Jones liquid results from a quantifiable correlation effect

On short time scales an underlying Arrhenius temperature dependence of the diffusion constant can be extracted from the fragile, super-Arrhenius diffusion of a binary Lennard-Jones mixture. This Arrhenius diffusion is related to the true super-Arrhenius behavior by a factor that depends on the average angle between steps in successive time windows. The correction factor accounts for the fact that on average, successive displacements are negatively correlated, and this effect can therefore be linked directly with the higher apparent activation energy for diffusion at low temperature.

Diagnosing broken ergodicity using an energy fluctuation metric

The Mountain and Thirumalai energy fluctuation metric, Omega(t), has been used to study the effective ergodicity of 60- and 256-atom binary Lennard-Jones mixtures in order to determine the reliability of the calculated diffusion constants at different energies. A plot of Omega(t) against 1time allows the identification of two distinct regimes: ergodic supercooled liquids, where Omega(t) approaches zero, and nonergodic glassy states, where Omega(t) asymptotically approaches a nonzero value on the molecular dynamics time scale. This approach seems to be more appropriate than attempting to define a threshold value for Omega(t)/Omega(0). The behavior of systems between these two limits, which are nonergodic on the time scale considered but may be approaching ergodicity, was examined for a range of simulation times. The calculated diffusion constants change as effective ergodicity is approached, moving closer to the Vogel-Tammann-Fulcher fit defined by higher-energy systems that are already considered to be effectively ergodic. Using the form of the decay of the metric as a measure of ergodicity, we have been able to reproduce the trend in fragility obtained by Sastry for a 256-atom system [Nature (London) 409, 164 (2001)], correcting some of our earlier results [J. Chem. Phys. 120, 8314 (2004)].