The Lagrange multiplier method for manipulating geometries. Implementation and applications in molecular mechanics

It is shown that a Lagrange multiplier method to constrain one or several internal coordinates, or averages and combinations of these, is easily implemented in a molecular mechanics computer program that uses Newton–Raphson (NR ) minimization. Results are given for constraints on nonbonded distances and torsion angles. When a potential energy surface is to be explored, it is much better to constrain the average of three torsion angles around a bond than to constrain a single torsion angle. Certain conversions can only be achieved when averages of torsion angles around different bonds are constrained. Combinations of constraints have been applied to evaluate differences between calculated and observed geometries and to obtain transition states for relatively large molecules from results for smaller molecules at relatively low costs. The efficieny of the combination of the Lagrange multiplier method and NR minimization in terms of computing time can be rated as good.     

B-spline method for energy minimization in grid-based molecular mechanics calculations

A method is described for molecular mechanics calculations based on a cubic B-spline approximation of the potential energy. This method is useful when parts of the system are allowed to remain fixed in position so that a potential energy grid can be precalculated and used to approximate the interaction energy between parts of a molecule or between molecules. We adapted and modified the conventional B-spline method to provide an approximation of the Empirical Conformational Energy Program for Peptides (ECEPP) potential energy function. The advantage of the B-spline method over simpler approximations is that the resulting B-spline function is C2 continuous, which allows minimization of the potential energy by any local minimization algorithm. The standard B-spline method provides a good approximation of the electrostatic energy; but in order to reproduce the Lennard–Jones and hydrogen-bonding functional forms accurately, it was necessary to modify the standard B-spline method. This modification of the B-spline method can also be used to improve the accuracy of trilinear interpolation for simulations that do not require continuous derivatives. As an example, we apply the B-spline method to rigid-body docking energy calculations using the ECEPP potential energy function. Energies are calculated for the complex of Phe-Pro-Arg with thrombin. For this system, we compare the performance of the B-spline method to that of the standard pairwise summation in terms of speed, accuracy, and overhead costs for a variety of grid spacings. In our rigid-body docking calculations, the B-spline method provided an accurate approximation of the total energy of the system, and it resulted in an 180-fold reduction in the time required for a single energy and gradient calculation for this system. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 71–85, 1998     

The molecular mechanics of valinomycin. I. Energy minimization calculations on the uncomplexed ionophore

We have used energy minimization calculations to study a number of conformations of uncomplexed valinomycin. In certain cases, x-ray diffraction atomic coordinates were used directly as input coordinates, while in other cases, conformations were found by altering the x-ray coordinates prior to minimization. Five calculated conformations are reported along with their relative energies. The conformation found theoretically to be the most stable is in agreement with earlier, cruder calculations, but does not correspond to the predominant conformation observed in nonpolar solvents. A possible rationale is presented.     

An efficient newton-like method for molecular mechanics energy minimization of large molecules

Techniques from numerical analysis and crystallographic refinement have been combined to produce a variant of the Truncated Newton nonlinear optimization procedure. The new algorithm shows particular promise for potential energy minimization of large molecular systems. Usual implementations of Newton's method require storage space proportional to the number of atoms squared (i.e., O(N²)) and computer time of O(N³). Our suggested implementation of the Truncated Newton technique requires storage of less than O(N^1.5) and CPU time of less than O(N²) for structures containing several hundred to a few thousand atoms. The algorithm exhibits quadratic convergence near the minimum and is also very tolerant of poor initial structures. A comparison with existing optimization procedures is detailed for cyclohexane, arachidonic acid, and the small protein crambin. In particular, a structure for crambin (662 atoms) has been refined to an RMS gradient of 3.6 × 10^?6 kcal/mol/Å per atom on the MM2 potential energy surface. Several suggestions are made which may lead to further improvement of the new method.     

Determination of the major solution conformation of tyrocidine A,using molecular mechanics energy minimization and NMR-derived distance and torsion angle constraints

Molecular mechanics energy calculations coupled with nuclear magnetic resonance-determined distance and torsion angle constraints have been used to determine the three-dimensional structure of tyrocidine A, a cyclic decapeptide which exists largely as a single conformation in solution. Two open-chain polyalanine models were used to represent separate halves of the peptide backbone and a combinatorial method of searching conformation space used to generate candidate structures consistent with experimental distance constraints. These structures were energy-minimized using the AMBER molecular mechanics forcefield and the resulting conformations classified by factor analysis of their Cartesian coordinates. Representative low-energy conformers of the two halves of the backbone were fused together and two candidate conformations of the completed backbone refined by further minimization using both distance and torsional constraints. Side chains were then added as their experimentally preferred rotamers and the whole molecule minimized without constraints to give the final model structure. This shows type II' and III ß turns at residues 4–5 and 9–10, respectively, coupled by twisted antiparallel strands which show hydrogen bonds between all four pairs of opposing peptide groups. The backbone conformation of residues 2–6 closely resembles that found in the crystal structure of gramicidin S.     

Implementation and use of the method of prudent ascent in conformational analysis using molecular mechanics

In two-dimensional conformational analysis the current practice is to perform an energy minimization for all possible combinations of two dihedral angles in the molecule, in a fixed order, and apply a certain dihedral angle step-size. A newly developed method is presented in which the order of the evaluation points on the energy-surface is not fixed, but is dependent on all previous results in a way which we call “the method of prudent ascent.” In this method the most promising calculation is carried out first, thus minimizing the risk of atomic collisions. In order to be able to take care of the many additional degrees of conformational freedom present in, e.g., carbohydrate molecules, all minimizations are performed using a set of different promising starting conformations on the basis of previous calculations, and only the lowest energy result for each point is saved. An application of the method to conformational analysis of methyl-cellobiose and the artificial sweetener trichlorogalactosucrose is also presented.     

Global methods in classical mechanics. The Euler–Lagrange equation

Starting from a smooth manifold Q as configurational space, the intrinsec form of Euler–Lagrange equation is derived using a differential geometrical approach in order to obtain a relation valid on the whole tangent bundle TQ that constitutes the phase space of a generical mechanical system. This revised version was published online in July 2006 with corrections to the Cover Date.     

On the use of electrostatic potential derived charges in molecular mechanics force fields. The relative solvation free energy of cis- and trans-N-methyl-acetamide

We have carried out free energy perturbation calculations on the relative solvation free energy of cis- and trans-N-methyl-acetamide (NMA). Experimentally, the solvation free energy difference has been found to be near zero. Using 6-31G* ab initio electrostatic potential derived charges for both the cis and trans conformations, we calculate a solvation free energy difference of 0.1 ± 0.1 kcal/mol. Using the 6-31G* charges derived for the trans conformation for both the cis and trans models leads to a solvation free energy difference of 0.9 ± 0.1 kcal/mol, compared to the value of 2.2 kcal/mol determined for the OPLS model for trans-NMA.     

The molecular mechanics of valinomycin. II: Comparative studies of alkali ion binding

We have carried out a series of molecular mechanics calculations on the alkali ion complexes of valinomycin. For the ions Na⁺, K⁺, Rb⁺, and Cs⁺ we have found three-fold rotationally symmetric conformations as the lowest energy structures, while for Li⁺ a markedly asymmetric configuration is preferred. The relative free energies of the complexes show that Li⁺ is by far the poorest binding partner in solution, followed by Na⁺, which is in turn far poorer than any of the three larger ions. The binding selectivity derives from the slower variation of the complexation free energy with ionic size than the ionic solvation free energy, so that the ionophore is unable to compete with the solvent for the smaller ions. Our calculated strain energies suggest that valinomycin's failure to form complexes with the smaller ions in solution is due partially to the rigidity of the ionophore structure, which prevents the central cavity from contracting to accommodate them. Certain geometric criteria indicate that K⁺ provides the best fit to the binding site, although there is some inconsistency between the energetic and geometric criteria of binding ability.     

Symmetry in the Heitler-London method. II. The determination of the allowed molecular multiplets

In the framework of the Heitler-London method, a method for determining the allowed molecular multiplets is proposed. The method is based on the connection of the total molecular spin with the permutation symmetry of the coordinate wave function and on the isomorphism of the molecular point group with a certain subgroup of the electron permutation group. The method does not depend on the approximation in which the molecular ions are considered.     

Optimizing the performance of the multiconfiguration molecular mechanics method

Multiconfiguration molecular mechanics (MCMM) is a general algorithm for constructing potential energy surfaces for reactive systems (Kim, Y.; Corchado, J. C.; Villà, J.; Xing, J.; Truhlar, D. G. J. Chem. Phys. 2000, 112, 2718). This paper illustrates how the performance of the MCMM method can be improved by adopting improved molecular mechanics parameters. We carry out calculations of reaction rate constants using variational transition state theory with optimized multidimensional tunneling on the MCMM PESs for three hydrogen transfer reactions, and we compare the results to direct dynamics. We find that the MCMM method with as little as one electronic structure Hessian can describe the dynamically important regions of the ground-electronic state PES, including the corner-cutting-tunneling region of the reaction swath, with practical accuracy.     

Atomic and molecular quantum mechanics by the path integral molecular dynamics method

The quantum path integral molecular dynamics method was applied to studies of excess electron localization by a Na⁺ ion and by a NaCl molecule. Spatial and energetic characterization of the ground state of the excess electron compare favorably with results of model potential calculations for Na and with SCF Cl calculations for NaCl⁻.     

Assessing the performance of the molecular mechanics/Poisson Boltzmann surface area and molecular mechanics/generalized Born surface area methods. II. The accuracy of ranking poses generated from docking

In molecular docking, it is challenging to develop a scoring function that is accurate to conduct high-throughput screenings. Most scoring functions implemented in popular docking software packages were developed with many approximations for computational efficiency, which sacrifices the accuracy of prediction. With advanced technology and powerful computational hardware nowadays, it is feasible to use rigorous scoring functions, such as molecular mechanics/Poisson Boltzmann surface area (MM/PBSA) and molecular mechanics/generalized Born surface area (MM/GBSA) in molecular docking studies. Here, we systematically investigated the performance of MM/PBSA and MM/GBSA to identify the correct binding conformations and predict the binding free energies for 98 protein-ligand complexes. Comparison studies showed that MM/GBSA (69.4%) outperformed MM/PBSA (45.5%) and many popular scoring functions to identify the correct binding conformations. Moreover, we found that molecular dynamics simulations are necessary for some systems to identify the correct binding conformations. Based on our results, we proposed the guideline for MM/GBSA to predict the binding conformations. We then tested the performance of MM/GBSA and MM/PBSA to reproduce the binding free energies of the 98 protein-ligand complexes. The best prediction of MM/GBSA model with internal dielectric constant 2.0, produced a Spearman's correlation coefficient of 0.66, which is better than MM/PBSA (0.49) and almost all scoring functions used in molecular docking. In summary, MM/GBSA performs well for both binding pose predictions and binding free-energy estimations and is efficient to re-score the top-hit poses produced by other less-accurate scoring functions.     

Quantum mechanics/molecular mechanics study of the catalytic cycle of water splitting in photosystem II

This paper investigates the mechanism of water splitting in photosystem II (PSII) as described by chemically sensible models of the oxygen-evolving complex (OEC) in the S0-S4 states. The reaction is the paradigm for engineering direct solar fuel production systems since it is driven by solar light and the catalyst involves inexpensive and abundant metals (calcium and manganese). Molecular models of the OEC Mn3CaO4Mn catalytic cluster are constructed by explicitly considering the perturbational influence of the surrounding protein environment according to state-of-the-art quantum mechanics/molecular mechanics (QM/MM) hybrid methods, in conjunction with the X-ray diffraction (XRD) structure of PSII from the cyanobacterium Thermosynechococcus elongatus. The resulting models are validated through direct comparisons with high-resolution extended X-ray absorption fine structure spectroscopic data. Structures of the S3, S4, and S0 states include an additional mu-oxo bridge between Mn(3) and Mn(4), not present in XRD structures, found to be essential for the deprotonation of substrate water molecules. The structures of reaction intermediates suggest a detailed mechanism of dioxygen evolution based on changes in oxidization and protonation states and structural rearrangements of the oxomanganese cluster and surrounding water molecules. The catalytic reaction is consistent with substrate water molecules coordinated as terminal ligands to Mn(4) and calcium and requires the formation of an oxyl radical by deprotonation of the substrate water molecule ligated to Mn(4) and the accumulation of four oxidizing equivalents. The oxyl radical is susceptible to nucleophilic attack by a substrate water molecule initially coordinated to calcium and activated by two basic species, including CP43-R357 and the mu-oxo bridge between Mn(3) and Mn(4). The reaction is concerted with water ligand exchange, swapping the activated water by a water molecule in the second coordination shell of calcium.     

On the use of the interaction picture in classical mechanics

The interaction picture is formulated as a canonical transformation in classical mechanics. It is shown that often great computational saving may be achieved by integrating Hamilton's equation in this representation.     

The ellipsoid algorithm as a method for the determination of polypeptide conformations from experimental distance constraints and energy minimization

A new method for constrained nonlinear optimization known as the ellipsoid algorithm is evaluated as a means of determining and refining the conformations of peptides. Advantages of the ellipsoid algorithm over conventional optimization methods include that it avoids many local minima that other methods would be trapped by, and that it is sometimes able to find optimum solutions in which the constraints are satisfied exactly. The dihedral angles about single bonds were used as variables to keep the dimensionality low (the rate of convergence decreases rapidly with increasing dimensionality of the problem). The method is evaluated on problems involving distance constraints, and for minimization of conformational energy functions. In an initial application, conformations consistent with an experimental set of NMR distance constraints were obtained in a problem involving 48 variable dihedral angles.     

An efficient method for the calculation of quantum mechanics/molecular mechanics free energies

The combination of quantum mechanics (QM) with molecular mechanics (MM) offers a route to improved accuracy in the study of biological systems, and there is now significant research effort being spent to develop QM/MM methods that can be applied to the calculation of relative free energies. Currently, the computational expense of the QM part of the calculation means that there is no single method that achieves both efficiency and rigor; either the QM/MM free energy method is rigorous and computationally expensive, or the method introduces efficiency-led assumptions that can lead to errors in the result, or a lack of generality of application. In this paper we demonstrate a combined approach to form a single, efficient, and, in principle, exact QM/MM free energy method. We demonstrate the application of this method by using it to explore the difference in hydration of water and methane. We demonstrate that it is possible to calculate highly converged QM/MM relative free energies at the MP2/aug-cc-pVDZ/OPLS level within just two days of computation, using commodity processors, and show how the method allows consistent, high-quality sampling of complex solvent configurational change, both when perturbing hydrophilic water into hydrophobic methane, and also when moving from a MM Hamiltonian to a QM/MM Hamiltonian. The results demonstrate the validity and power of this methodology, and raise important questions regarding the compatibility of MM and QM/MM forcefields, and offer a potential route to improved compatibility.     

On the space of eigenvectors of molecular quantum mechanics

As the representative point on an electronic energy surface moves along any closed level curve, the corresponding electronic eigenvectors behave as does the normal to a Moebius strip when its foot is moved around the strip. This introduces a topology into the space of eigenvectors that is not envisaged in the definition of Hilbert space, but is fundamental to molecular quantum mechanics. Salient properties of these “non-Hilbertian” vector spaces are established. © 1994 John Wiley & Sons, Inc.     

Some properties of the Lagrange multiplier μ in density functional theory

Some new properties of the Lagrange multiplier μ introduced through the normalization constraint on ρ in the variations of energy density functionals are determined. Through arguments concerning the homogeneity properties of these functionals with respect to μ, it is demonstrated that at the point of variation μ = μ₀ = E₀/N, where E₀ is the ground state energy and N is the total particle number. It is also shown that the value of μ₀ is independent of the normalization imposed on ρ. The interpretation of μ₀ as a chemical potential is discussed in the light of these findings.