Event horizon kinetic Monte Carlo.

In this study, we present a novel approach for modeling the dynamics of stochastic processes. The fundamental concept involves constructing a stochastic Markov chain comprising states separated by more than one stochastic event. Initially, the method explores the network of neighboring states connected by stochastic events. This exploration results in a "horizon" of events leading to a set of "boundary" states at the periphery of each local network. Subsequently, the next member in the Markov chain is selected from the "boundary" states based on the probability of reaching each of the "boundary" states for the first time. Meanwhile, the simulation clock is updated according to the time required to reach the boundary for the first time. This can be achieved using an analytical approach, where the probability of reaching each boundary state for the first time is estimated using absorbing conditions for all boundary states in the analytical solution of a master equation describing the local network of states around each current state. The proposed method is demonstrated in modeling the dynamics in networks of stochastic reactions but can be easily applied in any stochastic system whose dynamics can be expressed via the solution of a master equation. It is expected to enhance the efficiency of event-driven Monte Carlo simulations, originally introduced by Gillespie and widely regarded as the gold standard in the field, especially in cases where the presence of events is characterized by different timescales.

Unscrambling micro-solvation of -COOH and -NH groups in neat dimethyl sulfoxide: insights from 1H-NMR spectroscopy and computational studies.

Dimethyl sulfoxide (DMSO) has a significant, multi-faceted role in medicine, pharmacy, and biology as well as in biophysical chemistry and catalysis. Its physical properties and impact on biomolecular structures still attract major scientific interest, especially the interactions of DMSO with biomolecular functional groups. In the present study, we shed light on the "isolated" carboxylic (-COOH) and amide (-NH) interactions in neat DMSO via1H NMR studies along with extensive theoretical approaches, i.e. molecular dynamics (MD) simulations, density functional theory (DFT), and ab initio calculations, applied on model compounds (i.e. acetic and benzoic acid, ethyl acetamidocyanoacetate). Both experimental and theoretical results show excellent agreement, thereby permitting the calculation of the association constants between the studied compounds and DMSO molecules. Our coupled MD simulations, DFT and ab initio calculations, and NMR spectroscopy results indicated that complex formation is entropically driven and DMSO molecules undergo multiple strong interactions with the studied molecules, particularly with the -COOH groups. The combined experimental and theoretical techniques unraveled the interactions of DMSO with the most abundant functional groups of peptides (i.e. peptide bonds, side chain and terminal carboxyl groups) in high detail, providing significant insights on the underlying thermodynamics driving these interactions. Moreover, the developed methodology for the analysis of the simulation results could serve as a template for future thermodynamic and kinetic studies of similar systems.

Free energy calculations, enhanced by a Gaussian ansatz, for the "chemical work" distribution.

The evaluation of the free energy is essential in molecular simulation because it is intimately related with the existence of multiphase equilibrium. Recently, it was demonstrated that it is possible to evaluate the Helmholtz free energy using a single statistical ensemble along an entire isotherm by accounting for the "chemical work" of transforming each molecule, from an interacting one, to an ideal gas. In this work, we show that it is possible to perform such a free energy perturbation over a liquid vapor phase transition. Furthermore, we investigate the link between a general free energy perturbation scheme and the novel nonequilibrium theories of Crook's and Jarzinsky. We find that for finite systems away from the thermodynamic limit the second law of thermodynamics will always be an inequality for isothermal free energy perturbations, resulting always to a dissipated work that may tend to zero only in the thermodynamic limit. The work, the heat, and the entropy produced during a thermodynamic free energy perturbation can be viewed in the context of the Crooks and Jarzinsky formalism, revealing that for a given value of the ensemble average of the "irreversible" work, the minimum entropy production corresponded to a Gaussian distribution for the histogram of the work. We propose the evaluation of the free energy difference in any free energy perturbation based scheme on the average irreversible "chemical work" minus the dissipated work that can be calculated from the variance of the distribution of the logarithm of the work histogram, within the Gaussian approximation. As a consequence, using the Gaussian ansatz for the distribution of the "chemical work," accurate estimates for the chemical potential and the free energy of the system can be performed using much shorter simulations and avoiding the necessity of sampling the computational costly tails of the "chemical work." For a more general free energy perturbation scheme that the Gaussian ansatz may not be valid, the free energy calculation can be expressed in terms of the moment generating function of the "chemical work" distribution.

Multidimensional direct free energy perturbation.

In this work we propose a multidimensional free energy perturbation scheme that allows the evaluation of the free energy difference between a state sampled based on importance sampling and almost any state that can be constructed by the reduction of the number of molecules in the system and the change of either the interaction energy or the thermodynamic state variable (e.g., the temperature) of the system. We show that via this approach it is possible to evaluate any thermodynamic property included but not limited to free energy, chemical potential, and pressure, along a series of isotherms from a single molecular simulation.

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)].

Accessible Molecular System Creator: Building Molecular Configurations Based on the Inaccessible Molecular Volume and Accessible Molecular Surface via Static Monte Carlo Sampling.

Monte Carlo (MC) stochastic sampling is a powerful tool in classical molecular simulations that directly connects the observable macroscopic properties of matter and the underlying atomistic interactions. This connection operates within the framework of the statistical mechanics proposed by Gibbs. Most MC simulations are "dynamic," creating statistical ensembles of microstates via a Markovian chain, where each microstate in the ensemble depends only on its previous microstate. Herein, we re-examine an alternative form of MC that generates ensemble members through a "static" approach, building molecular systems stepwise. The basic theory for such an approach traces back to Rosenbluth and Rosenbluth, who proposed "static" stepwise sampling of a polymeric chain. It is almost as old as the Metropolis importance sampling approach used in dynamic MC, although the latter has been considerably more popular than the former. Herein, we address the main obstacle in static MC that has hindered the widespread adoption of Rosenbluth-based approaches in atomistic simulations. The obstacle lies in mapping the molecular accessible volume for adding a molecule in a Rosenbluth-like static sampling of atomistic configurations. We demonstrate a breakthrough by leveraging the ability to analytically map the inaccessible molecular volume and the accessible molecular surface owing to interatomically excluded volume interactions. This advance substantially enhances the ability to create molecular samples using a Rosenbluth-like static building process. The proposed approach can be used as a tool for creating initial configurations in MC or molecular dynamics simulations─a field where Rosenbluth-like static building has been applied. Additionally, this approach can be used as the first step in a perturbation scheme that accurately estimates free energy differences by estimating the chemical work related to molecule addition, removal, or reinsertion within the context of free energy perturbation schemes employed in molecular simulations.

Statistical Inference of Rate Constants in Chemical and Biochemical Reaction Networks Using an "Inverse" Event-Driven Kinetic Monte Carlo Method.

The use of rate models for networks of stochastic reactions is frequently used to comprehend the macroscopically observed dynamic properties of finite size reactive systems as well as their relationship to the underlying molecular events. Τhis particular approach usually stumbles on parameter derivation associated with stochastic kinetics, a quite demanding procedure. The present study incorporates a novel algorithm, which infers kinetic parameters from the system's time evolution, manifested as changes in molecular species populations. The proposed methodology reconstructs distributions required to infer kinetic parameters of a stochastic process pertaining to either a simulation or experimental data. The suggested approach accurately replicates rate constants of the stochastic reaction networks, which have evolved over time by event-driven Monte Carlo (MC) simulations using the Gillespie algorithm. Furthermore, our approach has been successfully used to estimate rate constants of association and dissociation events between molecular species developing during molecular dynamics (MD) simulations. We certainly believe that our method will be remarkably helpful for considering the macroscopic characteristic molecular roots related to stochastic physical and biological processes.

Lumping analysis for the prediction of long-time dynamics: from monomolecular reaction systems to inherent structure dynamics of glassy materials.

In this work we develop, test, and implement a methodology that is able to perform, in an automated manner, "lumping" of a high-dimensional, discrete dynamical system onto a lower-dimensional space. Our aim is to develop an algorithm which, without any assumption about the nature of the system's slow dynamics, is able to reproduce accurately the long-time dynamics with minimal loss of information. Both the original and the lumped systems conform to master equations, related via the "lumping" analysis introduced by Wei and Kuo [Ind. Eng. Chem. Fundam. 8, 114 (1969)], and have the same limiting equilibrium probability distribution. The proposed method can be used in a variety of processes that can be modeled via a first order kinetic reaction scheme. Lumping affords great savings in the computational cost and reveals the characteristic times governing the slow dynamics of the system. Our goal is to approach the best lumping scheme with respect to three criteria, in order for the lumped system to be able to fully describe the long-time dynamics of the original system. The criteria used are: (a) the lumping error arising from the reduction process; (b) a measure of the magnitude of singular values associated with long-time evolution of the lumped system; and (c) the size of the lumped system. The search for the optimum lumping proceeds via Monte Carlo simulation based on the Wang-Landau scheme, which enables us to overcome entrapment in local minima in the above criteria and therefore increases the probability of encountering the global optimum. The developed algorithm is implemented to reproduce the long-time dynamics of a glassy binary Lennard-Jones mixture based on the idea of "inherent structures," where the rate constants for transitions between inherent structures have been evaluated via hazard plot analysis of a properly designed ensemble of molecular dynamics trajectories.

On the hierarchical design of biochemical-based digital computations.

The understanding of the biochemical processes underpinning various biological systems has significantly increased in recent decades and has even prompted reverse engineering of certain of life's more complex processes. The most prominent example is modern computers designed to mimic neuron activity. This work forms part of growing endeavors to return advances in the theory of computation and electronics to biology. In this context, we present a set of requirements sufficient for the design of biochemical analogs of modern electronics in a hierarchical, modular fashion that mimics the design of modern computational devices. This theoretical approach is based on a simple enzymatic analog of the transistor and supported by numerical simulations of biochemical models of enzymatic networks equivalent to complex, and modular, interconnecting electronic circuitry (including clocks, Flip-Flops, adders, decoders, and multiplexers). Furthermore, the proposed approach has been implemented in the form of a Python library capable of creating and testing models of complex bio-analog digital computations based on the execution of an elementary universal logic gate. In tribute to Claude Shannon, our biochemical network materializes his example of a "password" recognition that moves the language of the modern theory of automata beyond combinatorial logic and towards sequential logic.

Probing subglass relaxation in polymers via a geometric representation of probabilities, observables, and relaxation modes for discrete stochastic systems.

The dynamics of many physical, chemical, and biological systems can be reduced to a succession of infrequent transitions in a network of discrete states representing low energy regions in configuration space. This enables accessing long-time dynamics and predicting macroscopic properties. Here we develop a new, perfectly general statistical mechanical/geometric formulation that expresses both state probabilities and all observables in the same Euclidean space, spanned by the eigenvectors of the symmetrized time evolution operator. Our formalism leads to simple expressions for nonequilibrium and equilibrium ensemble averages, variances, and time correlation functions of any observable and allows a rigorous decomposition of the dynamics into relaxation modes. Applying it to subglass segmental relaxation in atactic polystyrene up to times on the order of 10 micros, we probe the molecular mechanism of the gamma and delta processes and unequivocally identify the delta process with rotation of a single phenyl group around its stem.

On the role of inherent structures in glass-forming materials: II. Reconstruction of the mean square displacement by rigorous lifting of the inherent structure dynamics.

In a previous paper, we investigated the role of inherent structures in the vitrification process of glass-forming materials, showing that the dynamical transitions between inherent structures (states) can be well predicted by a first-order kinetic scheme based on infrequent-event theory at low temperatures. In this work, we utilize and extend that methodology in order to completely reconstruct the system dynamics in the form of the mean square atomic displacement as a function of time at finite temperatures on the basis of the succession of transitions in a network of states, the vibrational contribution being evaluated on the basis of short molecular dynamics runs artificially trapped within each one of the states. In order to do so, we provide the mathematical formulation for lifting the coarse-grained Poisson process model of transitions between states back to the atomistic level and thereby reproducing the full dynamics of the atomistic system within the Poisson approximation. Our result shows excellent agreement for temperatures around and below the glass transition temperature of our model Lennard-Jones two-component mixtures. Clearly, our approach is able to reproduce the full dynamics of the atomistic system at low temperatures, where the Poisson approximation is valid.

On the role of inherent structures in glass-forming materials: I. The vitrification process.

In this work, we investigate the role of inherent structures in the vitrification process of glass-forming materials by using a two-component Lennard-Jones mixture. We start with a simplified model that describes the dynamics of the atomistic system as a Poisson process consisting of a series of transitions from one potential energy minimum (inherent structure) to another, the rate of individual transitions being described by a first-order kinetic law. We investigate the validity of this model by comparing the mean square displacement resulting from atomistic molecular dynamics (MD) trajectories with the corresponding mean square displacement based on inherent structures. Furthermore, in the case of vitrification via stepwise cooling, we identify the role of the potential energy landscape in determining the properties of the resulting glass. Interestingly, the cooling rate is not sufficient to define the resulting glass in a stepwise cooling process, because the time spent by the system at different temperatures (length of the steps) has a highly nonlinear impact on the properties of the resulting glass. In contrast to previous investigations of supercooled liquids, we focus on a range of temperatures close to and below the glass transition temperature, where the use of MD is incapable of producing equilibrated samples of the metastable supercooled state. Our aim is to develop a methodology that enables mapping the dynamics under these conditions to a coarse-grained first-order kinetic model based on the Poisson process approximation. This model can be used in order to extend our dynamical sampling ability to much broader time scales and therefore allow us to create computer glasses with cooling rates closer to those used experimentally. In a continuation to this work, we provide the mathematical formulation for lifting the coarse-grained Poisson process model and reproducing the full dynamics of the atomistic system.

On the estimation of the free energy, from a single equilibrium statistical ensemble, via particle reinsertion.

In this work, we propose the evaluation of the free energy in molecular systems, in a "single" step, by "deleting" all the molecules in the system. The approach can be considered as the statistical mechanics analogue of the evaluation of the potential energy in classical mechanics by accounting for the necessary work to transfer all particles one by one to infinite distance. As a result, the free energy of an atomistic system can now be expressed as an ensemble average over a configurational function that corresponds to the contribution of each microstate to the free energy of the ensemble. Moreover, the proposed method is capable of evaluating the free energy as a function of the density, from the simulated density down to zero. Finally, the proposed method is related to the Rosenbluth sampling of the inverse process, that of inserting (instead of deleting) and provide the analogous theorems to Bennett's and Crooks' work (Bennett, C. H. J. Comput. Phys. 1976, 22, 245; Crooks, G. E. Phys. Rev. E 1999, 60, 2721). When the proposed process is envisioned as the transformation of an interacting to a noninteracting system, the proposed scheme reduces to the Jarzynski identity linking the free energy of the system to the chemical work related to this transformation.

Prediction of infinite dilution benzene solubility in linear polyethylene melts via the direct particle deletion method.

The solubility of benzene in linear polyethylene melts was estimated via Monte Carlo simulations using a united-atom molecular model at temperatures between 373 and 573 K, in the infinite dilution limit. The excess chemical potential of the solute was evaluated with the direct particle deletion (DPD) method, whose rigorous derivation is presented here in detail: in this scheme, the benzene molecule united atoms are converted to hard spheres and then removed from the polymer system. The simulations were carried out in the N(1)N(2)PT ensemble using advanced Monte Carlo moves to equilibrate the polymeric phase. The evaluation of the accessible volume fraction for the "hard sphere" solute molecule required by the DPD method was performed analytically. The effect of the value of the arbitrary hard sphere diameter, d, on the computed thermodynamic quantities was determined, allowing us to establish an optimal range for the system considered. The values of Henry's law constant are in good agreement with experimental data from the literature in the temperature range considered and are comparable to those obtained with the lattice fluid and PC(SAFT) equations of state for the same system.

Temperature accelerated dynamics in glass-forming materials.

In this work we propose a methodology for improving dynamical sampling in molecular simulations via temperature acceleration. The proposed approach combines the novel methods of Voter for temperature-accelerated dynamics with the multiple histogram reweighting method of Ferrenberg and Swendsen, its dynamical extension by Nieto-Draghi et al., and with hazard plot analysis, allowing optimal sampling with small computational cost over time scales inaccessible to classical molecular dynamics simulations by utilizing the "inherent structure" idea. The time evolution of the system is viewed as a succession of transitions between "basins" in its potential energy landscape, each basin containing a local minimum of the energy (inherent structure). Applying the proposed algorithm to a glass-forming material consisting of a mixture of spherical atoms interacting via Lennard-Jones potentials, we show that it is possible to perform an exhaustive search and evaluate rate constants for basin-to-basin transitions that cover several orders of magnitude on the time scale, far beyond the abilities of any competitive dynamical study, revealing an extreme ruggedness of the potential energy landscape in the vicinity of the glass transition temperature. By analyzing the network of inherent structures, we find that there are characteristic distances and rate constants related to the dynamical entrapment of the system in a neighborhood of basins (a metabasin), whereas evidence to support a random energy model is provided. The multidimensional configurational space displays a self-similar character depicted by a fractal dimension around 2.7 (+/-0.5) for the set of sampled inherent structures. Only transitions with small Euclidean measure can be considered as localized.

Dynamical integration of a Markovian web: a first passage time approach.

In this work we address the dynamics of Markovian systems by tracking the evolution of the probability distribution, utilizing mean first passage time theory to augment the set of states considered. The method is validated on a lattice system and is applied, in conjunction with landscape analysis (saddle point searches) and multidimensional transition-state theory, to an atomistic model of glassy atactic polystyrene, in order to follow its time evolution over more than ten orders of magnitude on the time scale, from less than 10(-15) up to 10(-5) s. Frequencies extracted from the eigenvalues of the rate constant matrix are in favorable agreement with experimental measurements of subglass relaxation transitions at 250 K.

Monte Carlo Sampling of a Markov Web.

The efficiency of Markov-Chain Monte Carlo simulations can be enhanced by exploiting information about trial moves that would normally be rejected. The original presentation of this approach was limited to a specific MC sampling scheme. Here we present a general derivation of a method to improve the sampling efficiency of Monte Carlo simulations by collecting information about the microstates that can be linked directly to the sampled point via an independent Markov transition matrix. As an illustration, we show that our approach greatly enhances the efficiency of a scheme to compute the density of states of a square-well fluid.

Novel Monte Carlo scheme for systems with short-ranged interactions.

We propose a Monte Carlo (MC) sampling algorithm to simulate systems of particles interacting via very short-ranged discontinuous potentials. Such models are often used to describe protein solutions or colloidal suspensions. Most normal MC algorithms fail for such systems because, at low temperatures, they tend to get trapped in local potential-energy local minima due to the short range of the pair potential. To circumvent this problem, we have devised a scheme that changes the construction of trial moves in such a way that the potential-energy difference between initial and final states drops out of the acceptance rule for the Monte Carlo trial moves. This approach allows us to simulate systems with short-ranged attraction under conditions that were unreachable up to now.