Determination of free energy profiles by repository based adaptive umbrella sampling: bridging nonequilibrium and quasiequilibrium simulations

We propose a new adaptive sampling approach to determine free energy profiles with molecular dynamics simulations, which is called as "repository based adaptive umbrella sampling" (RBAUS). Its main idea is that a sampling repository is continuously updated based on the latest simulation data, and the accumulated knowledge and sampling history are then employed to determine whether and how to update the biasing umbrella potential for subsequent simulations. In comparison with other adaptive methods, a unique and attractive feature of the RBAUS approach is that the frequency for updating the biasing potential depends on the sampling history and is adaptively determined on the fly, which makes it possible to smoothly bridge nonequilibrium and quasiequilibrium simulations. The RBAUS method is first tested by simulations on two simple systems: a double well model system with a variety of barriers and the dissociation of a NaCl molecule in water. Its efficiency and applicability are further illustrated in ab initio quantum mechanics/molecular mechanics molecular dynamics simulations of a methyl-transfer reaction in aqueous solution.     

CPMD/GULP QM/MM interface for modeling periodic solids: Implementation and its application in the study of Y‐zeolite supported Rhn clusters

We report here the development of hybrid quantum mechanics/molecular mechanics (QM/MM) interface between the plane‐wave density functional theory based CPMD code and the empirical force‐field based GULP code for modeling periodic solids and surfaces. The hybrid QM/MM interface is based on the electrostatic coupling between QM and MM regions. The interface is designed for carrying out full relaxation of all the QM and MM atoms during geometry optimizations and molecular dynamics simulations, including the boundary atoms. Both Born–Oppenheimer and Car–Parrinello molecular dynamics schemes are enabled for the QM part during the QM/MM calculations. This interface has the advantage of parallelization of both the programs such that the QM and MM force evaluations can be carried out in parallel to model large systems. The interface program is first validated for total energy conservation and parallel scaling performance is benchmarked. Oxygen vacancy in α‐cristobalite is then studied in detail and the results are compared with a fully QM calculation and experimental data. Subsequently, we use our implementation to investigate the structure of rhodium cluster (Rh_n; n = 2 to 6) formed from Rh(C₂H₄)₂ complex adsorbed within a cavity of Y‐zeolite in a reducible atmosphere of H₂ gas. © 2016 Wiley Periodicals, Inc.     

Prospects of transition interface sampling simulations for the theoretical study of zeolite synthesis

The transition interface sampling (TIS) technique allows large free energy barriers to be overcome within reasonable simulation time, which is impossible for straightforward molecular dynamics. Still, the method does not impose an artificial driving force, but it surmounts the timescale problem by an importance sampling of true dynamical pathways. Recently, it was shown that the efficiency of TIS when calculating reaction rates is less sensitive to the choice of reaction coordinate than those of the standard free energy based techniques. This could be an important advantage in complex systems for which a good reaction coordinate is usually very difficult to find. We explain the principles of this method and discuss some of the promising applications related to zeolite formation.     

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.     

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     

Repeated‐annealing sampling combined with multicanonical algorithm for conformational sampling of bio‐molecules

A novel conformational sampling method (repeated‐annealing sampling method) is proposed to execute an efficient conformational sampling at a reasonable computational cost. In the method, a molecular dynamics simulation is done with repeating an elemental process. An elemental process consists of four subprocesses: high‐temperature run, annealing, room‐temperature run, and fast heating. The sampling is done automatically according to a temperature‐control schedule. The room‐temperature run is treated with the multicanonical algorithm, and the other subprocesses are done with the conventional molecular dynamics algorithm. The method, differing from the generalized ensemble methods recently developed, is not warrantable to give the canonical ensemble because of the nonphysical process in the annealing. However, we observed that the slower the annealing and the longer the high‐temperature run, the closer the sampled conformations to those of the canonical ensemble. A test was performed with tri‐N‐acetyl‐D ‐glucosamine in vacuo, and the results were compared with those from the conventional multicanonical simulation. Not only the reweighted canonical distribution function but also the energy landscape were in good agreement with those from the conventional multicanonical simulation. The potential of mean force also showed a fairly good agreement with that from the conventional multicanonical simulation in the room‐temperature region. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1098–1106, 2001     

An improved replica-exchange sampling method: temperature intervals with global energy reassignment

In a molecular dynamics (MD) simulation, representative sampling over the entire phase space is desired to obtain an accurate canonical distribution at a given temperature. For large molecules, such as proteins, this is problematic because systems tend to become trapped in local energy minima. The extensively used replica-exchange molecular dynamics (REMD) simulation technique overcomes this kinetic-trapping problem by allowing Boltzmann-weighted configuration exchange processes to occur between numerous thermally adjacent and compositionally identical simulations that are thermostated at sequentially higher temperatures. While the REMD method provides much better sampling than conventional MD, there are two substantial difficulties that are inherent in its application: (1) the large number of replicas that must be used to span a designated temperature range and (2) the subsequent long time required for configurations sampled at high temperatures to exchange down for potential inclusion within the low-temperature ensemble of interest. In this work, a new method based on temperature intervals with global energy reassignment (TIGER) is presented that overcomes both of these problems. A TIGER simulation is conducted as a series of short heating-sampling-quenching cycles. At the end of each cycle, the potential energies of all replicas are simultaneously compared at the same temperature using a Metropolis sampling method and then globally reassigned to the designated temperature levels. TIGER is compared with regular MD and REMD methods for the alanine dipeptide in water. The results indicate that TIGER increases sampling efficiency while substantially reducing the number of central processing units required for a comparable conventional REMD simulation.     

Numerical examination of the extended phase‐space volume‐preserving integrator by the Nosé‐Hoover molecular dynamics equations

This article illustrates practical applications to molecular dynamics simulations of the recently developed numerical integrators [Phys Rev E 2006, 73, 026703] for ordinary differential equations. This method consists of extending any set of ordinary differential equations in order to define a time invariant function, and then use the techniques of divergence‐free solvable decomposition and symmetric composition to obtain volume‐preserving integrators in the extended phase space. Here, we have developed the technique by constructing multiple extended‐variable formalism in order to enhance the handling in actual simulation, and by constituting higher order integrators to obtain further accuracies. Using these integrators, we perform constant temperature molecular dynamics simulations of liquid water, liquid argon and peptide in liquid water droplet. The temperature control is obtained through an extended version of the Nosé‐Hoover equations. Analyzing the effects of the simulation conditions including time step length, initial values, boundary conditions, and equation parameters, we investigate local accuracy, global accuracy, computational cost, and sensitivity along with the sampling validity. According to the results of these simulations, we show that the volume‐preserving integrators developed by the current method are more effective than traditional integrators that lack the volume‐preserving property. © 2008 Wiley Periodicals, Inc. J Comput Chem, 2009     

Sampling free energy surfaces as slices by combining umbrella sampling and metadynamics

Metadynamics (MTD) is a very powerful technique to sample high‐dimensional free energy landscapes, and due to its self‐guiding property, the method has been successful in studying complex reactions and conformational changes. MTD sampling is based on filling the free energy basins by biasing potentials and thus for cases with flat, broad, and unbound free energy wells, the computational time to sample them becomes very large. To alleviate this problem, we combine the standard Umbrella Sampling (US) technique with MTD to sample orthogonal collective variables (CVs) in a simultaneous way. Within this scheme, we construct the equilibrium distribution of CVs from biased distributions obtained from independent MTD simulations with umbrella potentials. Reweighting is carried out by a procedure that combines US reweighting and Tiwary–Parrinello MTD reweighting within the Weighted Histogram Analysis Method (WHAM). The approach is ideal for a controlled sampling of a CV in a MTD simulation, making it computationally efficient in sampling flat, broad, and unbound free energy surfaces. This technique also allows for a distributed sampling of a high‐dimensional free energy surface, further increasing the computational efficiency in sampling. We demonstrate the application of this technique in sampling high‐dimensional surface for various chemical reactions using ab initio and QM/MM hybrid molecular dynamics simulations. Further, to carry out MTD bias reweighting for computing forward reaction barriers in ab initio or QM/MM simulations, we propose a computationally affordable approach that does not require recrossing trajectories. © 2016 Wiley Periodicals, Inc.     

Unbinding of the streptavidin-biotin complex by atomic force microscopy: a hybrid simulation study

A hybrid molecular simulation technique, which combines molecular dynamics and continuum mechanics, was used to study the single-molecule unbinding force of a streptavidin-biotin complex. The hybrid method enables atomistic simulations of unbinding events at the millisecond time scale of atomic force microscopy (AFM) experiments. The logarithmic relationship between the unbinding force of the streptavidin-biotin complex and the loading rate (the product of cantilever spring constant and pulling velocity) in AFM experiments was confirmed by hybrid simulations. The unbinding forces, cantilever and tip positions, locations of energy barriers, and unbinding pathway were analyzed. Hybrid simulation results from this work not only interpret unbinding AFM experiments but also provide detailed molecular information not available in AFM experiments.     

Enabling grand‐canonical Monte Carlo: Extending the flexibility of GROMACS through the GromPy python interface module

We report on a python interface to the GROMACS molecular simulation package, GromPy (available at https://github.com/GromPy ). This application programming interface (API) uses the ctypes python module that allows function calls to shared libraries, for example, written in C. To the best of our knowledge, this is the first reported interface to the GROMACS library that uses direct library calls. GromPy can be used for extending the current GROMACS simulation and analysis modes. In this work, we demonstrate that the interface enables hybrid Monte‐Carlo/molecular dynamics (MD) simulations in the grand‐canonical ensemble, a simulation mode that is currently not implemented in GROMACS. For this application, the interplay between GromPy and GROMACS requires only minor modifications of the GROMACS source code, not affecting the operation, efficiency, and performance of the GROMACS applications. We validate the grand‐canonical application against MD in the canonical ensemble by comparison of equations of state. The results of the grand‐canonical simulations are in complete agreement with MD in the canonical ensemble. The python overhead of the grand‐canonical scheme is only minimal. © 2012 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.     

Nucleophilic substitution reactions at liquid/liquid interfaces: molecular dynamics simulation of a model S(N)1 dissociation reaction at the water/carbon tetrachloride interface

The ionic dissociation step of the nucleophilic substitution reaction t-BuCl --> t-Bu(+) + Cl(-) is studied at the water/carbon tetrachloride interface using molecular dynamics computer simulations. The empirical valence bond approach is used to couple two diabatic states, covalent and ionic, in the electronically adiabatic limit. The umbrella sampling technique is used to calculate the potential of mean force along the reaction coordinate (defined as the t-Bu to Cl distance) at several interface regions of varying distances from the Gibbs dividing surface. We find a significant increase of the ionic dissociation barrier height and of the reaction free energy at the interface relative to bulk water. This is shown to be due to the reduced polarity of the interface which causes a destabilization of the pure ionic state. However, deformation to the neat interface structure in the form of water protrusions into the organic phase may provide partial stabilization of the ionic species. The importance of these structural effects is examined by repeating the calculations with an artificially smooth interface. The destabilization of the ionic state at the interface also manifests itself with a rapid (picosecond time scale) recombination dynamics of the ions to form the parent molecule followed by a slow vibrational relaxation.     

Modeling the Tautomeric Equilibrium and Absorption Spectrum of 4,5-Dimethyl-2-(2'-hydroxyphenyl)imidazole

Russian Journal of Physical Chemistry A - The potential of a new interface for molecular dynamics simulations with potentials of combined quantum mechanics/molecular mechanics is demonstrated...     

"Up" versus "down" alignment and hydration structures of solutes at the air/water interface revealed by heterodyne-detected electronic sum frequency generation with classical molecular dynamics simulation

We unambiguously demonstrate the "up" versus "down" alignment of a pair of prototypical solute molecules adsorbed at the air/water interface for the first time using heterodyne-detected electronic sum frequency generation spectroscopy. This molecular alignment is also reproduced by classical molecular dynamics (MD) simulation theoretically. Furthermore, the MD simulation indicates distinctly different interface-specific hydration structures around the two solute molecules, which dictate the molecular alignment at the interface. It is concluded that the hydrophilicity difference between the terminal functional groups of the solute governs the molecular orientation and surrounding hydration structures at the interface.     

Improved sampling for simulations of interfacial membrane proteins: application of generalized shadow hybrid Monte Carlo to a peptide toxin/bilayer system

The computational costs associated with performing molecular dynamics (MD) simulations are still somewhat prohibitive and therefore limit the time and length scales that can be currently achieved. One approach to overcoming the limited size and duration of a simulation is to reduce the amount of detail when representing a system of interest, generally termed "coarse-graining". An alternative approach is via more efficient sampling methods that offer an enhanced search of a complex multidimensional energy landscape. One could also combine enhanced sampling methods with a coarse-grained (CG) force field. Here, we apply generalized shadow hybrid Monte Carlo (GSHMC), a recently proposed simulation protocol, to a biomolecular system of moderate size and show that GSHMC offers improved sampling compared to standard MD simulation. Our test system is a CG representation of a small peptide toxin interacting with a phospholipid bilayer. Specifically, we show that GSHMC allows for a quicker localization of the toxin to its equilibrium location of interaction at the headgroup/water interface of the bilayer. GSHMC therefore potentially allows for future exploration of larger and more complex systems over longer periods, which would otherwise be impractical to perform using conventional simulation methodology.     

Quantum mechanics/molecular mechanics strategies for docking pose refinement: distinguishing between binders and decoys in cytochrome C peroxidase

We investigate the effect of systematically applying molecular dynamics (MD) and quantum mechanics/molecular mechanics (QM/MM) to docked poses in an attempt to improve the correspondence between theoretical prediction and experimental observation. The proposed scheme involves running a short time scale MD simulation on a docked ligand pose (and any known structurally important crystal structure waters in the active site), followed by QM/MM minimization. Both of these steps are relatively fast for moderately sized ligands; longer time scale MD involving the protein is not found to improve the results. The final binding energy is given in terms of the QM/MM total energy, a van der Waals correction, and a term to account for desolvation effects. This methodology is first tested with a trypsin inhibitor, for which we establish the importance of running MD before reoptimizing with QM/MM. The method is then applied to cytochrome c peroxidase using a set of binders and decoys. In this example, the proposed methodology affords much better discrimination between binders and decoys than the traditional docking approach used. For both systems presented, application of this protocol results in a significantly better energetic ranking and a smaller root mean squared deviation from known crystallographic ligand poses. This work highlights the importance of including polarization effects through QM/MM and of sampling with MD to refine a set of initial docked poses.     

Ab initio simulations of two‐dimensional electronic spectra: The SOS//QM/MM approach

Two‐dimensional electronic spectroscopy (2DES) is a cutting‐edge technique for investigating with high temporal resolution energy transfer, structure, and dynamics in a wide range of systems in physical chemistry, energy sciences, biophysics, and biocatalysis. However, the interpretation of 2DES is challenging and requires computational modeling. This perspective provides a roadmap for the development of computational tools that could be routinely applied to simulate 2DES spectra of multichromophoric systems active in the UV region (2DUV) using state‐of‐the‐art ab initio electronic structure methods within a quatum mechanics/molecular mechanics (QM/MM) scheme and the sum‐over‐states (SOS) approach (here called SOS//QM/MM). Multiconfigurational and multireference perturbative methods, such as the complete active space self‐consistent field and second‐order multireference perturbation theory (CASPT2) techniques, can be applied to reliably calculate the electronic properties of multichromophoric systems. Hybrid QM/MM method and molecular dynamics techniques can be used to assess environmental and conformational effects, respectively, that shape the 2D electronic spectra. DNA and proteins are important biological targets containing UV chromophores. We report ab initio simulation of 2DUV spectra of a cyclic tetrapeptide containing two interacting aromatic side chains, a model system for the study of protein structure and dynamics by means of 2DUV spectroscopy. © 2013 Wiley Periodicals, Inc.     

On the Binding of Congo Red to Amyloid Fibrils

Amyloids are characterized by their capacity to bind Congo red (CR), one of the most used amyloid‐specific dyes. The structural features of CR binding were unknown for years, mainly because of the lack of amyloid structures solved at high resolution. In the last few years, solid‐state NMR spectroscopy enabled the determination of the structural features of amyloids, such as the HET‐s prion forming domain (HET‐s PFD), which also has recently been used to determine the amyloid–CR interface at atomic resolution. Herein, we combine spectroscopic data with molecular docking, molecular dynamics, and excitonic quantum/molecular mechanics calculations to examine and rationalize CR binding to amyloids. In contrast to a previous assumption on the binding mode, our results suggest that CR binding to the HET‐s PFD involves a cooperative process entailing the formation of a complex with 1:1 stoichiometry. This provides a molecular basis to explain the bathochromic shift in the maximal absorbance wavelength when CR is bound to amyloids.