Measurement of nonequilibrium entropy from space-time thermodynamic integration

The entropy of a system transiently driven out of equilibrium by a time-inhomogeneous stochastic dynamics is first expressed as a transient response function generalizing the nonlinear Kawasaki-Crooks response. This function is then reformulated into three statistical averages defined over ensembles of nonequilibrium trajectories. The first average corresponds to a space-time thermodynamic perturbation relation, while the two following ones correspond to space-time thermodynamic integration relations. Provided that trajectories are initiated starting from a distribution of states that is analytically known, the ensemble averages are computationally amenable to Markov chain Monte Carlo methods. The relevance of importance sampling in path ensembles is confirmed in practice by computing the nonequilibrium entropy of a driven toy system. We finally study a situation where the dynamics produces entropy. In this case, we observe that space-time thermodynamic integration still yields converged estimates, while space-time thermodynamic perturbation turns out to converge very slowly.     

A targeted reweighting method for accelerating the exploration of high-dimensional configuration space

Time scales available to biomolecular simulations are limited by barriers among states in a high-dimensional configuration space. If equilibrium averages are to be computed, methods that accelerate barrier passage can be carried out by non-Boltzmann sampling. Barriers can be reduced by modifying the potential-energy function and running dynamics on the modified surface. The Boltzmann average can be restored by reweighting each point along the trajectory. We introduce a targeted reweighting scheme where some barriers are reduced, while others are not modified. If only equilibrium properties are desired, trajectories in configuration space can be generated by Langevin dynamics. Once past a transient time, these trajectories guarantee equilibrium sampling when reweighted. A relatively high-order stochastic integration method can be used to generate trajectories. The targeted reweighting scheme is illustrated by a series of double-well models with varying degrees of freedom and shown to be a very efficient method to provide the correct equilibrium distributions, in comparison with analytic results. The scheme is applied to a protein model consisting of a chain of connected beads characterized by dihedral angles and the van der Waals interactions among the beads. We investigate the sampling of configuration space for a model of a helix-turn-helix motif. The targeted reweighting is found to be essential to permit the original all-helical conformation to bend and generate turn structures while still maintaining the alpha-helical segments.     

Molecular latent space simulators

Small integration time steps limit molecular dynamics (MD) simulations to millisecond time scales. Markov state models (MSMs) and equation-free approaches learn low-dimensional kinetic models from MD simulation data by performing configurational or dynamical coarse-graining of the state space. The learned kinetic models enable the efficient generation of dynamical trajectories over vastly longer time scales than are accessible by MD, but the discretization of configurational space and/or absence of a means to reconstruct molecular configurations precludes the generation of continuous atomistic molecular trajectories. We propose latent space simulators (LSS) to learn kinetic models for continuous atomistic simulation trajectories by training three deep learning networks to (i) learn the slow collective variables of the molecular system, (ii) propagate the system dynamics within this slow latent space, and (iii) generatively reconstruct molecular configurations. We demonstrate the approach in an application to Trp-cage miniprotein to produce novel ultra-long synthetic folding trajectories that accurately reproduce atomistic molecular structure, thermodynamics, and kinetics at six orders of magnitude lower cost than MD. The dramatically lower cost of trajectory generation enables greatly improved sampling and greatly reduced statistical uncertainties in estimated thermodynamic averages and kinetic rates.

Latent space simulators learn kinetic models for atomistic simulations and generate novel trajectories at six orders of magnitude lower cost.     

New insights into the meaning and usefulness of principal component analysis of concatenated trajectories

A comparison between different conformations of a given protein, relating both structure and dynamics, can be performed in terms of combined principal component analysis (combined‐PCA). To that end, a trajectory is obtained by concatenating molecular dynamics trajectories of the individual conformations under comparison. Then, the principal components are calculated by diagonalizing the correlation matrix of the concatenated trajectory. Since the introduction of this approach in 1995 it has had a large number of applications. However, the interpretation of the eigenvectors and eigenvalues so obtained is based on intuitive foundations, because analytical expressions relating the concatenated correlation matrix with those of the individual trajectories under consideration have not been provided yet. In this article, we present such expressions for the cases of two, three, and an arbitrary number of concatenated trajectories. The formulas are simple and show what is to be expected and what is not to be expected from a combined‐PCA. Their correctness and usefulness is demonstrated by discussing some representative examples. The results can be summarized in a simple sentence: the correlation matrix of a concatenated trajectory is given by the average of the individual correlation matrices plus the correlation matrix of the individual averages. From this it follows that the combined‐PCA of trajectories belonging to different free energy basins provides information that could also be obtained by alternative and more straightforward means. © 2014 Wiley Periodicals, Inc.     

Biased sampling of nonequilibrium trajectories: can fast switching simulations outperform conventional free energy calculation methods?

We have investigated the maximum computational efficiency of reversible work calculations that change control parameters in a finite amount of time. Because relevant nonequilibrium averages are slow to converge, a bias on the sampling of trajectories can be beneficial. Such a bias, however, can also be employed in conventional methods for computing reversible work, such as thermodynamic integration or umbrella sampling. We present numerical results for a simple one-dimensional model and for a Widom insertion in a soft sphere liquid, indicating that, with an appropriately chosen bias, conventional methods are in fact more efficient. We describe an analogy between nonequilibrium dynamics and mappings between equilibrium ensembles, which suggests that the practical inferiority of fast switching is quite general. Finally, we discuss the relevance of adiabatic invariants in slowly driven Hamiltonian systems for the application of Jarzynski's theorem.     

On the relationships between kinetic schemes and two-state single molecule trajectories

Trajectories of a signal that fluctuates between two states which originate from single molecule activities have become ubiquitous. Common examples are trajectories of ionic flux through individual membrane channels and of photon counts collected from diffusion, activity, and conformational changes of biopolymers. By analyzing the trajectory, one wishes to deduce the underlying mechanism, which is usually described by a multisubstate kinetic scheme. In previous works [O. Flomenborn, J. Klafter, and A. Szabo, Biophys. J. 88, 3780 (2005); O. Flomenbom and J. Klafter, Acta Phys. Pol. B 36, 1527 (2005)], we divided kinetic schemes that generate two-state trajectories into two types: reducible schemes and irreducible schemes. A full characterization of the reducible ones was given. We showed that all the information in trajectories generated from reducible schemes is contained in the waiting time probability density functions (PDFs) of the two states. It follows that reducible schemes with the same waiting time PDFs are not distinguishable; namely, such schemes lead to identical two-state trajectories in the statistical sense. In this work, we further characterize the topologies of kinetic schemes, now of irreducible schemes, and further study two-state trajectories from the two types of scheme. We suggest various methods for extracting information about the underlying kinetic scheme from the trajectory (e.g., calculate the binned successive waiting times PDFs and analyze the ordered waiting time trajectories), and point out the advantages and disadvantages of each. We show that the binned successive waiting times PDFs are not only more robust than other functions when analyzing finite trajectories, but contain, in most cases, more information about the underlying kinetic scheme than other functions in the limit of infinitely long trajectories. For some cases, however, analyzing the ordered waiting times trajectory may supply unique information about the underlying kinetic scheme.     

Fluctuations in the ensemble of reaction pathways

The dominant reaction pathway is a rigorous framework to microscopically compute the most probable trajectories, in nonequilibrium transitions. In the low-temperature regime, such dominant pathways encode the information about the reaction mechanism and can be used to estimate nonequilibrium averages of arbitrary observables. On the other hand, at sufficiently high temperatures, the stochastic fluctuations around the dominant paths become important and have to be taken into account. In this work, we develop a technique to systematically include the effects of such stochastic fluctuations, to order k(B)T. This method is used to compute the probability for a transition to take place through a specific reaction channel and to evaluate the reaction rate.     

Theoretical unimolecular kinetics for CH4 + M ⇄ CH3 + H + M in eight baths, M = He, Ne, Ar, Kr, H2, N2, CO, and CH4

Ensembles of classical trajectories are used to study collisional energy transfer in highly vibrationally excited CH(4) for eight bath gases. Several simplifying assumptions for the CH(4) + M interaction potential energy surface are tested against full dimensional direct dynamics trajectory calculations for M = He, Ne, and H(2). The calculated energy transfer averages are confirmed to be sensitive to the shape of the repulsive wall of the intermolecular potential, with an exponential repulsive wall required for quantitative predictions. For the diatomic baths, the usual "separable pairwise" approximation for the interaction potential is unable to describe the orientation dependence of the interaction potential accurately, and the ambiguity in the resulting parametrizations contributes an additional uncertainty to the predicted energy transfer averages of 20-40%. On the other hand, the energy transfer averages are shown to be insensitive to the level of theory used to describe the intramolecular CH(4) potential, with a computationally efficient semiempirical tight binding potential for hydrocarbons performing equally well as an MP2 potential. The relative collisional energy transfer efficiencies of the eight bath gases are discussed and shown to be a function of temperature. The ensemble-averaged energy transferred in deactivating collisions <ΔE(d)> for each bath is used to parametrize a single-exponential-down model for collisional energy transfer in master equation calculations. The predicted decomposition rate coefficients for CH(4) agree well with available experimental rate coefficients for M = He, Ar, Kr, and CH(4). The effect of vibrational anharmonicity on the predicted rate coefficients is considered briefly.     

Dynamical reweighting: improved estimates of dynamical properties from simulations at multiple temperatures

Dynamical averages based on functionals of dynamical trajectories, such as time-correlation functions, play an important role in determining kinetic or transport properties of matter. At temperatures of interest, the expectations of these quantities are often dominated by contributions from rare events, making the precise calculation of these quantities by molecular dynamics simulation difficult. Here, we present a reweighting method for combining simulations from multiple temperatures (or from simulated or parallel tempering simulations) to compute an optimal estimate of the dynamical properties at the temperature of interest without the need to invoke an approximate kinetic model (such as the Arrhenius law). Continuous and differentiable estimates of these expectations at any temperature in the sampled range can also be computed, along with an assessment of the associated statistical uncertainty. For rare events, aggregating data from multiple temperatures can produce an estimate with the desired precision at greatly reduced computational cost compared with simulations conducted at a single temperature. Here, we describe use of the method for the canonical (NVT) ensemble using four common models of dynamics (canonical distribution of Hamiltonian trajectories, Andersen thermostatting, Langevin, and overdamped Langevin or Brownian dynamics), but it can be applied to any thermodynamic ensemble provided the ratio of path probabilities at different temperatures can be computed. To illustrate the method, we compute a time-correlation function for solvated terminally-blocked alanine peptide across a range of temperatures using trajectories harvested using a modified parallel tempering protocol.     

Temporal correlation functions of concentration fluctuations: an anomalous case

We calculate, within the framework of the continuous time random walk (CTRW) model, multiparticle temporal correlation functions of concentration fluctuations (CCF) in systems that display anomalous subdiffusion. The subdiffusion stems from the nonstationary nature of the CTRW waiting times, which also lead to aging and ergodicity breaking. Due to aging, a system of diffusing particles tends to slow down as time progresses, and therefore, the temporal correlation functions strongly depend on the initial time of measurement. As a consequence, time averages of the CCF differ from ensemble averages, displaying therefore ergodicity breaking. We provide a simple example that demonstrates the difference between these two averages, a difference that might be amenable to experimental tests. We focus on the case of ensemble averaging and assume that the preparation time of the system coincides with the starting time of the measurement. Our analytical calculations are supported by computer simulations based on the CTRW model.     

TRAVIS - a free analyzer and visualizer for Monte Carlo and molecular dynamics trajectories

We present TRAVIS ("TRajectory Analyzer and VISualizer"), a free program package for analyzing and visualizing Monte Carlo and molecular dynamics trajectories. The aim of TRAVIS is to collect as many analyses as possible in one program, creating a powerful tool and making it unnecessary to use many different programs for evaluating simulations. This should greatly rationalize and simplify the workflow of analyzing trajectories. TRAVIS is written in C++, open-source freeware and licensed under the terms of the GNU General Public License (GPLv3). It is easy to install (platform independent, no external libraries) and easy to use. In this article, we present some of the algorithms that are implemented in TRAVIS - many of them widely known for a long time, but some of them also to appear in literature for the first time. All shown analyses only require a standard MD trajectory as input data.     

Extension of the quantum virial and Hellmann–Feynman theorems to the quantum statistical averages

In the present article, we extended the quantum virial and Hellmann–Feynman theorems to the quantum statistical averages, that is, to the thermal states. We obtained some new formulas which make possible expressing the thermodynamical observables of the system as functions of the moments of coordinates, as we see in a short example relating to the pseudoharmonical oscillator. © 1998 John Wiley & Sons, Inc. Int J Quant Chem 69: 159–165, 1998     

Using Selectively Applied Accelerated Molecular Dynamics to Enhance Free Energy Calculations

Accelerated molecular dynamics (aMD) has been shown to enhance conformational space sampling relative to classical molecular dynamics; however, the exponential reweighting of aMD trajectories, which is necessary for the calculation of free energies relating to the classical system, is oftentimes problematic, especially for systems larger than small poly peptides. Here, we propose a method of accelerating only the degrees of freedom most pertinent to sampling, thereby reducing the total acceleration added to the system and improving the convergence of calculated ensemble averages, which we term selective aMD. Its application is highlighted in two biomolecular cases. First, the model system alanine dipeptide is simulated with classical MD, all-dihedral aMD, and selective aMD, and these results are compared to the infinite sampling limit as calculated with metadynamics. We show that both forms of aMD enhance the convergence of the underlying free energy landscape by 5-fold relative to classical MD; however, selective aMD can produce improved statistics over all-dihedral aMD due to the improved reweighting. Then we focus on the pharmaceutically relevant case of computing the free energy of the decoupling of oseltamivir in the active site of neuraminidase. Results show that selective aMD greatly reduces the cost of this alchemical free energy transformation, whereas all-dihedral aMD produces unreliable free energy estimates.     

The beta-peptide hairpin in solution: conformational study of a beta-hexapeptide in methanol by NMR spectroscopy and MD simulation.

The structural and thermodynamic properties of a 6-residue beta-peptide that was designed to form a hairpin conformation have been studied by NMR spectroscopy and MD simulation in methanol solution. The predicted hairpin would be characterized by a 10-membered hydrogen-bonded turn involving residues 3 and 4, and two extended antiparallel strands. The interproton distances and backbone torsional dihedral angles derived from the NMR experiments at room temperature are in general terms compatible with the hairpin conformation. Two trajectories of system configurations from 100-ns molecular-dynamics simulations of the peptide in solution at 298 and 340 K have been analyzed. In both simulations reversible folding to the hairpin conformation is observed. Interestingly, there is a significant conformational overlap between the unfolded state of the peptide at each of the temperatures. As already observed in previous studies of peptide folding, the unfolded state is composed of a (relatively) small number of predominant conformers and in this case lacks any type of secondary-structure element. The trajectories provide an excellent ground for the interpretation of the NMR-derived data in terms of ensemble averages and distributions as opposed to single-conformation interpretations. From this perspective, a relative population of the hairpin conformation of 20% to 30% would suffice to explain the NMR-derived data. Surprisingly, however, the ensemble of structures from the simulation at 340 K reproduces more accurately the NMR-derived data than the ensemble from the simulation at 298 K, a question that needs further investigation.     

In–out asymmetry and interference effects in plasmon excitation by swift charged particles traversing a surface

We investigate the excitation of plasmons by a fast charged particle moving in the vicinity or traversing the surface of a solid along an arbitrary trajectory. We use both quantum‐mechanical and semiclassical dielectric formulations to study how the particle couples with the bulk and surface excitations. We pay special attention to the differences and similarities between incoming and outgoing trajectories, finding some novel oscillatory structures that can be ascribed to an interference effect between direct and reflected plasmon excitations. Copyright © 2013 John Wiley & Sons, Ltd.     

Temperature-induced structural transitions and vibrational properties of microclusters: ab-initio molecular dynamics studies

We present the results of ab-initio molecular dynamics studies of selected microclusters of sodium, silicon and magnesium at finite temperatures, and especially discuss those obtained around room temperature. In particular, from the analysis of the atomic trajectories we can identify in some cases the existence of different isomers and the isomerization pathways. We have also calculated vibrational spectra at low temperatures and find that they can be used as a very sensitive structural probe also in sodium clusters, where the electronic properties are quite insensitive to the geometry.     

Some properties of ion trajectories in the radial plane of an axially symmetric ion trap

In this article we analyse the trajectories of externally generated ions injected in the radial plane into the ion trap. The shape of the envelope curves for two-dimensional (2D) ion trajectories is determined. Conditions under which these envelope curves can be transformed into circles are found. We show that the amplitude of ion oscillations is a minimum in this case and that this mode corresponds to optimised ion trapping conditions. Also we discuss a ring-shaped ion trap mass spectrometer electrode system which consists of two ring electrodes, and which utilises ion trajectories with circular envelope curves.     

Dynamic effects on j-couplings across hydrogen bonds in proteins

We combine molecular dynamics simulation (MD) with density functional theory (DFT)/finite perturbation theory (FPT) to obtain the Fermi contact contributions to the 3hJNC' couplings in the SMN Tudor domain. Our results show that the effect of conformational motion needs to be considered for the accurate prediction of these scalar couplings. For hydrogen bonds in the beta-sheet regions, the calculated cumulative J-coupling averages remain constant after the first 200 ps. In the more flexible regions at the edges of the beta-sheets, the cumulative J-coupling averages vary over the 500-ps trajectory, providing a qualitative insight into the degree of conformational motion in different regions of the protein.     

Modal dynamics of proteins in water

We studied a dynamical model for the motion of the large scales of proteins in water. The model was obtained by projecting the (averaged) Newton equations onto some set of harmonic modes. We compared the statistics of the so‐obtained trajectories with those obtained by standard techniques, and concluded that our dynamical model is able to fairly reproduce the average properties of the large‐scale motion of the protein, and at the same time allow time steps one order of magnitude larger than the standard ones. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1274–1282, 2000