Showing papers in "Journal of Computational Chemistry in 1998"

Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function

Abstract: A novel and robust automated docking method that predicts the bound conformations of flexible ligands to macromolecular targets has been developed and tested, in combination with a new scoring function that estimates the free energy change upon binding. Interestingly, this method applies a Lamarckian model of genetics, in which environmental adaptations of an individual's phenotype are reverse transcribed into its genotype and become . heritable traits sic . We consider three search methods, Monte Carlo simulated annealing, a traditional genetic algorithm, and the Lamarckian genetic algorithm, and compare their performance in dockings of seven protein)ligand test systems having known three-dimensional structure. We show that both the traditional and Lamarckian genetic algorithms can handle ligands with more degrees of freedom than the simulated annealing method used in earlier versions of AUTODOCK, and that the Lamarckian genetic algorithm is the most efficient, reliable, and successful of the three. The empirical free energy function was calibrated using a set of 30 structurally known protein)ligand complexes with experimentally determined binding constants. Linear regression analysis of the observed binding constants in terms of a wide variety of structure-derived molecular properties was performed. The final model had a residual standard y1 y1 .

Geometry optimization of molecular structures in solution by the polarizable continuum model

Abstract: A new implementation of analytical gradients for the polarizable continuum model is presented, which allows Hartree‐Fock and density functional calculations taking into account both electrostatic and nonelectrostatic contributions to energies and gradients for closed and open shell systems. Simplified procedures neglecting the derivatives of the cavity surface and/or using single spheres for XHn groups have also been implemented and tested. The solvent‐induced geometry relaxation has been studied for a number of representative systems in order to test the efficiency of the procedure and to investigate the role of different contributions. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 404–417, 1998

Exact and efficient analytical calculation of the accessible surface areas and their gradients for macromolecules

Abstract: A new method for exact analytical calculation of the accessible surface areas and their gradients with respect to atomic coordinates is described. The new surface routine, GETAREA, finds solvent-exposed vertices of intersecting atoms, and thereby avoids calculating buried vertices which are not needed to determine the accessible surface area by the Gauss)Bonnet theorem. The surface routine was implemented in FANTOM, a program for energy minimization and Monte Carlo simulation, and tested for accuracy and efficiency in extensive energy minimizations of Met-enkephalin, the a-amylase inhibitor . tendamistat, and avian pancreatic polypeptide APP . The CPU time for the exact calculation of the accessible surface areas and their gradients has been . . reduced by factors of 2.2 Met-enkephalin and 3.2 tendamistat compared with our previous approach. The efficiency of our exact method is similar to the recently described approximate methods MSEED and SASAD. The performance of several atomic solvation parameter sets was tested in searches for low energy conformations of APP among conformations near the native X-ray crystal structure and highly distorted structures. The protein solvation parameters from w . x

Natural resonance theory: I. General formalism

Abstract: We present a new quantum-mechanical resonance theory based on the first-order reduced density matrix and its representation in terms of natural bond orbitals (NBOs). This “natural” resonance theory (NRT) departs in important respects from the classical Pauling-Wheland formulation, yet it leads to quantitative resonance weights that are in qualitative accord with conventional resonance theory and chemical intuition. The NRT variational functional leads to an optimal resonance-weighted approximation to the full density matrix, combining the “single reference” limit of weak delocalization (incorporating diagonal population changes only) with the full “multireference” limit of strong delocalization (incorporating off-diagonal couplings between resonance structures. The NRT variational functional yields an error measure that serves as an intrinsic criterion of accuracy of the resonance-theoretic description. The NRT program structure, algorithms, and numerical characteristics are described in supplementary material, and detailed chemical applications are presented in two companion papers. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 593–609, 1998

Performance of parallel TURBOMOLE for density functional calculations

Abstract: The parallelization of density functional treatments of molecular electronic energy and first-order gradients is described, and the performance is documented. The quadrature required for exchange correlation terms and the treatment of exact Coulomb interaction scales virtually linearly up to 100 nodes. The RI-J technique to approximate Coulomb interactions (by means of an auxiliary basis set approximation for the electron density) even shows superlinear speedup on distributed memory architectures. The bottleneck is then linear algebra. Demonstrative application examples include molecules with up to 300 atoms and 3000 basis functions that can now be treated in a few hours per geometry optimization cycle in C1 symmetry. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1746–1757, 1998

Natural resonance theory: II. Natural bond order and valency

Abstract: Resonance weights derived from the Natural Resonance Theory .NRT , introduced in the preceding paper are used to calculate ''natural bond order,'' ''natural atomic valency,'' and other atomic and bond indices reflecting the resonance composition of the wave function. These indices are found to give significantly better agreement with observed properties empirical valency, . bond lengths than do corresponding MO-based indices. A characteristic feature of the NRT treatment is the description of bond polarity by a ''bond ionicity'' . index resonance-averaged NBO polarization ratio , which replaces the ''covalent-ionic resonance'' of Pauling-Wheland theory and explicity exhibits the complementary relationship of covalency and electrovalency that underlies empirical assignments of atomic valency. We present ab initio NRT applications . to prototype saturated and unsaturated molecules methylamine, butadiene , . polar compounds fluoromethanes , and open-shell species: hydroxymethyl . radical to demonstrate the numerical stability, convergence, and chemical reasonableness of the NRT bond indices in comparison to other measures of valency and bond order in current usage. Q 1998 John Wiley & Sons, Inc. J Comput Chem 19: 610)627, 1998

Natural resonance theory: III. Chemical applications

Abstract: We describe quantitative numerical applications of the natural resonance theory (NRT) to a variety of chemical bonding types, in order to demonstrate the generality and practicality of the method for a wide range of chemical systems. Illustrative applications are presented for (1) benzene and polycyclic aromatics; (2) CO2, formate, and related acyclic species; (3) ionic and polar compounds; (4) coordinate covalent compounds and complexes; (5) hypervalent and electron-deficient species; (6) noncovalent H-bonded complex; and (7) a model Diels-Alder chemical reaction surface. The examples exhibit the general harmony of NRT weightings with qualitative resonance-theoretic concepts and illustrate how these concepts can be extended to many new types of chemical phenomena at a quanitative ab initio level. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 628–646, 1998

Temperature dependence of TIP3P, SPC, and TIP4P water from NPT Monte Carlo simulations: Seeking temperatures of maximum density

Abstract: Monte Carlo statistical mechanics simulations have been carried out with the TIP3P, SPC, and TIP4P models for liquid water at 13 temperatures from −50°C to 100°C at 1 atm. Long runs with 512 water molecules provided definitive results for densities. Although the TIP4P model yields a flat region in the density profile and a temperature of maximum density near −15°C, the SPC and TIP3P models show monotonically increasing density with decreasing temperature. Results for heats of vaporization, isothermal compressibilities, and coefficients of thermal expansion and their convergence characteristics are also reported. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1179–1186, 1998

Parametrization of aliphatic CHn united atoms of GROMOS96 force field

Abstract: The derivation of the van der Waals parameters for the aliphatic CHn united atoms of the GROMOS96 force field is presented. The parameters have been adjusted to reproduce the experimental enthalpies of vaporization and vapor pressures or densities of a set of nine alkanes in the liquid state at 298 K (or at the boiling point in the case of methane), using a cutoff radius for the van der Waals interactions of 1.6 nm. Force fields to be used in molecular simulations are bound to the conditions chosen for their parametrization, for example, the temperature, the densities of the systems included in the calibration set, or the cutoff radius used for the nonbonded interactions. Van der Waals parameters for the CHn united atoms of earlier GROMOS force fields were developed using a cutoff radius of 0.8 nm for the van der Waals interactions. Because the van der Waals interaction energy between aliphatic groups separated by distances between 0.8 and 1.4 nm is not negligible at liquid densities, the use of these parameters in combination with longer cutoffs leads to an overestimation of the attractive van der Waals interaction energy. The relevance of this excess attraction depends on the size of the groups that are interacting, as well as on their local densities. Free energies of hydration have been calculated for five alkanes. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 535–547, 1998

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

Abstract: 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

Comparison of methods for deriving atomic charges from the electrostatic potential and moments

Abstract: Four methods for deriving partial atomic charges from the quantum chemical electrostatic potential CHELP, CHELPG, Merz-Kollman, and . RESP have been compared and critically evaluated. It is shown the charges strongly depend on how and where the potential points are selected. Two alternative methods are suggested to avoid the arbitrariness in the point-selection schemes and van der Waals exclusion radii: CHELP-BOW, which also estimates the charges from the electrostatic potential, but with potential points that are Boltzmann-weighted after their occurrence in actual simulations using the energy function of the program in which the charges will be used, and CHELMO, which estimates the charges directly from the electrostatic multipole moments. Different criteria for the quality of the charges are discussed. The CHELMO method gives the best multipole moments for small and medium-sized polar systems, whereas the CHELP-BOW charges reproduce best the total interaction energy in actual simulations. Among the standard methods, the Merz-Kollman charges give the best moments and potentials, but they show an appreciable dependence on the orientation of the molecule. We have also examined the recent warning that charges derived by a least-squares fit to the electrostatic potential normally are not statistically valid. It is shown that no rank-deficiency problems are encountered for molecules with up to 84 atoms if the least-squares fit is performed using pseudoinverses calculated by singular value decomposition and if constraints are treated by elimination. Q 1998 John Wiley & Sons, Inc. J Comput Chem 19: 377)395, 1998

TheRate: Program for Ab Initio Direct Dynamics Calculations of Thermal and Vibrational-State-Selected Rate Constants

Abstract: We introduce TheRate THEoretical RATEs , a complete . application program with a graphical user interface GUI for calculating rate constants from first principles. It is based on canonical variational transition-state . theory CVT augmented by multidimensional semiclassical zero and small . curvature tunneling approximations. Conventional transition-state theory TST with one-dimensional Wigner or Eckart tunneling corrections is also available. Potential energy information needed for the rate calculations are obtained from ab initio molecular orbital andror density functional electronic structure theory. Vibrational-state-selected rate constants may be calculated using a diabetic model. TheRate also introduces several technical advancements, namely the focusing technique and energy interpolation procedure. The focusing technique minimizes the number of Hessian calculations required by distributing more Hessian grid points in regions that are critical to the CVT and tunneling calculations and fewer Hessian grid points elsewhere. The energy interpolation procedure allows the use of a computationally less demanding electronic structure theory such as DFT to calculate the Hessians and geometries, while the energetics can be improved by performing a small number of single-point energy calculations along the MEP at a more accurate level of theory. The CH q H l CH q H reaction is used as a model to demonstrate usage of the 43 2 program, and the convergence of the rate constants with respect to the number of electronic structure calculations. Q 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1039)1052, 1998

Combined ab initio/empirical approach for optimization of Lennard–Jones parameters

Abstract: Obtaining accurate Lennard–Jones (L-J) parameters is a vital part of the optimization of empirical force fields due to their significant contribution to condensed-phase properties. We present a novel approach to optimize L-J parameters via the use of ab initio data on interactions between rare gas atoms and model compounds combined with the reproduction of experimental pure solvent properties. Relative values of ab initio minimum interaction energies and geometries between helium or neon and model compounds were used to optimize the relative magnitude of the L-J parameters. Absolute values of the L-J parameters were determined by reproducing experimental heats of vaporization and molecular volumes for pure solvents. Application of the approach was performed on methane, ethane, and propane. Free energies of aqueous solvation and butane pure solvent and aqueous solvation calculations were used to test the developed L-J parameters. The new alkane parameters are similar or improved as compared with current empirical force field parameters with respect to experimental pure solvent properties and free energies of aqueous solvation. Also included is a description of the internal portion of the force field. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 334–348, 1998

Searching for saddle points of potential energy surfaces by following a reduced gradient

Abstract: The old coordinate driving procedure to find transition structures . in chemical systems is revisited. The well-known gradient criterion, =E x s 0, . which defines the stationary points of the potential energy surface PES , is reduced by one equation corresponding to one search direction. In this manner, abstract curves can be defined connecting stationary points of the PES. Starting at a given minimum, one follows a well-selected coordinate to reach the saddle of interest. Usually, but not necessarily, this coordinate will be related to the . reaction progress. The method, called reduced gradient following RGF , locally has an explicit analytical definition. We present a predictor)corrector method for tracing such curves. RGF uses the gradient and the Hessian matrix or updates of the latter at every curve point. For the purpose of testing a whole surface, the six-dimensional PES of formaldehyde, H CO, was explored by RGF 2 . using the restricted Hartree)Fock RHF method and the STO-3G basis set. Forty-nine minima and saddle points of different indices were found. At least seven stationary points representing bonded structures were detected in addition to those located using another search algorithm on the same level of theory. Further examples are the localization of the saddle for the HCN | CNH . isomerization used for steplength tests and for the ring closure of azidoazo- methine to 1 H-tetrazole. The results show that following the reduced gradient may represent a serious alternative to other methods used to locate saddle points in quantum chemistry. Q 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1087)1100, 1998

United-residue force field for off-lattice protein-structure simulations: III. Origin of backbone hydrogen-bonding cooperativity in united-residue potentials

Abstract: Based on the dipole model of peptide groups developed in our earlier work [Liwo et al., Prot. Sci.,2, 1697 (1993)], a cumulant expansion of the average free energy of the system of freely rotating peptide-group dipoles tethered to a fixed α-carbon trace is derived. A graphical approach is presented to find all nonvanishing terms in the cumulants. In particular, analytical expressions for three- and four-body (correlation) terms in the averaged interaction potential of united peptide groups are derived. These expressions are similar to the cooperative forces in hydrogen bonding introduced by Kolinski and Skolnick [J. Chem. Phys.,97, 9412 (1992)]. The cooperativity arises here naturally from the higher order terms in the power-series expansion (in the inverse of the temperature) for the average energy. Test calculations have shown that addition of the derived four-body term to the statistical united-residue potential of our earlier work [Liwo et al., J. Comput. Chem.,18, 849, 874 (1997)] greatly improves its performance in folding poly-l-alanine into an α-helix. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 259–276, 1998

Comparative study of BSSE correction methods at DFT and MP2 levels of theory

Abstract: A comparative study of intermolecular potential energy curves is performed on the complexes H2O(SINGLE BOND)HF, H2O(SINGLE BOND)H2O, H2O(SINGLE BOND)H2S, and H2S(SINGLE BOND)H2S using nine different basis sets at the MP2 and DFT (BLYP and B3LYP) levels of theory. The basis set superposition error is corrected by means of the counterpoise scheme and based on the “chemical Hamiltonian approach.” The counterpoise and CHA-corrected DFT curves are generally close to each other. Using small basis sets, the B3LYP functional cannot be favored against the BLYP one because the BLYP results sometimes get closer to the MP2 values than those of B3LYP. From the results—including the available literature data—we suggest that one has to use at least polarized-valence triple-zeta-quality basis sets (TZV, 6-311G) for the investigation of hydrogen-bonded complexes. Special attention must be paid to the physical nature of the binding. If the dispersion forces become significant DFT methods are not able to describe the interaction. Proper correction for the basis set superposition error is found to be mandatory in all cases. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 575–584, 1998

Assessing energy functions for flexible docking

Abstract: A good docking algorithm requires an energy function that is selective, in that it clearly differentiates correctly docked structures from misdocked ones, and that is efficient, meaning that a correctly docked structure can be identified quickly. We assess the selectivity and efficiency of a broad spectrum of energy functions, derived from systematic modifications of the CHARMM param19/toph19 energy function. In particular, we examine the effects of the dielectric constant, the solvation model, the scaling of surface charges, reduction of van der Waals repulsion, and nonbonded cutoffs. Based on an assessment of the energy functions for the docking of five different ligand–receptor complexes, we find that selective energy functions include a variety of distance-dependent dielectric models together with truncation of the nonbonded interactions at 8 A. We evaluate the docking efficiency, the mean number of docked structures per unit of time, of the more selective energy functions, using a simulated annealing molecular dynamics protocol. The largest improvements in efficiency come from a reduction of van der Waals repulsion and a reduction of surface charges. We note that the most selective potential is quite inefficient, although a hierarchical approach can be employed to take advantage of both selective and efficient energy functions. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1612–1622, 1998

Derivation of class II force fields: V. Quantum force field for amides, peptides, and related compounds

Abstract: As the field of biomolecular structure advances, there is an ever-growing need for accurate modeling of molecular energy surfaces to simulate and predict the properties of these important systems. To address this need, a second generation amide force field for use in simulations of small organics as well as proteins and peptides has been derived. The critical question of what accuracy can be expected from calculations in general, and with this class II force field in particular, is addressed for structural, dynamic, and energetic properties. The force field is derived from a recent methodology we have developed that involves the systematic use of quantum mechanical observables. Systematic ab initio calculations were carried out for numerous configurations of 17 amide and related compounds. Relative energies and first and second derivatives of the energy of 638 structures of these compounds resulted in 140,970 ab initio quantum mechanical observables. The class II peptide quantum mechanical force field (QMFF), containing 732 force constants and reference values, was parameterized against these observables. A major objective of this work is to help establish the role of anharmonicity and coupling in improving the accuracy of molecular force fields, as these terms have not yet become an agreed upon standard in the ever more extensive simulations being used to probe biomolecular properties. This has been addressed by deriving a class I harmonic diagonal force field (HDFF), which was fit to the same energy surface as the QMFF, thus providing an opportunity to quantify the effects of these coupling and anharmonic contributions. Both force field representations are assessed in terms of their ability to fit the observables. They have also been tested by calculating the properties of 11 stationary states of these amide molecules. Optimized structures, vibrational frequencies, and conformational energies obtained from the quantum calculations and from both the QMFF and the HDFF are compared. Several strained and derivatized compounds including urea, formylformamide, and butyrolactam are included in these tests to assess the range of applicability (transferability) of the force fields. It was found that the class II coupled anharmonic force field reproduced the structures, energies, and vibrational frequencies significantly more faithfully than the class I harmonic diagonal force field. An important measure, rms energy deviation, was found to be 1.06 kcal/mol with the class II force field, and 2.30 kcal/mol with the harmonic diagonal force field. These deviations represent the error in relative configurational energy differences for strained and distorted structures calculated with the force fields compared with quantum mechanics. This provides a measure of the accuracy that might be expected in applications where strain may be important such as calculating the energy of a system as it approaches a (rotational) barrier, in ligand binding to a protein, or effects of introducing substituents into a molecule that may induce strain. Similar results were found for structural properties. Protein dynamics is becoming of ever-increasing interest, and, to simulate dynamic properties accurately, the dynamic behavior of model compounds needs to be well accounted for. To this end, the ability of the class I and class II force fields to reproduce the vibrational frequencies obtained from the quantum energy surface was assessed. An rms deviation of 43 cm−1 was achieved with the coupled anharmonic force field, as compared to 105 cm−1 with the harmonic diagonal force field. Thus, the analysis presented here of the class II force field for the amide functional group demonstrates that the incorporation of anharmonicity and coupling terms in the force field significantly improves the accuracy and transferability with regard to the simulation of structural, energetic, and dynamic properties of amides. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 430–458, 1998

Protonation of Thymine, Cytosine, Adenine, and Guanine DNA Nucleic Acid Bases: Theoretical Investigation into the Framework of Density Functional Theory

Abstract: Gradient-corrected density functional computations with triple-zeta-type basis sets were performed to determine the preferred protonation site and the absolute gas-phase proton affinities of the most stable tautomer of the DNA bases thymine (T), cytosine (C), adenine (A), and guanine (G). Charge distribution, bond orders, and molecular electrostatic potentials were considered to rationalize the obtained results. The vibrational frequencies and the contribution of the zero-point energies were also computed. Significant geometrical changes in bond lengths and angles near the protonation sites were found. At 298 K, proton affinities values of 208.8 (T), 229.1 (C), 225.8 (A), and 230.3 (G) kcal/mol were obtained in agreement with experimental results. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 989–1000, 1998

Factors controlling relative stability of anomers and hydroxymethyl conformers of glucopyranose

Abstract: The relative energies of 11 different conformers of D-glucose, including different exo-anomeric orientations and at least one of each hydroxymethyl conformer (G−, G+, and T) for each of the two anomeric forms (α and β), were calculated at much more complete levels of quantum mechanical (QM) electronic structure theory than previously available, and relative free energies in solution were calculated by adding rotational, vibrational, and aqueous solvation effects. The gas-phase results are based on very large basis sets (up to 624 contracted basis functions) and the coupled cluster method for electron correlation. Solvation Model 5.4/ AM1 was used to calculate the effects of aqueous solvation. Factors contributing to the relative energies of these conformers have been analyzed. Relative energies varied considerably (up to 4.5 kcal/mol), depending on the theoretical level, and different levels of theory disagreed as to which anomer was lower in energy. The highest-level gas-phase calculations predicted the α-anomer to be lower in free energy by 0.4 kcal/mol (Boltzmann average). Gas-phase energies from several different classical force fields were compared to QM results. The QM structures optimized at the MP2/cc-pVDZ level of theory compared well with experiment for three different crystal structures. In water, the β-anomers were better solvated than the α-anomers by 0.6 kcal/mol (Boltzmann average). Contributions of individual hydrophilic groups to the solvation free energies were analyzed. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1111–1129, 1998

Docking small ligands in flexible binding sites

Abstract: A novel procedure for docking ligands in a flexible binding site is presented. It relies on conjugate gradient minimization, during which nonbonded interactions are gradually switched on. Short Monte Carlo minimization runs are performed on the most promising candidates. Solvation is implicitly taken into account in the evaluation of structures with a continuum model. It is shown that the method is very accurate and can model induced fit in the ligand and the binding site. The docking procedure has been successfully applied to three systems. The first two are the binding of progesterone and 5b-androstane-3,17-dione to the antigen binding fragment of a steroid binding antibody. A comparison of the crystal structures of the free and the two complexed forms reveals that any attempt to model binding must take protein rearrangements into account. Furthermore, the two ligands bind in two different orientations, posing an additional challenge. The third test case is the docking of a . N- 2-naphthyl-sulfonyl-glycyl -D-para-amidino-phenyl-alanyl-piperidine .NAPAP to human a-thrombin. In contrast to steroids, NAPAP is a very flexible ligand, and no information of its conformation in the binding site is used. All docking calculations are started from X-ray conformations of proteins with the uncomplexed binding site. For all three systems the best minima in terms of free energy have a root mean square deviation from the X-ray structure ˚ smaller than 1.5 A for the ligand atoms. Q 1998 John Wiley & Sons, Inc. J Comput Chem 19 :2 1 )37, 1998

Directional hydrogen bonding in the mm3 force field : ii

Abstract: Extensive calculations on hydrogen bonded systems were carried out using the improved MM3 directional hydrogen bond potential. The resulting total function was reoptimized. Comparisons of the hydrogen bonding potential function from ab initio calculations (MP2/6-31G**); the original MM3(89); and the reoptimized MM3 force field MM3(96), for a variety of C, N, O, and Cl systems including the formamide dimer and formamide–water complex, are described herein. Hydrogen bonding is shown to be a far more complicated and ubiquitous phenomenon than is generally recognized. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1001–1016, 1998

Assessing search strategies for flexible docking

Abstract: We assess the efficiency of molecular dynamics (MD), Monte Carlo (MC), and genetic algorithms (GA) for docking five representative ligand–receptor complexes. All three algorithms employ a modified CHARMM-based energy function. The algorithms are also compared with an established docking algorithm, AutoDock. The receptors are kept rigid while flexibility of ligands is permitted. To test the efficiency of the algorithms, two search spaces are used: an 11-A-radius sphere and a 2.5-A-radius sphere, both centered on the active site. We find MD is most efficient in the case of the large search space, and GA outperforms the other methods in the small search space. We also find that MD provides structures that are, on average, lower in energy and closer to the crystallographic conformation. The GA obtains good solutions over the course of the fewest energy evaluations. However, due to the nature of the nonbonded interaction calculations, the GA requires the longest time for a single energy evaluation, which results in a decreased efficiency. The GA and MC search algorithms are implemented in the CHARMM macromolecular package. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1623–1631, 1998

Reference program for molecular calculations with slater-type orbitals

Abstract: A program for computing all the integrals appearing in molecular calculation with Slater-type orbitals is reported. The program is mainly intended as a reference for testing and comparing other algorithms and techniques. An analysis of the performance of the program is presented, paying special attention to the computational cost and the accuracy of the results. Results are also compared with others obtained with Gaussian basis sets of similar quality. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1284–1293, 1998

Radical performance enhancements for combinatorial optimization algorithms based on the dead‐end elimination theorem

Abstract: Recent advances in protein design have demonstrated the effectiveness of optimization algorithms based on the dead-end elimination theorem. The algorithms solve the combinatorial problem of finding the optimal placement of side chains for a set of backbone coordinates. Although they are powerful tools, these algorithms have severe limitations when the number of side chain rotamers is large. This is due to the high-order time dependence of the aspect of the calculation that deals with rotamer doubles. We present three independent algorithmic enhancements that significantly increase the speed of the doubles computation. These methods work by using quantities that are inexpensive to compute as a basis for forecasting which expensive calculations are worthwhile. One of the methods, the comparison of extrema, is derived from analytical considerations, and the remaining two, the “magic-bullet” and the “qrs” and “quv” metrics, are based on empirical observation of the distribution of energies in the system. When used together, these methods effect an overall speed improvement of as much as a factor of 47, and for the doubles aspect of the calculation, a factor of 95. Together, these enhancements extend the envelope of inverse folding to larger proteins by making formerly intractable calculations attainable in reasonable computer time. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1505–1514, 1998

Fast summation boundary element method for calculating solvation free energies of macromolecules

Abstract: We present a boundary element method (BEM) for calculating the reaction field energy of a macromolecule embedded in a high-dielectric medium such as water. In a BEM calculation, the key computational task is the calculation of the induced surface charge distribution at the dielectric boundary. This is obtained by solving a system of linear equations whose dimension can run into the tens of thousands for a macromolecule. In this work, we use a fast summation hierarchical multipole method to solve for the induced surface charge densities. By careful analysis of the levels of approximation required for the various terms in the calculation, we avoid the unnecessary computation of terms that contribute negligibly to the final outcome and, consequently, achieve high computational efficiency. For a protein such as BPTI with 890 atoms, the calculation of the induced surface charge density distribution and the reaction field energy was completed in 7.9 s on an SGI workstation with an R10000 CPU. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1494–1504, 1998

Automated molecular mechanics parameterization with simultaneous utilization of experimental and quantum mechanical data

Abstract: A general set of procedures for automated parameterization of molecular mechanics force fields is presented. The current implementation for MacroModel force fields can easily be adapted to other programs running under Unix. Several variants of numerical Newton–Raphson and simplex methodologies are used to fit simultaneously both experimental and quantum mechanical data. In a sample parameterization of ethane, 24 different modes of combining various sources of reference data are compared. The best cost/performance ratio for generating quantum mechanical source data was obtained with the B3LYP/6‐31G* method. The best optimization strategy consists of initial subset optimizations with a modified simplex method, followed by Newton–Raphson optimization using Lagrange multipliers. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1146–1166, 1998

Solvation effects on reaction profiles by the polarizable continuum model coupled with the Gaussian density functional method

Abstract: An efficient version of the polarizable continuum model for solvation has been implemented in the Gaussian density-functional-based code called deMon. Solvation free energies of representative compounds have been calculated as a preliminary test. The hydration effects on the reaction profile of the Cl−+CH3ClClCH3+Cl− reaction and the thermodynamics of the Menschutkin reaction have also been investigated. Finally, the conformational behavior of the 1,2-diazene cis–trans isomerization process in water was examined. Comparisons between the results obtained and the available experimental data and previous theoretical computations have been made. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 290–299, 1998

Implementation and validation of the Lacks-Gordon exchange functional in conventional density functional and adiabatic connection methods

Abstract: We present an analysis of the numerical performances of the exchange functional proposed by Lacks and Gordon, which we have implemented in the Gaussian series of programs. This functional has been built with the double aim of respecting most of the known scaling and asymptotic properties and of giving good numerical performances, especially as concerns noncovalent interactions. We have found that the coupling of the Lacks-Gordon exchange and Lee-Yang-Parr correlation functionals provides a reliable conventional density functional approach. The corresponding parameter-free adiabatic connection model, in which the ratio between Hartree-Fock and Lacks-Gordon exchange is determined a priori from purely theoretical considerations, allows us to obtain remarkable results for both covalent and noncovalent interactions in a satisfactory theoretical scheme, encompassing the free electron gas limit and most of the known scaling conditions. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 418–429, 1998

Molecular dynamics simulations of sorption of organic compounds at the clay mineral / aqueous solution interface

Abstract: The adsorption of trichloroethene, C2HCl3, on clay mineral surfaces in the presence of water has been modeled as an example describing a general program that uses molecular dynamics simulations to study the sorption of organic materials at the clay mineral/aqueous solution interface. Surfaces of the clay minerals kaolinite and pyrophyllite were hydrated at different water levels corresponding to partial and complete monolayers of water. In agreement with experimental trends, water was found to outcompete C2HCl3 for clay surface sites. The simulations suggest that at least three distinct mechanisms coexist for C2HCl3 on clay minerals in the environment. The most stable interaction of C2HCl3 with clay surfaces is by full molecular contact, coplanar with the basal surface. This kind of interaction is suppressed by increasing water loads. A second less stable and more reversible interaction involves adsorption through single-atom contact between one Cl atom and the surface. In a third mechanism, adsorbed C2HCl3 never contacts the clay directly but sorbs onto the first water layer. To test the efficacy of existing force field parameters of organic compounds in solid state simulations, molecular dynamics simulations of several representative organic crystals were also performed and compared with the experimental crystal structures. These investigations show that, in general, in condensed-phase studies, parameter evaluations are realistic only when thermal motion effects are included in the simulations. For chlorohydrocarbons in particular, further explorations are needed of atomic point charge assignments. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 144–153, 1998