Exploring chemical reactivity of complex systems with path‐based coordinates: Role of the distance metric

Path‐based reaction coordinates constitute a valuable tool for free‐energy calculations in complex processes. When a reference path is defined by means of collective variables, a nonconstant distance metric that incorporates the nonorthonormality of these variables should be taken into account. In this work, we show that, accounting for the correct metric tensor, these kind of variables can provide iso‐hypersurfaces that coincide with the iso‐committor surfaces and that activation free energies equal the value that would be obtained if the committor function itself were used as reaction coordinate. The advantages of the incorporation of the variable metric tensor are illustrated with the analysis of the enzymatic reaction catalyzed by isochorismate‐pyruvate lyase. Hybrid QM/MM techniques are used to obtain the free energy profile and to analyze reactive trajectories initiated at the transition state. For this example, the committor histogram is peaked at 0.5 only when a variable metric tensor is incorporated in the definition of the path‐based coordinate. © 2014 Wiley Periodicals, Inc.     

Automated Search of Minimum Free‐Energy Path by Umbrella Integration

The free‐energy landscape is an important factor for understanding the conformational equilibria of chemical reactions, and many techniques have been developed to calculate the potential of the mean force. Unfortunately, these methods require a previous knowledge of the system for calculations because the results depend on the reaction coordinates. In this study, we combine the scaled hypersphere search method with the umbrella integration method to obtain the transition states on free‐energy landscapes and minimum‐free‐energy paths (MFEPs). With this approach, the MFEP connections between known and unknown equilibrium points are constructed without the prior knowledge of the free‐energy landscape. The problem of reaction coordinates can be solved by using a multidimensional, fully automated interrogation of MFEPs for acquiring the potential of mean force. The efficiency of the proposed method is demonstrated by applying it to alanine dipeptide and alanine tripeptide. © 2018 Wiley Periodicals, Inc.     

Free energy profiles for monomer capture in Grubbs- and SHOP-type olefin polymerization catalysts: a constraint ab initio molecular dynamics study

Density functional theory together with Car-Parrinello ab initio molecular dynamics simulation has been used to investigate the free energy profiles (FEP) of monomer capture in Grubbs- and SHOP-type olefin polymerization catalysts. The FEPs along the reaction coordinates at 300 K were determined directly by a point wise thermodynamic integration technique. Comparison between potential energy profile (PEP) and the FEP has been made. The results show that, for both catalysts, the PEP for the monomer ethylene uptake by the metal center is a typical Morse curve without energy barrier. However, a small barrier (1.8 kcal/mol for Grubbs catalyst and 2.4 kcal/mol for SHOP catalyst) exists on the FEP. The pi complexation energy on the FES at 300 K is higher by 10-12 kcal/mol over that on the PES. The differences between FES and PES are due to entropy contribution. Slow growth simulations on the ethylene capture process show that the ethylene attacks the metal center by an asynchronous mode. This indicates that the forming of the pi-bonding between the metal and ethylene is initiated by electrophilic attack of the metal to one of the ethylene carbons.     

Increasing the time step with mass scaling in Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations

Born‐Oppenheimer ab initio QM/MM molecular dynamics simulation with umbrella sampling is a state‐of‐the‐art approach to calculate free energy profiles of chemical reactions in complex systems. To further improve its computational efficiency, a mass‐scaling method with the increased time step in MD simulations has been explored and tested. It is found that by increasing the hydrogen mass to 10 amu, a time step of 3 fs can be employed in ab initio QM/MM MD simulations. In all our three test cases, including two solution reactions and one enzyme reaction, the resulted reaction free energy profiles with 3 fs time step and mass scaling are found to be in excellent agreement with the corresponding simulation results using 1 fs time step and the normal mass. These results indicate that for Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, the mass‐scaling method can significantly reduce its computational cost while has little effect on the calculated free energy profiles. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Effects of system net charge and electrostatic truncation on all‐atom constant pH molecular dynamics

Constant pH molecular dynamics offers a means to rigorously study the effects of solution pH on dynamical processes. Here, we address two critical questions arising from the most recent developments of the all‐atom continuous constant pH molecular dynamics (CpHMD) method: (1) What is the effect of spatial electrostatic truncation on the sampling of protonation states? (2) Is the enforcement of electrical neutrality necessary for constant pH simulations? We first examined how the generalized reaction field and force‐shifting schemes modify the electrostatic forces on the titration coordinates. Free energy simulations of model compounds were then carried out to delineate the errors in the deprotonation free energy and salt‐bridge stability due to electrostatic truncation and system net charge. Finally, CpHMD titration of a mini‐protein HP36 was used to understand the manifestation of the two types of errors in the calculated pK_a values. The major finding is that enforcing charge neutrality under all pH conditions and at all time via cotitrating ions significantly improves the accuracy of protonation‐state sampling. We suggest that such finding is also relevant for simulations with particle mesh Ewald, considering the known artifacts due to charge‐compensating background plasma. © 2014 Wiley Periodicals, Inc.     

Molecular dynamics coupled with a virtual system for effective conformational sampling

An enhanced conformational sampling method is proposed: virtual‐system coupled canonical molecular dynamics (VcMD). Although VcMD enhances sampling along a reaction coordinate, this method is free from estimation of a canonical distribution function along the reaction coordinate. This method introduces a virtual system that does not necessarily obey a physical law. To enhance sampling the virtual system couples with a molecular system to be studied. Resultant snapshots produce a canonical ensemble. This method was applied to a system consisting of two short peptides in an explicit solvent. Conventional molecular dynamics simulation, which is ten times longer than VcMD, was performed along with adaptive umbrella sampling. Free‐energy landscapes computed from the three simulations mutually converged well. The VcMD provided quicker association/dissociation motions of peptides than the conventional molecular dynamics did. The VcMD method is applicable to various complicated systems because of its methodological simplicity. © 2018 Wiley Periodicals, Inc.     

Finite-Temperature Dimer Method for Finding Saddle Points on Free Energy Surfaces

The dimer method and its variants have been shown to be efficient in finding saddle points on potential surfaces. In the dimer method, the most unstable direction is approximately obtained by minimizing the total potential energy of the dimer. Then, the force in this direction is reversed to move the dimer toward saddle points. When the finite-temperature effect is important for a high-dimensional system, one usually needs to describe the dynamics in a low-dimensional space of reaction coordinates. In this case, transition states are collected as saddle points on the free energy surface. The traditional dimer method cannot be directly employed to find saddle points on a free energy surface since the surface is not known a priori. Here, we develop a finite-temperature dimer method for searching saddle points on the free energy surface. In this method, a constrained rotation dynamics of the dimer system is used to sample dimer directions and an efficient average method is used to obtain a good approximation of the most unstable direction. This approximated direction is then used in reversing the force component and evolving the dimer toward saddle points. Our numerical results suggest that the new method is efficient in finding saddle points on free energy surfaces. © 2019 Wiley Periodicals, Inc.     

Use of Cartesian Rotation Vectors in Brownian Dynamics Algorithms: Theory and Simulation Results

Summary: We have shown that the components of Cartesian rotation vectors can be used successfully as generalized coordinates describing angular orientation in Brownian dynamics simulations of non‐spherical nanoparticles. For this particular choice of generalized coordinates, we rigorously derived the conformation‐space diffusion equations from kinetic theory for both free nanoparticles and nanoparticles interconnected by springs or holonomic constraints into polymer chains. The equivalent stochastic differential equations were used as a foundation for the Brownian dynamics algorithms. These new algorithms contain singularities only for points in the conformation‐space where both the probability density and its first coordinate derivative equal zero (weak singularities). In addition, the coordinate values after a single Brownian dynamics time step are throughout the conformation‐space equal to the old coordinate values plus the respective increments. For some parts of the conformation‐space these features represent a major improvement compared to the situation when Eulerian angles describe rotational dynamics. The presented simulation results of the equilibrium probability density for free nanoparticles are in perfect agreement with the results from kinetic theory.

Simulation of p^(eq)(Φ) for free nanoparticles.     

Polymer reversal rate calculated via locally scaled diffusion map

A recent study on the dynamics of polymer reversal inside a nanopore by Huang and Makarov [J. Chem. Phys. 128, 114903 (2008)] demonstrated that the reaction rate cannot be reproduced by projecting the dynamics onto a single empirical reaction coordinate, a result suggesting the dynamics of this system cannot be correctly described by using a single collective coordinate. To further investigate this possibility we have applied our recently developed multiscale framework, locally scaled diffusion map (LSDMap), to obtain collective reaction coordinates for this system. Using a single diffusion coordinate, we obtain a reversal rate via Kramers expression that is in good agreement with the exact rate obtained from the simulations. Our mathematically rigorous approach accounts for the local heterogeneity of molecular configuration space in constructing a diffusion map, from which collective coordinates emerge. We believe this approach can be applied in general to characterize complex macromolecular dynamics by providing an accurate definition of the collective coordinates associated with processes at different time scales.     

Enhanced sampling simulations of DNA step parameters

A novel approach for the selection of step parameters as reaction coordinates in enhanced sampling simulations of DNA is presented. The method uses three atoms per base and does not require coordinate overlays or idealized base pairs. This allowed for a highly efficient implementation of the calculation of all step parameters and their Cartesian derivatives in molecular dynamics simulations. Good correlation between the calculated and actual twist, roll, tilt, shift, and slide parameters is obtained, while the correlation with rise is modest. The method is illustrated by its application to the methylated and unmethylated 5′‐CATGTGACGTCACATG‐3′ double stranded DNA sequence. One‐dimensional umbrella simulations indicate that the flexibility of the central CG step is only marginally affected by methylation. © 2014 Wiley Periodicals, Inc.     

On‐the‐fly reconstruction of free‐energy profiles using logarithmic mean‐force dynamics

Mean‐force dynamics (MFD), which is a fictitious dynamics for a set of collective variables on a potential of mean‐force, is a powerful algorithm to efficiently explore free‐energy landscapes. Recently, we have introduced logarithmic MFD (LogMFD) (Morishita et al., Phys. Rev. E 2012, 85, 066702) which overcomes difficulties encounterd in free‐energy calculations using standard approaches such as thermodynamic integration. Here, we present a guide to implementing LogMFD calculations paying attention to the practical issues in choosing the parameters in LogMFD. A primary focus is given to the effect of the parameters on the accuracy of the reconstructed free‐energy profiles. A recipe for reducing the errors due to energy dissipation is presented. We also demonstrate that multidimensional free‐energy landscapes can be reconstructed on‐the‐fly using LogMFD, which cannot be accomplished using any other free‐energy calculation techniques. © 2013 Wiley Periodicals, Inc.     

Net charge changes in the calculation of relative ligand‐binding free energies via classical atomistic molecular dynamics simulation

The calculation of binding free energies of charged species to a target molecule is a frequently encountered problem in molecular dynamics studies of (bio‐)chemical thermodynamics. Many important endogenous receptor‐binding molecules, enzyme substrates, or drug molecules have a nonzero net charge. Absolute binding free energies, as well as binding free energies relative to another molecule with a different net charge will be affected by artifacts due to the used effective electrostatic interaction function and associated parameters (e.g., size of the computational box). In the present study, charging contributions to binding free energies of small oligoatomic ions to a series of model host cavities functionalized with different chemical groups are calculated with classical atomistic molecular dynamics simulation. Electrostatic interactions are treated using a lattice‐summation scheme or a cutoff‐truncation scheme with Barker–Watts reaction‐field correction, and the simulations are conducted in boxes of different edge lengths. It is illustrated that the charging free energies of the guest molecules in water and in the host strongly depend on the applied methodology and that neglect of correction terms for the artifacts introduced by the finite size of the simulated system and the use of an effective electrostatic interaction function considerably impairs the thermodynamic interpretation of guest‐host interactions. Application of correction terms for the various artifacts yields consistent results for the charging contribution to binding free energies and is thus a prerequisite for the valid interpretation or prediction of experimental data via molecular dynamics simulation. Analysis and correction of electrostatic artifacts according to the scheme proposed in the present study should therefore be considered an integral part of careful free‐energy calculation studies if changes in the net charge are involved. © 2013 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Systematic and statistical error in histogram-based free energy calculations

A common technique for the numerical calculation of free energies involves estimation of the probability density along a given coordinate from a set of configurations generated via simulation. The process requires discretization of one or more reaction coordinates to generate a histogram from which the continuous probability density is inferred. We show that the finite size of the intervals used to construct the histogram leads to quantifiable systematic error. The width of these intervals also determines the statistical error in the free energy, and the choice of the appropriate interval is therefore driven by the need to balance the two sources of error. We present a method for the construction of the optimal histogram for a given system, and show that the use of this technique requires little additional computational expense. We demonstrate the efficacy of the technique for a model system, and discuss how the principles governing the choice of discretization interval could be used to improve extended sampling techniques.     

CAST: A new program package for the accurate characterization of large and flexible molecular systems

The presented program package, Conformational Analysis and Search Tool (CAST) allows the accurate treatment of large and flexible (macro) molecular systems. For the determination of thermally accessible minima CAST offers the newly developed TabuSearch algorithm, but algorithms such as Monte Carlo (MC), MC with minimization, and molecular dynamics are implemented as well. For the determination of reaction paths, CAST provides the PathOpt, the Nudge Elastic band, and the umbrella sampling approach. Access to free energies is possible through the free energy perturbation approach. Along with a number of standard force fields, a newly developed symmetry‐adapted perturbation theory‐based force field is included. Semiempirical computations are possible through DFTB+ and MOPAC interfaces. For calculations based on density functional theory, a Message Passing Interface (MPI) interface to the Graphics Processing Unit (GPU)‐accelerated TeraChem program is available. The program is available on request. © 2014 Wiley Periodicals, Inc.     

Efficient calculation of two-dimensional adiabatic and free energy maps: Application to the isomerization of the C13C14 and C15N16 bonds in the retinal of bacteriorhodopsin

Accurate calculation of potential energy and free-energy profiles along reaction coordinates of biological processes such as enzymatic reactions or conformational changes is fundamental to the obtention of theoretical insight into protein function. We describe here the practical implementation of the Automatic Map Refinement Procedure (AMRP) and two-dimensional Weighted Histogram Analysis Method (WHAM) for efficient computation of adiabatic potential energy and free-energy maps, respectively. Methods for efficiently sampling configuration space with high-energy barriers and for removing hysteresis in the case of periodic reaction coordinates are presented. The application of these techniques to the isomerization of the C13C14 and C15N16 bonds in the retinal of bacteriorhodopsin is described. In dark-adapted bacteriorhodopsin (bR), the retinal moiety exists in two conformers, all-trans and (13,15)cis, with the latter making ≃67% of the population. This experimental free energy difference is reproduced here to within k_BT. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1644–1658, 1999     

Multidimensional umbrella sampling and replica‐exchange molecular dynamics simulations for structure prediction of transmembrane helix dimers

Structural information of a transmembrane (TM) helix dimer is useful in understanding molecular mechanisms of important biological phenomena such as signal transduction across the cell membrane. Here, we describe an umbrella sampling (US) scheme for predicting the structure of a TM helix dimer in implicit membrane using the interhelical crossing angle and the TM–TM relative rotation angles as the reaction coordinates. This scheme conducts an efficient conformational search on TM–TM contact interfaces, and its robustness is tested by predicting the structures of glycophorin A (GpA) and receptor tyrosine kinase EphA1 (EphA1) TM dimers. The nuclear magnetic resonance (NMR) structures of both proteins correspond to the global free‐energy minimum states in their free‐energy landscapes. In addition, using the landscape of GpA as a reference, we also examine the protocols of temperature replica‐exchange molecular dynamics (REMD) simulations for structure prediction of TM helix dimers in implicit membrane. A wide temperature range in REMD simulations, for example, 250–1000 K, is required to efficiently obtain a free‐energy landscape consistent with the US simulations. The interhelical crossing angle and the TM–TM relative rotation angles can be used as reaction coordinates in multidimensional US and be good measures for conformational sampling of REMD simulations. © 2013 Wiley Periodicals, Inc.     

Surface free energy of ladderlike polyepoxysiloxane–diamine reaction systems

We studied three kinds of ladderlike polyepoxysiloxanes, which have different side groups grafted on the ladderlike backbones. 1,3‐Bis(aminopropyl)tetramethyl disiloxane (diamine) was used as the curing agent. The reaction between ladderlike polyepoxysiloxanes and diamine was investigated by contact angle measurements and surface free energy study. Several factors such as diamine amount, reaction time, and temperature can affect the systems' surface tension (or surface free energy), which were determined by two‐liquid geometric and three‐liquid acid‐base methods. The experimental results showed that an increase in the diamine amount in the reaction systems results in an increase in the polar part of surface free energy because of electron donate characteristics of the diamine. However, because epoxy (electron acceptor) and diamine (electron donor) react fast at elevated temperatures, increasing reaction temperature decreases the polar part of the surface free energy, while increases the nonpolar part of the surface free energy. The evolution of surface free energy with time for various epoxy–diamine reaction systems at various temperatures has also been studied. It was found that it took a relatively long time (50–60 h) to reach the equilibrium state. The experimental results can be well interpreted by the epoxy–diamine reaction mechanism and van Oss–Good's Lewis acid‐base theory. © 2000 John Wiley & Sons, Inc. J Polym Sci B: Polym Phys 38: 1449–1460, 2000     

Application of free‐volume theory to self diffusion of solvents in polymers below the glass transition temperature: A review

Theories based on free‐volume concepts have been developed to characterize the self and mutual‐diffusion coefficients of low molecular weight penetrants in rubbery and glassy polymer‐solvent systems. These theories are applicable over wide ranges of temperature and concentration. The capability of free‐volume theory to describe solvent diffusion in glassy polymers is reviewed in this article. Two alternative free‐volume based approaches used to evaluate solvent self‐diffusion coefficients in glassy polymer‐solvent systems are compared in terms of their differences and applicability. The models can correlate/predict temperature and concentration dependencies of the solvent diffusion coefficient. With the appropriate accompanying thermodynamic factors they can be used to model concentration profiles in mutual diffusion processes that are Fickian such as drying of coatings. The free‐volume methodology has been found to be consistent with molecular dynamics simulations. © 2011 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys, 2011     

Local normal modes and vibrational adiabatic potentials

Local normal-mode analysis for a collinear potential energy surface generates a system of curvilinear coordinates, which are orthogonal in the mass-skewed system. The motion is locally separable in these coordinates. We compare the utility of one of the normal modes as a transversal-vibration coordinate, with the conventional choice of the direction perpendicular to the reaction coordinate in the mass-skewed system. The comparison is done for two commonly used reaction coordinates: BEBO and the steepest-descent path. Results differ for different choices of directions and reaction coordinates. Future work should concentrate on a choice of a reaction coordinate which is itself one of the normal modes.