Highly parallelizable path sampling with minimal rejections using asynchronous replica exchange and infinite swaps.

Capturing rare yet pivotal events poses a significant challenge for molecular simulations. Path sampling provides a unique approach to tackle this issue without altering the potential energy landscape or dynamics, enabling recovery of both thermodynamic and kinetic information. However, despite its exponential acceleration compared to standard molecular dynamics, generating numerous trajectories can still require a long time. By harnessing our recent algorithmic innovations-particularly subtrajectory moves with high acceptance, coupled with asynchronous replica exchange featuring infinite swaps-we establish a highly parallelizable and rapidly converging path sampling protocol, compatible with diverse high-performance computing architectures. We demonstrate our approach on the liquid-vapor phase transition in superheated water, the unfolding of the chignolin protein, and water dissociation. The latter, performed at the ab initio level, achieves comparable statistical accuracy within days, in contrast to a previous study requiring over a year.

PyRETIS 3: Conquering rare and slow events without boundaries.

We present and discuss the advancements made in PyRETIS 3, the third instalment of our Python library for an efficient and user-friendly rare event simulation, focused to execute molecular simulations with replica exchange transition interface sampling (RETIS) and its variations. Apart from a general rewiring of the internal code towards a more modular structure, several recently developed sampling strategies have been implemented. These include recently developed Monte Carlo moves to increase path decorrelation and convergence rate, and new ensemble definitions to handle the challenges of long-lived metastable states and transitions with unbounded reactant and product states. Additionally, the post-analysis software PyVisa is now embedded in the main code, allowing fast use of machine-learning algorithms for clustering and visualising collective variables in the simulation data.

Path sampling with memory reduction and replica exchange to reach long permeation timescales.

Assessing kinetics in biological processes with molecular dynamics simulations remains a computational and conceptual challenge, given the large time and length scales involved. For kinetic transport of biochemical compounds or drug molecules, the permeability through the phospholipid membranes is a key kinetic property, but long timescales are hindering the accurate computation. Technological advances in high-performance computing therefore need to be accompanied by theoretical and methodological developments. In this contribution, the replica exchange transition interface sampling (RETIS) methodology is shown to give perspective toward observing longer permeation pathways. It is first reviewed how RETIS, a path-sampling methodology that gives in principle exact kinetics, can be used to compute membrane permeability. Next, recent and current developments in three RETIS aspects are discussed: several new Monte Carlo moves in the path-sampling algorithm, memory reduction by reducing pathlengths, and exploitation of parallel computing with CPU-imbalanced replicas. Finally, the memory reduction presenting a new replica exchange implementation, coined REPPTIS, is showcased with a permeant needing to pass a membrane with two permeation channels, either representing an entropic or energetic barrier. The REPPTIS results showed clearly that inclusion of some memory and enhancing ergodic sampling via replica exchange moves are both necessary to obtain correct permeability estimates. In an additional example, ibuprofen permeation through a dipalmitoylphosphatidylcholine membrane was modeled. REPPTIS succeeded in estimating the permeability of this amphiphilic drug molecule with metastable states along the permeation pathway. In conclusion, the presented methodological advances allow for deeper insight into membrane biophysics even if the pathways are slow, as RETIS and REPPTIS push the permeability calculations to longer timescales.

Enhanced path sampling using subtrajectory Monte Carlo moves.

Path sampling allows the study of rare events, such as chemical reactions, nucleation, and protein folding, via a Monte Carlo (MC) exploration in path space. Instead of configuration points, this method samples short molecular dynamics (MD) trajectories with specific start- and end-conditions. As in configuration MC, its efficiency highly depends on the types of MC moves. Since the last two decades, the central MC move for path sampling has been the so-called shooting move in which a perturbed phase point of the old path is propagated backward and forward in time to generate a new path. Recently, we proposed the subtrajectory moves, stone-skipping (SS) and web-throwing, that are demonstrably more efficient. However, the one-step crossing requirement makes them somewhat more difficult to implement in combination with external MD programs or when the order parameter determination is expensive. In this article, we present strategies to address the issue. The most generic solution is a new member of subtrajectory moves, wire fencing (WF), that is less thrifty than the SS but more versatile. This makes it easier to link path sampling codes with external MD packages and provides a practical solution for cases where the calculation of the order parameter is expensive or not a simple function of geometry. We demonstrate the WF move in a double-well Langevin model, a thin film breaking transition based on classical force fields, and a smaller ruthenium redox reaction at the ab initio level in which the order parameter explicitly depends on the electron density.

Exchanging Replicas with Unequal Cost, Infinitely and Permanently.

We developed a replica exchange method that is effectively parallelizable even if the computational cost of the Monte Carlo moves in the parallel replicas are considerably different, for instance, because the replicas run on different types of processor units or because of the algorithmic complexity. To prove detailed-balance, we make a paradigm shift from the common conceptual viewpoint in which the set of parallel replicas represents a high-dimensional superstate, to an ensemble-based criterion in which the other ensembles represent an environment that might or might not participate in the Monte Carlo move. In addition, based on a recent algorithm for computing permanents, we effectively increase the exchange rate to infinite without the steep factorial scaling as a function of the number of replicas. We illustrate the effectiveness of this replica exchange methodology by combining it with a quantitative path sampling method, replica exchange transition interface sampling (RETIS), in which the costs for a Monte Carlo move can vary enormously as paths in a RETIS algorithm do not have the same length and the average path lengths tend to vary considerably for the different path ensembles that run in parallel. This combination, coined ∞RETIS, was tested on three model systems.

The role of pressure and defects in the wurtzite to rock salt transition in cadmium selenide.

Using molecular dynamics and path sampling techniques we investigated the effect of pressure and defects in the wurtzite to rock salt transition in cadmium selenide (CdSe). In the pressure range 2-10 GPa, rate constants of transition are in the order of 10-23 to 105 s-1 for the transformation of a relatively small wurtzite crystal consisting of 1024 atoms with periodic boundary conditions. The transition paths predominantly evolve through an intermediate 5-coordinated structure, as reported before, though its typical lifetime within the transition paths is particularly long in the intermediate pressure range (4-6 GPa). The defects were created by removing Cd-Se pairs from an otherwise perfect crystal. The removals were either selected fully randomized or grouped in clusters (cavity creation). We find that the rate of transition due to the defects increases by several orders of magnitude even for a single pair removal. This is caused by a change in the transition mechanism that no longer proceeds via the intermediate 5-coordinated structure, when defects are present. Further, the cavity creation yields a lower rate than the fully randomized removal.

Permeation Rates of Oxygen through a Lipid Bilayer Using Replica Exchange Transition Interface Sampling.

Several simulations strategies have emerged to predict the permeability of solutes across membranes, which is important for many biological or industrial processes such as drug design. The widespread inhomogeneous solubility-diffusion (ISD) model is based on the Smoluchowski equation and describes permeation as purely diffusive. The counting method, which counts membrane transitions in a long molecular dynamics (MD) trajectory, is free of this diffusive assumption, but it lacks sufficient statistics when the permeation involves high free energy barriers. Metadynamics and variations thereof can overcome such barriers, but they generally lack the kinetics information. The milestoning framework has been used to describe permeation as a rare event, but it still relies on the Markovian assumption between the milestones. Replica Exchange Transition Interface Sampling (RETIS) has been shown to be an effective method for sampling rare events while simultaneously describing the kinetics without assumptions. This paper is the first permeation application of RETIS on an all-atom lipid bilayer consisting of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) to compute the entrance, escape and complete transition of molecular oxygen. Conventional MD was performed as a benchmark, and the MD rates from counting were converted to rate constants, giving good agreement with the RETIS values. Moreover, a correction factor was derived to convert the collective order parameter in RETIS, which was aimed to improve efficiency, to a single-particle order parameter. With this work, we showed how the exact kinetics of drug molecules permeation can be assessed with RETIS even if the permeation is truly a rare event or if the permeation is non-Markovian. RETIS will therefore be a valuable tool for future permeation studies.

PyRETIS 2: An improbability drive for rare events.

The algorithmic development in the field of path sampling has made tremendous progress in recent years. Although the original transition path sampling method was mostly used as a qualitative tool to sample reaction paths, the more recent family of interface-based path sampling methods has paved the way for more quantitative rate calculation studies. Of the exact methods, the replica exchange transition interface sampling (RETIS) method is the most efficient, but rather difficult to implement. This has been the main motivation to develop the open-source Python-based computer library PyRETIS that was released in 2017. PyRETIS is designed to be easily interfaced with any molecular dynamics (MD) package using either classical or ab initio MD. In this study, we report on the principles and the software enhancements that are now included in PyRETIS 2, as well as the recent developments on the user interface, improvements of the efficiency via the implementation of new shooting moves, easier initialization procedures, analysis methods, and supported interfaced software. © 2019 The Authors. Journal of Computational Chemistry published by Wiley Periodicals, Inc.

Multiscale partial charge estimation on graphene for neutral, doped and charged flakes.

The minimal-basis iterative stockholder (MBIS) and restrained electrostatic potential (RESP) methods were applied to examine the effects of edges and of nitrogen and boron dopants on the atomic partial charges of neutral and charged graphene flakes. The results provided the parameters to fit a second-order atom-condensed Kohn-Sham DFT model (ACKS2), accurately determining the partial charges, the dipole and local electric fields in large graphene flakes with negligible cost. Our approach can lead to improvements of graphene force fields in charged conditions and guide the design of media for catalytic applications.

Diffusion of gas mixtures in the sI hydrate structure.

Replacing methane with carbon dioxide in gas hydrates has been suggested as a way of harvesting methane, while at the same time storing carbon dioxide. Experimental evidence suggests that this process is facilitated if gas mixtures are used instead of pure carbon dioxide. We studied the free energy barriers for diffusion of methane, carbon dioxide, nitrogen, and hydrogen in the sI hydrate structure using molecular simulation techniques. Cage hops between neighboring cages were considered with and without a water vacancy and with a potential inclusion of an additional gas molecule in either the initial or final cage. Our results give little evidence for enhanced methane and carbon dioxide diffusion if nitrogen is present as well. However, the inclusion of hydrogen seems to have a substantial effect as it diffuses rapidly and can easily enter occupied cages, which reduces the barriers of diffusion for the gas molecules that co-occupy a cage with hydrogen.

Local initiation conditions for water autoionization.

The pH of liquid water is determined by the infrequent process in which water molecules split into short-lived hydroxide and hydronium ions. This reaction is difficult to probe experimentally and challenging to simulate. One of the open questions is whether the local water structure around a slightly stretched OH bond is actually initiating the eventual breakage of this bond or whether this event is driven by a global ordering that involves many water molecules far away from the reaction center. Here, we investigated the self-ionization of water at room temperature by rare-event ab initio molecular dynamics and obtained autoionization rates and activation energies in good agreement with experiments. Based on the analysis of thousands of molecular trajectories, we identified a couple of local order parameters and show that if a bond stretch occurs when all these parameters are around their ideal range, the chance for the first dissociation step (double-proton jump) increases from [Formula: see text] to 0.4. Understanding these initiation triggers might ultimately allow the steering of chemical reactions.

Improving the mesoscopic modeling of DNA denaturation dynamics.

Although previously developed mesoscopic DNA models have successfully reproduced thermodynamic denaturation data, recent studies show that these overestimate the rate of denaturation by orders of magnitude. Using adapted Peyrard-Bishop-Dauxois (PBD) models, we have calculated the denaturation rates of several DNA hairpins and made comparison with experimental data. We show that the addition of a barrier at the onsite potential of the PBD model gives a more accurate description of the unzipping dynamics of short DNA sequences. The new models provide a refined theoretical insight on the dynamical mechanisms of unzipping which can have implications for the understanding of transcription and replication processes. Still, this class of adapted PBD models seems to have a fundamental limitation which implies that it is not possible to get agreement with available experimental results on the dynamics of long DNA sequences and at the same time maintain the good agreement regarding its thermodynamics. The reason for this is that the denaturation rate of long DNA chains is not dramatically lowered by the additional barrier-as the base-pairs that open are more likely to remain open, facilitating the opening of the full DNA molecule. Some care has to be taken, since experimental techniques suitable to the study of denaturation rates of long sequences seem not to agree with other experimental data on short DNA sequences. Further research, both theoretical and experimental, is therefore needed to resolve these inconsistencies-which will be a starting point for new minimalistic models that are able to describe both thermodynamics and dynamics at a predictive level.

Foundations and latest advances in replica exchange transition interface sampling.

Nearly 20 years ago, transition path sampling (TPS) emerged as an alternative method to free energy based approaches for the study of rare events such as nucleation, protein folding, chemical reactions, and phase transitions. TPS effectively performs Monte Carlo simulations with relatively short molecular dynamics trajectories, with the advantage of not having to alter the actual potential energy surface nor the underlying physical dynamics. Although the TPS approach also introduced a methodology to compute reaction rates, this approach was for a long time considered theoretically attractive, providing the exact same results as extensively long molecular dynamics simulations, but still expensive for most relevant applications. With the increase of computer power and improvements in the algorithmic methodology, quantitative path sampling is finding applications in more and more areas of research. In particular, the transition interface sampling (TIS) and the replica exchange TIS (RETIS) algorithms have, in turn, improved the efficiency of quantitative path sampling significantly, while maintaining the exact nature of the approach. Also, open-source software packages are making these methods, for which implementation is not straightforward, now available for a wider group of users. In addition, a blooming development takes place regarding both applications and algorithmic refinements. Therefore, it is timely to explore the wide panorama of the new developments in this field. This is the aim of this article, which focuses on the most efficient exact path sampling approach, RETIS, as well as its recent applications, extensions, and variations.

Fast Decorrelating Monte Carlo Moves for Efficient Path Sampling.

Many relevant processes in chemistry, physics, and biology are rare events from a computational perspective as they take place beyond the accessible time scale of molecular dynamics (MD). Examples are chemical reactions, nucleation, and conformational changes of biomolecules. Path sampling is an approach to break this time scale limit via a Monte Carlo (MC) sampling of MD trajectories. Still, many trajectories are needed for accurately predicting rate constants. To improve the speed of convergence, we propose two new MC moves, stone skipping and web throwing. In these moves, trajectories are constructed via a sequence of subpaths obeying superdetailed balance. By a reweighting procedure, almost all paths can be accepted. Whereas the generation of a single trajectory becomes more expensive, the reduced correlation results in a significant speedup. For a study on DNA denaturation, the increase was found to be a factor 12.

PyRETIS: A well-done, medium-sized python library for rare events.

Transition path sampling techniques are becoming common approaches in the study of rare events at the molecular scale. More efficient methods, such as transition interface sampling (TIS) and replica exchange transition interface sampling (RETIS), allow the investigation of rare events, for example, chemical reactions and structural/morphological transitions, in a reasonable computational time. Here, we present PyRETIS, a Python library for performing TIS and RETIS simulations. PyRETIS directs molecular dynamics (MD) simulations in order to sample rare events with unbiased dynamics. PyRETIS is designed to be easily interfaced with any molecular simulation package and in the present release, it has been interfaced with GROMACS and CP2K, for classical and ab initio MD simulations, respectively. © 2017 Wiley Periodicals, Inc.

Rare event simulations reveal subtle key steps in aqueous silicate condensation.

A replica exchange transition interface sampling (RETIS) study combined with Born-Oppenheimer molecular dynamics (BOMD) is used to investigate the dynamics, thermodynamics and the mechanism of the early stages of the silicate condensation process. In this process, two silicate monomers, of which one is an anionic species, form a negatively charged five-coordinated silicate dimer. In a second stage, this dimer can fall apart again, forming the original monomers, or release a water molecule into the solution. We studied the association and dissociation reaction in the gas phase, and the dissociation and water removal step in the aqueous phase. The results on the aqueous phase dissociation suggest two possible mechanisms. The breakage of the bond between the intermediate oxygen and the five-coordinated silicon is sometimes accompanied by a proton transfer. After dissociation into silicate monomers, the anionic monomer is either the previously four-coordinated silicon or the previously five-coordinated silicon depending on whether the hydrogen transfer occurs or not. Our results show that the mechanism of proton transfer is highly predominant. Water removal simulations also show two possible mechanisms distinguished by the proton transfer reaction path. Proton transfer can occur either via a direct or via a water mediated reaction step. The calculations reveal that although both mechanisms contribute to the water removal process, the direct proton transfer is slightly favorable and occurs roughly in six out of ten occasions.

Analyzing Complex Reaction Mechanisms Using Path Sampling.

We introduce an approach to analyze collective variables (CVs) regarding their predictive power for a reaction. The method is based on already available path sampling data produced by, for instance, transition interface sampling or forward flux sampling, which are path sampling methods used for efficient computation of reaction rates. By a search in CV space, a measure of predictiveness can be optimized and, in addition, the number of CVs can be reduced using projection operations which keep this measure invariant. The approach allows testing hypotheses on the reaction mechanism but could, in principle, also be used to construct the phase-space committor surfaces without the need of additional trajectory sampling. The procedure is illustrated for a one-dimensional double-well potential, a theoretical model for an ion-transfer reaction in which the solvent structure can lower the barrier, and an ab initio molecular dynamics study of water auto-ionization. The analysis technique enhances the quantitative interpretation of path sampling data which can provide clues on how chemical reactions can be steered in desired directions.

Gluing potential energy surfaces with rare event simulations.

We develop a new method combining replica exchange transition interface sampling with two distinct potential energy surfaces. The method can be used to combine different levels of theory in a simulation of a molecular process (e.g., a chemical reaction), and it can serve as a dynamical version of QM-MM, connecting classical dynamics with Ab Initio dynamics in the time domain. This new method, which we coin QuanTIS, could be applied to use accurate but expensive density functional theory based molecular dynamics for the breaking and making of chemical bonds, while the diffusion of reactants in the solvent are treated with classical force fields. We exemplify the method by applying it to two simple model systems (an ion dissociation reaction and a classical hydrogen model), and we discuss a possible extension of the method in which classical force field parameters for chemical reactions can be optimized on the fly.

A test on reactive force fields for the study of silica dimerization reactions.

We studied silica dimerization reactions in the gas and aqueous phase by density functional theory (DFT) and reactive force fields based on two parameterizations of ReaxFF. For each method (both ReaxFF force fields and DFT), we performed constrained geometry optimizations, which were subsequently evaluated in single point energy calculations using the other two methods. Standard fitting procedures typically compare the force field energies and geometries with those from quantum mechanical data after a geometry optimization. The initial configurations for the force field optimization are usually the minimum energy structures of the ab initio database. Hence, the ab initio method dictates which structures are being examined and force field parameters are being adjusted in order to minimize the differences with the ab initio data. As a result, this approach will not exclude the possibility that the force field predicts stable geometries or low transition states which are realistically very high in energy and, therefore, never considered by the ab initio method. Our analysis reveals the existence of such unphysical geometries even at unreactive conditions where the distance between the reactants is large. To test the effect of these discrepancies, we launched molecular dynamics simulations using DFT and ReaxFF and observed spurious reactions for both ReaxFF force fields. Our results suggest that the standard procedures for parameter fitting need to be improved by a mutual comparative method.

Density Functional Theory Study on the Interactions of Metal Ions with Long Chain Deprotonated Carboxylic Acids.

In this work, interactions between carboxylate ions and calcium or sodium ions are investigated via density functional theory (DFT). Despite the ubiquitous presence of these interactions in natural and industrial chemical processes, few DFT studies on these systems exist in the literature. Special focus has been placed on determining the influence of the multibody interactions (with up to 4 carboxylates and one metal ion) on an effective pair-interaction potential, such as those used in molecular mechanics (MM). Specifically, DFT calculations are employed to quantify an effective pair-potential that implicitly includes multibody interactions to construct potential energy curves for carboxylate-metal ion pairs. The DFT calculated potential curves are compared to a widely used molecular mechanics force field (OPLS-AA). The calculations indicate that multibody effects do influence the energetic behavior of these ionic pairs and the extent of this influence is determined by a balance between (a) charge transfer from the carboxylate to the metal ions which stabilizes the complex and (b) repulsion between carboxylates, which destabilizes the complex. Additionally, the potential curves of the complexes with 1 and 2 carboxylates and one counterion have been examined to higher separation distance (20 Å) by the use of relaxed scan optimization and constrained density functional theory (CDFT). The results from the relaxed scan optimization indicate that near the equilibrium distance, the charge transfer between the metal ion and the deprotonated carboxylic acid group is significant and leads to non-negligible differences between the DFT and MM potential curves, especially for calcium. However, at longer separation distances the MM calculated interaction potential functions converge to those calculated with CDFT, effectively indicating the approximate domain of the separation distance coordinate where charge transfer between the ions is occurring.