自动化势能面搜索方法新进展

本文简要介绍了势能面搜索的两个主要问题,即过渡态搜索和全局势能面搜索,并针对这两个问题介绍了本研究组最近发展的3种方法,即限制最小化双子算法、偏置势函数辅助限制性最小化双子算法以及势能面随机行走方法．这3种方法均只需计算一阶梯度,能够用于快速自动化的搜索过渡态以及势能面．通过几个典型例子分别说明了3种方法的特点及优势．     

Randomized tree construction algorithm to explore energy landscapes

In this work, a new method for exploring conformational energy landscapes is described. The method, called transition-rapidly exploring random tree (T-RRT), combines ideas from statistical physics and robot path planning algorithms. A search tree is constructed on the conformational space starting from a given state. The tree expansion is driven by a double strategy: on the one hand, it is naturally biased toward yet unexplored regions of the space; on the other, a Monte Carlo-like transition test guides the expansion toward energetically favorable regions. The balance between these two strategies is automatically achieved due to a self-tuning mechanism. The method is able to efficiently find both energy minima and transition paths between them. As a proof of concept, the method is applied to two academic benchmarks and the alanine dipeptide.     

MC(JBW): Simple but smart Monte Carlo algorithm for free energy simulations of multiconformational molecules

Many of the most common molecular simulation methods, including Monte Carlo (MC) and molecular or stochastic dynamics (MD or SD), have significant difficulties in sampling the space of molecular potential energy surfaces characterized by multiple conformational minima and significant energy barriers. In such cases improved sampling can be obtained by special techniques that lower such barriers or somehow direct search steps toward different low energy regions of space. We recently described a hybrid MC/SD algorithm [MC(JBW)/SD] incorporating such a technique that directed MC moves of selected torsion and bond angles toward known low energy regions of conformational space. Exploration of other degrees of freedom was left to the SD part of the hybrid algorithm. In the work described here, we develop a related but simpler simulation algorithm that uses only MC to sample all degrees of freedom (e.g., stretch, bend, and torsion). We term this algorithm MC(JBW). Using simulations on various model potential energy surfaces and on simple molecular systems (n-pentane, n-butane, and cyclohexane), MC(JBW) is shown to generate ensembles of states that are indistinguishable from the canonical ensembles generated by classical Metropolis MC in the limit of very long simulations. We further demonstrate the utility of MC(JBW) by evaluating the room temperature free energy differences between conformers of various substituted cyclohexanes and the larger ring hydrocarbons cycloheptane, cyclooctane, cyclononane, and cyclodecane. The results compare favorably with available experimental data and results from previously reported MC(JBW)/SD conformational free energy calculations. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1736–1745, 1998     

A flexible nudged elastic band program for optimization of minimum energy pathways using ab initio electronic structure methods

A driver program for carrying out nudged elastic band optimizations of minimum energy reaction pathways is described. This approach allows for the determination of minimum energy pathways using only energies and gradient information. The driver code has been interfaced with the GAUSSIAN 98 program. Applications to two isomerization reactions and to a cluster model for H(2) desorption from the Si(100)-2 x 1 surface are presented.     

Binding affinity prediction with different force fields: examination of the linear interaction energy method

A systematic study of the linear interaction energy (LIE) method and the possible dependence of its parameterization on the force field and system (receptor binding site) is reported. We have calculated the binding free energy for nine different ligands in complex with P450cam using three different force fields (Amber95, Gromos87, and OPLS-AA). The results from these LIE calculations using our earlier parameterization give relative free energies of binding that agree remarkably well with the experimental data. However, the absolute energies are too positive for all three force fields, and it is clear that an additional constant term (gamma) is required in this case. Out of five examined LIE models, the same one emerges as the best for all three force fields, and this, in fact, corresponds to our earlier one apart from the addition of the constant gamma, which is almost identical for the three force fields. Thus, the present free energy calculations clearly indicate that the coefficients of the LIE method are independent of the force field used. Their relation to solvation free energies is also demonstrated. The only free parameter of the best model is gamma, which is found to depend on the hydrophobicity of the binding site. We also attempt to quantify the binding site hydrophobicity of four different proteins which shows that the ordering of gamma's for these sites reflects the fraction of hydrophobic surface area.     

Bifurcation of reaction pathways: the set of valley ridge inflection points of a simple three-dimensional potential energy surface

This paper serves for the better understanding of the branching phenomenon of reaction paths of potential energy hypersurfaces in more than two dimensions. We apply the recently proposed reduced gradient following (RGF) method for the analysis of potential energy hypersurfaces having valley-ridge inflection (VRI) points. VRI points indicate the region of possible reaction path bifurcation. The relation between RGF and the so-called global Newton search for stationary points (Branin method) is shown. Using a 3D polynomial test surface, a whole 1D manifold of VRI points is obtained. Its relation to RGF curves, steepest descent and gradient extremals is discussed as well as the relation of the VRI manifold to bifurcation points of these curves. Received: 8 July 1998 / Accepted: 24 August 1998 / Published online: 23 November 1998     

Effective polarization energy of the naphthalene molecular crystal: a study on the polarizable force field

polarization energy of the localized charge in organic solids consists of electronic polarization energy, permanent electrostatic interactions, and inter/intra molecular relaxation energies. The effective electronic polarization energies for an electron/hole carrier were successfully estimated by AMOEBA polarizable force field in naphthalene molecular crystals. Both electronic polarization energy and permanent electrostatic interaction were in agreement with the preview experimental values. In addition, the influence of the multipoles from different distributed mutipole analysis (DMA) fitting options on the electrostatic interactions are discussed in this paper. We found that the multipoles obtained from Gauss-Hermite quadrature without diffuse function or grid-based quadrature with 0.325 Å H atomic radius will give reasonable electronic polarization energies and permanent interactions for electron and hole carriers.     

Conformational transition pathway in the allosteric process of calcium-induced recoverin: molecular dynamics simulations

Recoverin is an important neuronal calcium sensor (NCS) protein, which have been implicated in a wide range of Ca(2+) signaling events in neurons and photoreceptors. To characterize the conformational transition of recoverin from the myristoyl sequestered state to the extrusion state, a series of conventional molecular dynamics (CMD) and targeted molecular dynamics (TMD) simulations were performed. The 36.8 ns long CMD and TMD simulations on recoverin revealed a reliably conformational transition pathway, which can be viewed as a sequential two-stage process. A very important mechanistic conclusion from the present TMD simulations is that the hydrophobic and hydrophilic interactions modulate the allostery cooperatively in the conformational transition pathway. In the first stage, three salt-bridges broken between Lys-84 and Gly-124, between Lys-5 and Glu-103 and between Gly-16 and Lys-97 are major components to destabilize the structure of state T and trigger the swivel of the N- and C-terminal domains. In the second stage, the rupture of H-bond Phe-56-O(...)H(O)-Thr-21 leads to the two helices of EF-1 apart from each other, destabilizing the hydrophobic interactions of the myristoyl group with its environment, whereas the making of H-bond Leu-108-O(...)H(O)-Ser-72 helps the interfacial domain maintain its structural integrity during the course of the myristoyl extrusion. The molecular dynamics simulations results are beneficial to understanding the mechanism of how and why NCS proteins make progress in the photo-signal transduction processes. Further experimental and theoretical studies are still very desirable.     

Conformational free energy surface of the linear DPDPE peptide: cost of pre-organization for disulfide bond formation

We have calculated the free energy differences between four conformers of the linear form of the opioid pentapeptide DPDPE in aqueous solution. The conformers are Cyc, representing the structure adopted by the linear peptide prior to disulfide bond formation, β _C and β _E , two slightly different β-turns previously identified in unconstrained molecular dynamics simulations, and Ext, an extended structure. Our simulations indicate that β _E is the most stable of the studied conformers of linear DPDPE in aqueous solution, with β _C , Cyc and Ext having free energies higher by 2.3, 6.3, and 28.2 kcal/mol, respectively. The free energy differences of 4.0 kcal/mol between β _C and Cyc, and 6.3 kcal/mol between β _E and Cyc, reflect the cost of pre-organizing the linear peptide into a conformation conducive for disulfide bond formation. Such a conformational change is a pre-requisite for the chemical reaction of S–S bond formation to proceed. The relatively low population of the cyclic-like structure agrees qualitatively with observed lower potency and different receptor specificity of the linear form relative to the cyclic peptide, and with previous unconstrained simulation results. Free energy component analysis indicates that the moderate stability difference of 4.0–6.3 kcal/mol between the β-turns and the cyclic-like structure results from cancellation of two large opposing effects. In accord with intuition, the relaxed β-turns have conformational strain 43–45 kcal/mol lower than the Cyc structure. However, the cyclic-like conformer interacts with water about 39 kcal/mol strongly than the open β-turns. Our simulations are the first application of the recently developed multidimensional conformational free energy thermodynamic integration (CFTI) protocol to a solvated system, with fast convergence of the free energy obtained by fixing all flexible dihedrals. Additionally, the availability of the CFTI multidimensional free energy gradient leads to a new decomposition scheme, giving the contribution of each fixed dihedral to the overall free energy change and providing additional insight into the microscopic mechanisms of the studied processes. Received: 20 April 1998 / Accepted: 9 September 1998 / Published online: 7 December 1998     

The flying ice cube: Velocity rescaling in molecular dynamics leads to violation of energy equipartition

This article describes an unexpected phenomenon encountered during MD simulations: velocity rescaling using standard protocols can systematically change the proportion of total kinetic energy (KE) found in motions associated with the various degrees of freedom. Under these conditions, the simulation violates the principle of equipartition of energy, which requires a mean kinetic energy of RT/2 in each degree of freedom. A particularly pathological form of this problem occurs if one does not periodically remove the net translation of (and rotation about) the center of mass. In this case, almost all of the kinetic energy is converted into these two kinds of motion, producing a system with almost no kinetic energy associated with the internal degrees of freedom. We call this phenomenon “the flying ice cube.” We present a mathematical analysis of a simple diatomic system with two degrees of freedom, to document the origin of the problem. We then present examples from three kinds of MD simulations, one being an in vacuo simulation on a diatomic system, one involving a low resolution model of DNA in vacuo, and the third using a traditional all-atom DNA model with full solvation, periodic boundary conditions, and the particle mesh Ewald method for treating long-range electrostatics. Finally, we discuss methods for avoiding the problem. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 726–740, 1998     

Binding of active site directed ligands to phospholipase A2: Implications on the molecular constraints and catalytic mechanism

Molecular constraints for the localization of active site directed ligands (competitive inhibitors and substrates) in the active site of phospholipase A2 (PLA2) are characterized. Structure activity relationships with known inhibitors suggest that the head group interactions dominate the selectivity as well as a substantial part of the affinity. Theab initio fitting of the amide ligands in the active site was carried out to characterize the head group interactions. Based on a systematic coordinate space search, formamide is docked with known experimental constraints such as coordination of the carbonyl group to Ca²⁺ and hydrogen bond between amide nitrogen and ND1 of His48. An optimal position for a bound water molecule is identified and its significance for the catalytic mechanism is postulated. Unlike the traditional “pseudo-triad” mechanism, the “Ca-coordinated-oxyanion” mechanism proposed here invokes activation of the catalytic water to form the oxyanion in the coordination sphere of calcium. As it attacks the carbonyl carbon of the ester, a near-tetrahedral intermediate is formed. As the second proton of the catalytic water is abstracted by the ester oxygen, its reorientation and simultaneous cleavage form hydrogen bond with ND1 of His48. In this mechanism of esterolysis, a catalytic role for the water co-ordinated to Ca²⁺ is recognised.     

The READY program: Building a global potential energy surface and reactive dynamic simulations for the hydrogen combustion

READY (REActive DYnamics) is a program for studying reactive dynamic systems using a global potential energy surface (PES) built from previously existing PESs corresponding to each of the most important elementary reactions present in the system. We present an application to the combustion dynamics of a mixture of hydrogen and oxygen using accurate PESs for all the systems involving up to four oxygen and hydrogen atoms. Results at the temperature of 4000 K and pressure of 2 atm are presented and compared with model based on rate constants. Drawbacks and advantages of this approach are discussed and future directions of research are pointed out. © 2014 Wiley Periodicals, Inc.     

The effect of confinement on the electronic energy and polarizability of a hydrogen molecular ion

The electronic energy and the polarizability of a confined hydrogen molecular ion in the ground state and the first excited state, for cavities of different volumes, are calculated using the variational method. In the treatment adopted an alternative molecular wave function is introduced with only one variational parameter and based on wave functions used for confined atoms. © 2016 Wiley Periodicals, Inc.     

Sampling enhancement for the quantum mechanical potential based molecular dynamics simulations: a general algorithm and its extension for free energy calculation on rugged energy surface

An approach is developed in the replica exchange framework to enhance conformational sampling for the quantum mechanical (QM) potential based molecular dynamics simulations. Importantly, with our enhanced sampling treatment, a decent convergence for electronic structure self-consistent-field calculation is robustly guaranteed, which is made possible in our replica exchange design by avoiding direct structure exchanges between the QM-related replicas and the activated (scaled by low scaling parameters or treated with high "effective temperatures") molecular mechanical (MM) replicas. Although the present approach represents one of the early efforts in the enhanced sampling developments specifically for quantum mechanical potentials, the QM-based simulations treated with the present technique can possess the similar sampling efficiency to the MM based simulations treated with the Hamiltonian replica exchange method (HREM). In the present paper, by combining this sampling method with one of our recent developments (the dual-topology alchemical HREM approach), we also introduce a method for the sampling enhanced QM-based free energy calculations.     

Efficient evaluation of the accuracy of molecular quantum dynamics on an approximate analytical or interpolated ab initio potential energy surface

Ab initio electronic structure methods have reached a satisfactory accuracy for the calculation of static properties, but remain too expensive for quantum dynamical calculations. Recently, an efficient semiclassical method was proposed to evaluate the accuracy of quantum dynamics on an approximate potential without having to perform the expensive quantum dynamics on the accurate potential. Here, this method is applied for the first time to evaluate the accuracy of quantum dynamics on an approximate analytical or interpolated potential in comparison to the quantum dynamics on an accurate potential obtained by an ab initio electronic structure method. Specifically, the vibrational dynamics of H₂ on a Morse potential is compared with that on the full CI potential, and the photodissociation dynamics of CO₂ on a LEPS potential with that on the excited ¹Π surface computed at the EOM‐CCSD/aug‐cc‐pVDZ level of theory. Finally, the effect of discretization of a potential energy surface on the quantum dynamics is evaluated. © 2010 Wiley Periodicals, Inc. Int J Quantum Chem 110:2426–2435, 2010     

Thermochemistry data from kinetics results: A test of quality of the potential energy surface

A straightforward test of the quality of potential energy surfaces in polyatomic systems is proposed based on obtaining thermochemical properties, in this case, standard enthalpies of formation, from kinetics results. Agreement with experiment lends confidence to the PESs, and the disagreement serves to help improve them.     

Localization and quantification of hydrophobicity: The molecular free energy density (MolFESD) concept and its application to sweetness recognition

A method for the localization, the quantification, and the analysis of hydrophobicity of a molecule or a molecular fragment is presented. It is shown that the free energy of solvation for a molecule or the transfer free energy from one solvent to another can be represented by a surface integral of a scalar quantity, the molecular free energy surface density (MolFESD), over the solvent accessible surface of that molecule. This MolFESD concept is based on a model approach where the solvent molecules are considered to be small in comparison to the solute molecule, and the solvent can be represented by a continuous medium with a given dielectric constant. The transfer energy surface density for a 1-octanol/water system is empirically determined employing a set of atomic increment contributions and distance dependent membership functions measuring the contribution of the increments to the surface value of the MolFESD. The MolFESD concept can be well used for the quantification of the purely hydrophobic contribution to the binding constants of molecule-receptor complexes. This is demonstrated with the sweeteners sucrose and sucralose and various halogen derivatives. Therein the relative sweetness, which is assumed to be proportional to the binding constant, nicely correlates to the surface integral over the positive, hydrophobic part of the MolFESD, indicating that the sweetness receptor can be characterized by a highly flexible hydrophobic pocket instead of a localized binding site.