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.

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.

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.

PyVisA: Visualization and Analysis of path sampling trajectories.

Rare event methods applied to molecular simulations are growing in popularity, accessible and customizable software solutions have thus been developed and released. One of the most recent is PyRETIS, an open Python library for performing path sampling simulations. Here, we introduce PyVisA, a postprocessing package for path sampling simulations, which includes visualization and analysis tools for interpreting path sampling outputs. PyVisA integrates PyRETIS functionalities and aims to facilitate the determination of: (a) the correlation of the order parameter with other descriptors; (b) the presence of latent variables; and (c) intermediate meta-stable states. To illustrate some of the main PyVisA features, we investigate the proton transfer reaction in a protonated water trimer simulated via a simple polarizable model (Stillinger-David).

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.

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.

Predicting the mechanism and rate of H-NS binding to AT-rich DNA.

Bacteria contain several nucleoid-associated proteins that organize their genomic DNA into the nucleoid by bending, wrapping or bridging DNA. The Histone-like Nucleoid Structuring protein H-NS found in many Gram-negative bacteria is a DNA bridging protein and can structure DNA by binding to two separate DNA duplexes or to adjacent sites on the same duplex, depending on external conditions. Several nucleotide sequences have been identified to which H-NS binds with high affinity, indicating H-NS prefers AT-rich DNA. To date, highly detailed structural information of the H-NS DNA complex remains elusive. Molecular simulation can complement experiments by modelling structures and their time evolution in atomistic detail. In this paper we report an exploration of the different binding modes of H-NS to a high affinity nucleotide sequence and an estimate of the associated rate constant. By means of molecular dynamics simulations, we identified three types of binding for H-NS to AT-rich DNA. To further sample the transitions between these binding modes, we performed Replica Exchange Transition Interface Sampling, providing predictions of the mechanism and rate constant of H-NS binding to DNA. H-NS interacts with the DNA through a conserved QGR motif, aided by a conserved arginine at position 93. The QGR motif interacts first with phosphate groups, followed by the formation of hydrogen bonds between acceptors in the DNA minor groove and the sidechains of either Q112 or R114. After R114 inserts into the minor groove, the rest of the QGR motif follows. Full insertion of the QGR motif in the minor groove is stable over several tens of nanoseconds, and involves hydrogen bonds between the bases and both backbone and sidechains of the QGR motif. The rate constant for the process of H-NS binding to AT-rich DNA resulting in full insertion of the QGR motif is in the order of 10(6) M-1s-1, which is rate limiting compared to the non-specific association of H-NS to the DNA backbone at a rate of 10(8) M-1s-1.

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.

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.

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.

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.

Molecular dynamics approach to locally resolve elastic constants in nanocomposites and thin films: mechanical description of solid-soft matter interphases via Young's modulus, Poisson's ratio and shear modulus.

A molecular dynamics approach based on a small-deformation mechanical response has been extended from the evaluation of locally resolved Poisson's ratios, νj, in nanocomposites to the calculation of local Young's moduli, Ej, (with j labelling a subvolume of the studied sample). On the basis of the νj and Ej, the local values of the shear modulus, Gj, can be derived as well. The capability of the developed method to derive locally resolved elastic constants of complex (nanocomposite) systems has been tested for an atomistic model of silica and atactic polystyrene. When measuring the interphase dimension of the composite in terms of local Ej, νj and Gj elements, a surface influence exceeding three times the polymer bulk radius of gyration (Rg ≈ 1 nm in the studied 20 mer composite) is predicted while for the majority of static quantities (e.g., polymer mass density, polymer orientation relative to the nanoparticle surface, radius of gyration, end-to-end distance) interphase dimensions only slightly larger than the polymer Rg are found. Calculated local values of mechanical descriptors can be adopted as input parameters in the micromechanical modelling of multicomponent nanocomposites.

Modeling the construction of polymeric adsorbent media: effects of counter-ions on ligand immobilization and pore structure.

Molecular dynamics modeling and simulations are employed to study the effects of counter-ions on the dynamic spatial density distribution and total loading of immobilized ligands as well as on the pore structure of the resultant ion exchange chromatography adsorbent media. The results show that the porous adsorbent media formed by polymeric chain molecules involve transport mechanisms and steric resistances which cause the charged ligands and counter-ions not to follow stoichiometric distributions so that (i) a gradient in the local nonelectroneutrality occurs, (ii) non-uniform spatial density distributions of immobilized ligands and counter-ions are formed, and (iii) clouds of counter-ions outside the porous structure could be formed. The magnitude of these counter-ion effects depends on several characteristics associated with the size, structure, and valence of the counter-ions. Small spherical counter-ions with large valence encounter the least resistance to enter a porous structure and their effects result in the formation of small gradients in the local nonelectroneutrality, higher ligand loadings, and more uniform spatial density distributions of immobilized ligands, while the formation of exterior counter-ion clouds by these types of counter-ions is minimized. Counter-ions with lower valence charges, significantly larger sizes, and elongated shapes, encounter substantially greater steric resistances in entering a porous structure and lead to the formation of larger gradients in the local nonelectroneutrality, lower ligand loadings, and less uniform spatial density distributions of immobilized ligands, as well as substantial in size exterior counter-ion clouds. The effects of lower counter-ion valence on pore structure, local nonelectroneutrality, spatial ligand density distribution, and exterior counter-ion cloud formation are further enhanced by the increased size and structure of the counter-ion. Thus, the design, construction, and functionality of polymeric porous adsorbent media will significantly depend, for a given desirable ligand to be immobilized and represent the adsorption active sites, on the type of counter-ion that is used during the ligand immobilization process. Therefore, the molecular dynamics modeling and simulation approach presented in this work could contribute positively by representing an engineering science methodology to the design and construction of polymeric porous adsorbent media which could provide high intraparticle mass transfer and adsorption rates for the adsorbate biomolecules of interest which are desired to be separated by an adsorption process.

Effect of nanoroughness on highly hydrophobic and superhydrophobic coatings.

The effect of nanoroughness on contact angles and pinning is investigated experimentally and numerically for low-energy surfaces. Nanoroughness is introduced by chemical vapor deposition of tetraethoxysilane and was quantified by scanning force microscopy. Addition of a root-mean-square roughness of 2 nm on a flat surface can increase the contact angle after fluorination by a semifluorinated silane by up to 30°. On the other hand, nanoroughness can improve or impair the liquid repellency of superhydrophobic surfaces that were made from assembled raspberry particles. Molecular dynamics simulations are performed in order to gain a microscopic understanding on how the length and the surface coating density of semifluorinated silanes influence the hydrophobicity. Solid-liquid surface free energy computations reveal that the wetting behavior strongly depends on the density and alignment of the semifluorinated silane. At coating densities in the range of experimental values, some water molecules can penetrate between the semifluorinated chains, thus increasing the surface energy. Combining the experimental and numerical data exhibits that a roughness-induced increase of the contact angle competes with increased pinning caused by penetration of liquid into nanopores or between neighboring semifluorinated molecules.

Molecular modeling of polymeric adsorbent media: the effects of counter-ions on ligand immobilization and pore structure.

Molecular dynamics modeling and simulations are employed to study the immobilization of ligands on the surface of the pores of a base porous polymeric matrix. The results show the significant effects that the counter-ions have on the spatial distribution of the density of immobilized ligands as well as on the pore size and pore connectivity distributions of the porous adsorbent medium being constructed. The results for the systems studied in this work indicate that by using doubly charged counter-ions whose numbers during ligand immobilization are half to those of singly charged counter-ions, the ligand immobilization process proceeds faster and the magnitude of local nonelectroneutrality becomes smaller. More importantly, the pore structures of the adsorbent media resulting from the system using doubly charged counter-ions have porous structures that are characterized by more mid-sized pores and higher pore connectivity than the porous adsorbent structures generated by the system employing singly charged counter-ions and, furthermore, the density distribution of the immobilized ligands in the porous structures where doubly charged counter-ions are employed tends to be more uniform laterally and the ligands are surrounded by fewer counter-ions. These characteristics affected by the use of doubly charged counter-ions could provide important advantages with respect to the transport and adsorption of adsorbate biomolecules of interest. Furthermore, the results of this work indicate that the type of counter-ions being used in the ligand immobilization process could represent an additional control variable for affecting the ligand density distribution as well as the pore size and pore connectivity distributions of the porous structure of the adsorbent medium being constructed.

Exploring the effects of approach velocity on depletion force and coalescence in oil-in-water emulsions.

An emulsion is a thermodynamically unstable system consisting of at least two immiscible liquid phases, one of which is dispersed in the other in the form of droplets of varying size. Most studies on emulsions have focused on the behaviour of emulsion droplets with diameter from ∼50 µm and upwards. However, the properties of smaller droplets may be highly relevant in order to understand the behaviour of emulsions, including their performance in numerous applications within the fields of food, industry, and medical science. The relatively long life-time and small size of these droplets compared to other emulsion droplets, make them suited for optical trapping and micromanipulation technologies. Optical tweezers have previously shown potential in the study of stabilized emulsions. Here we employ optical tweezers to examine unstable oil-in-water emulsions to determine the effects of system parameters on depletion force and coalescence times.

Chemistrees: Data-Driven Identification of Reaction Pathways via Machine Learning.

We propose to analyze molecular dynamics (MD) output via a supervised machine learning (ML) algorithm, the decision tree. The approach aims to identify the predominant geometric features which correlate with trajectories that transition between two arbitrarily defined states. The data-driven algorithm aims to identify these features without the bias of human "chemical intuition". We demonstrate the method by analyzing the proton exchange reactions in formic acid solvated in small water clusters. The simulations were performed with ab initio MD combined with a method to efficiently sample the rare event, path sampling. Our ML analysis identified relevant geometric variables involved in the proton transfer reaction and how they may change as the number of solvating water molecules changes.

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.

The immuneML ecosystem for machine learning analysis of adaptive immune receptor repertoires.

Adaptive immune receptor repertoires (AIRR) are key targets for biomedical research as they record past and ongoing adaptive immune responses. The capacity of machine learning (ML) to identify complex discriminative sequence patterns renders it an ideal approach for AIRR-based diagnostic and therapeutic discovery. To date, widespread adoption of AIRR ML has been inhibited by a lack of reproducibility, transparency, and interoperability. immuneML (immuneml.uio.no) addresses these concerns by implementing each step of the AIRR ML process in an extensible, open-source software ecosystem that is based on fully specified and shareable workflows. To facilitate widespread user adoption, immuneML is available as a command-line tool and through an intuitive Galaxy web interface, and extensive documentation of workflows is provided. We demonstrate the broad applicability of immuneML by (i) reproducing a large-scale study on immune state prediction, (ii) developing, integrating, and applying a novel deep learning method for antigen specificity prediction, and (iii) showcasing streamlined interpretability-focused benchmarking of AIRR ML.