ALBERT

Notes: The structures and energies of B+13, observed experimentally to be an unusually abundant species among cationic boron clusters, have been studied systematically with B3LYP/6-31G* density functional theory. The most thermodynamically stable B+12 and B+13 clusters are confirmed to have planar or quasiplanar rather than globular structures. However, the computed dissociation energies of the 3-dimensional B+13 clusters are much closer to the experimental values than those of the planar or quasiplanar structures. Hence, planar and 3-dimensional B+13 may both exist.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 203-214, 1998

Notes: The structural determination of endohedral metallofullerenes has attracted special attention in disclosing the formation mechanism and developing new routes to bulk production. Recent advances in the theoretical and experimental studies are summarized with representative mono- and dimetallofullerenes such as M@C82 (M=Ca, Sc, Y, and La), Sc2@C84, and La2@C80.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 232-239, 1998

Notes: 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 Koliński 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

Notes: Vibrational analysis of tellurium tetrachloride, TeCl4, was performed with Hartree-Fock (HF), MP2, and generalized gradient approximation density functional theory (DFT) methods supplemented with polarized double-zeta split valence (DZVP) basis sets and relativistic effective core potentials (RECP) of Hay and Wadt. The molecular geometry is best reproduced at the HF and MP2/RECP+DZVP [polarized Hay and Wadt RECP for Te and 6-31G(d) basis set for Cl] levels of theory. The DFT methods gave rise to poorer results, especially those using Becke's 1988 exchange functional. Generally, the vibrational frequencies calculated by the MP2 and B3-type DFT methods with the all electron and RECP+DZVP basis sets as well as at the HF/RECP level were in satisfactory accord with the experimental data. The agreement was good enough to assist the assignment of the measured vibrational spectra. The best agreement with the experimental vibrational frequencies was achieved with the scaled HF/RECP force field. Consistent results were obtained for the unobserved A2 (ν4) fundamental, where the results of the best methods were within 4 cm-1. The best force fields were obtained with the following methods: Becke3-Lee-Yang-Parr and Becke3-Perdew/all electron basis, MP2 and Becke3-Perdew/RECP+DZVP, and HF/RECP. The methods using RECPs are advantageous for large-scale computations. The RECP basis set effectively compensates the errors of the HF method for TeCl4; however, it provides poor results with correlated methods.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 308-318, 1998

Notes: Based on a study of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) update Hessian formula to locate minima on a hypersurface potential energy, we present an updated Hessian formula to locate and optimize saddle points of any order that in some sense preserves the initial structure of the BFGS update formula. The performance and efficiency of this new formula is compared with the previous updated Hessian formulae such as the Powell and MSP ones. We conclude that the proposed update is quite competitive but no more efficient than the normal updates normally used in any optimization of saddle points.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 349-362, 1998

Notes: A Monte-Carlo model for the simulation of alkene hydrogenation on metallic catalysts has been developed and implemented in Fortran language. We describe the model employed for ethylene hydrogenation on platinum and show the flow chart of the program. Computational characteristics such as number of necessary calculations to reach steady state, running times on different platforms, and effect of the size of the catalyst matrix, are presented. Good correlation between simulated and experimental data was observed. A subroutine allows for visual observation of the reaction. This approach is very useful for obtaining a personal impression of the important factors governing the reaction. By using this example the advantages of Monte-Carlo simulation to test the level of understanding of catalytic phenomena are discussed.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 396-403, 1998

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

Notes: An approximation to the Fukui function in atoms recently proposed in the form of a gradient correction to the local density approximation expression is here investigated. The spatial behavior of this function is analyzed, focusing on the gradient correction term. Physical information on the shell structure of atoms is shown to be conveyed by the radial distribution of that term. The analytically modeled densities (AMD) procedure is also implemented, and global atomic hardnesses are calculated with Hartree-Fock and AMD representations of atomic electron densities.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 488-503, 1998

Notes: A new molecular mechanics force field has been developed that takes into account the fact that, upon coordination to a transition metal ion, the redistribution of electron density leads to small but significant structural changes in the organic backbone of the ligand. Structural studies indicate that the perturbation by coordination to a metal ion extends to the α-carbon atom of the donor, the perturbation is roughly independent of the metal center for M2+ and M3+ and negligible for M+, and the perturbation of the Cα(SINGLE BOND)Cα′ bond is roughly independent of the donor atom. New parameter sets for oxalates, imidazoles, and pyrazoles are also presented. The refined parameters have been validated with a large number of monodentate, multidentate, and macrocyclic ligands.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 512-523, 1998

Notes: A protein energy surface is constructed. Validation is through applications of global energy minimization to surface loops of protein crystal structures. For 9 of 10 predictions, the native backbone conformation is identified correctly. Electrostatic energy is modeled as a pairwise sum of interactions between anisotropic atomic charge densities. Model repulsion energy has a softness similar to that seen in ab initio data. Intrinsic torsional energy is modeled as a sum over pairs of adjacent torsion angles of 2-dimensional Fourier series. Hydrophobic energy is that of a hydration shell model. The remainder of hydration free energy is obtained as the energetic effect of a continuous dielectric medium. Parameters are adjusted to reproduce the following data: a complete set of ab initio energy surfaces, meaning one for each pair of adjacent torsion angles of each blocked amino acid; experimental crystal structures and sublimation energies for nine model compounds; ab initio energies over 1014 conformations of 15 small-molecule dimers; and experimental hydration free energies for 48 model compounds. All ab initio data is at the Hartree-Fock/6-31G* level.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 548-573, 1998

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

Notes: 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.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 610-627, 1998

Notes: A parallel full configuration interaction (FCI) code, implemented on a distributed memory MPP computer, has been modified in order to use a direct algorithm to compute the lists of mono- and biexcitations each time they are needed. We were able to perform FCI calculations on the ground state of the acetylene molecule with two different basis sets, corresponding to more than 2.5 and 5 billion Slater determinants, respectively. The calculations were performed on a Cray-T3D and a Cray-T3E, both machines having 128 processors. Performance and comparison between the two computers are reported and discussed.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 658-672, 1998

Notes: Hydrogen-saturated cut-outs of hexagonal boron nitride have been used to model the solid state. Model compounds have been geometry optimized by means of density functional theory, whereas chemical shift calculations have been carried out at the coupled-perturbed Hartree-Fock level of theory employing gauge-including atomic orbital (GIAO) basis sets. The reliability of results has been tested against experimental values for chemical shifts in stable molecules with similar structural elements. With increasing cluster size, viz. a vanishing influence of the saturating hydrogens on the innermost nitrogen atoms, we find a convergence of 15N chemical shifts. A classification scheme for the chemical environment of a nitrogen atom has been set up according to its bonding graph including the second coordination sphere. For a given connectivity, chemical shifts vary within a few parts per million, thus enabling us to predict a 15N-NMR chemical shift of -285 ± 5 ppm for solid α-boron nitride.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 716-725, 1998

Notes: A new electrostatic model for the calculation of infrared intensities in molecular mechanics and molecular dynamics is presented. The model is based on atomic charges, atomic charge fluxes, and internal coordinate dipoles and their fluxes. The internal coordinate dipoles are used instead of atomic dipoles, thus simplifying the derivation of parameters. The model is designed to reproduce ab initio dipole derivatives, and the parameters can be obtained by (iterative) transformations from these, or by linear least squares fitting to them. A first application to linear alkanes has been made. For these molecules, the intensities can be predicted with an average accuracy of 30-40%.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 754-768, 1998

Notes: The procedure previously proposed for the parameterization of the exchange-repulsion energy using probe atom calculations is applied here to the study of the methane dimer and refined to give a very accurate anisotropic model potential in terms of atomic parameters. The most accurate model uses sites on C and H atoms and requires 12 parameters, but a description using just four isotropic sites shifted inward from the H atoms by 15% of the bond length is almost as accurate and is very simple in form.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 847-857, 1998

Notes: The approach of Karfunkel and Gdanitz has been used to predict possible crystal structures of acetic acid and three of its monohalogenated analogs starting with the molecular structure alone. The results demonstrate that this approach is capable of finding many, if not all, of the possible packing arrangements of molecules of this size, but that it is not currently capable of correctly ranking these structures in terms of their enthalpy. This deficiency is probably due to inadequacies in the force field used to minimize the structures. The inadequacies relate to the description of acidic hydrogen bonds and halogen-halogen interactions.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1-20, 1998

Notes: 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 5β-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 Nα-(2-naphthyl-sulfonyl-glycyl)-D-para-amidino-phenyl-alanyl-piperidine (NAPAP) to human α-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 Å for the ligand atoms.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 21-37, 1998

Notes: A computer code was prepared for the molecular dynamics (MD) simulation of a multimolecular system to atomic resolution. Based on a widely accepted force field and aided by simulated annealing of single molecules, $\tilde V$ and $\tilde E$ were computed for ortho-terphenyl and 1,3,5-tri-α-naphthyl benzene across an extended range of temperatures. Although neither the simulation time (40-100 ps) nor the system size (27 and 64 molecules) appeared to affect the computational results to an appreciable extent, it was clear that a longer simulation time or a larger system tended to yield a more consistent set of data. In comparison to experimental observations, simulation was capable of representing $\tilde V$ to within 2-7%, Tg to within 10 K, and ΔCp across Tg to within 10%.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 86-93, 1998

Notes: We present an iterative constraint algorithm, QSHAKE, for use with semirigid molecules in molecular dynamics simulations. The algorithm is based on “SHAKE-ing” bond constraints between rigid bodies, whose equations of motion are solved in the quaternion framework. The algorithm is derived and its performance compared with SHAKE for liquid octane. QSHAKE is significantly more efficient whenever SHAKE requires triangles (or tetrahedra) of constraints to maintain molecular shape. Efficiencies are gained because QSHAKE reduces the number of holonomic constraints that must be solved iteratively and requires fewer iterations to obtain convergence. The gains in efficiency are most noticeable when a high degree of precision is imposed on the constraint criteria. QSHAKE is also stable at larger time steps than SHAKE, thus allowing for even more efficient exploration of phase space for semirigid molecules.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 102-111, 1998

Notes: Theoretical calculations in combination with experimental gas phase structure research can be performed in two ways. The first is to support and improve experimental analyses by including additional data from theoretical calculations. This is to the advantage of the experiment. The second way is a comparison of geometric structures and conformational properties obtained with different theoretical methods with the experimental result. This comparison indicates which theoretical method or methods are suitable for a specific compound. This approach is to the advantage of the theory.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 123-128, 1998

Notes: We present a formal and numerical analysis of the errors related to the use of molecular cavities in the continuum description of solute-solvent electrostatic interactions. In this approximation the solvent response is fully described by an apparent charge distribution induced on the cavity surface. The latter is then discretized into a set of point charges that are generally affected by errors of different origin but all depending on the definition of the cavity boundaries, and the way its surface is partitioned. The numerical analysis is based on results obtained for a set of couples of neutral/anionic solutes obtained with various versions of PCM methods, exploiting different correction procedures, as well as with two recently developed continuum methods, complementing PCM.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 833-846, 1998

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

Notes: Ab initio studies at the HF/6-31G* and B3LYP/6-31G* levels are reported for two bowl-shaped hydrocarbons related to C60: C30H12 and C36H12, of C3 and C3v symmetry, respectively. The former has an approximate heat of formation of 211 kcal/mol. Bowl-to-bowl interconversion may occur through a planar (C3h) form of ca. 64 kcal/mol greater energy having one imaginary vibrational frequency. The larger C36H12 bowl has a calculated ΔH°f of 265 kcal/mol. Its HF/6-31G*, B3LYP/6-31G*, and MM3 bond lengths are in good agreement with a recent X-ray structure. Chemical shifts for both compounds calculated by the GIAO method are in good agreement with the measured NMR spectra. The observed 13C chemical shifts increase with the extent of pyramidalization.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 189-194, 1998

Notes: The thermal isomerization of vinylcyclopropane to cyclopentene was discovered in 1960 and soon recognized as the simplest known example of a [1,3] sigmatropic shift. Experimental observations for the parent rearrangement and for isomerizations shown by substituted systems suggest that diradical transition structures are involved; recent theoretical treatments of the reaction find no minima corresponding to diradical intermediates. The common dichotomy opposing concerted versus diradical and thus necessarily stepwise mechanisms appears inappropriate. The reaction of vinylcyclopropane may involve four energetically concerted paths traversed by different conformational forms of nearly isoenergetic diradical species leading through four isometric diradical transition structures to cyclopentene.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 222-231, 1998

Notes: Electrostatic solvation free energies were computed for several small neutral bases and their conjugate acids using a continuum solvation model called the self-consistent isodensity polarizable continuum model (SCIPCM). The solvation energies were computed at the restricted Hartree-Fock (RHF) and second-order Møller-Plesset (MP2) levels of theory, as well as with the Becke3-Lee-Yang-Parr (B3LYP) density functional theory, using the standard 6-31G** Gaussian basis set. The RHF solvation energies are similar to those computed at the correlated MP2 and B3LYP theoretical levels. A model for computing protonation enthalpies for neutral bases in fluorosulfonic acid solvent leads to the equation ΔHprot, HSO3F(B)=-PA(B)+ΔEt(BH+)-ΔEt(B)+β, where PA(B) is the gas phase proton affinity for base B, ΔEt(BH+) is the SCIPCM solvation energy for the conjugate acid, and ΔEt(B) is the solvation energy for the base. A fit to experimental values of ΔHprot, HSO3F(B) for 10 neutral bases (H2O, MeOH, Me2O, H2S, MeSH, Me2S, NH3, MeNH2, Me2NH, and PH3) gives β=238.4±2.9 kcal/mol when ΔΔEt is computed using the 0.0004 e⋅bohr-3 isodensity surface for defining the solute cavity at the RHF/6-31G** level. The model predicts that for carbon monoxide ΔHprot, HSO3F(CO)=10 kcal/mol. Thus, protonation of CO is endothermic, and the conjugate acid HCO+ (formyl cation) behaves as a strong acid in fluorosulfonic acid.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 250-257, 1998

Notes: Ab initio molecular orbital theory and density functional theory calculations have been carried out on the reactions of the trifluoromethyl radical with the hydroxyl and the hydrogen radicals. These reactions are key reactions that underlie a new fire extinguishing mechanism of non-bromine-containing halon replacements. The activation energies calculated by the MP2 and QCISD methods are in good agreement with the experimental values. The B3LYP, as well as MP2 and QCISD, give good results for the calculations of the heats of reactions. The GAUSSIAN-1 and GAUSSIAN-2 theory calculations present the most acxcurate results on both the activation energies and the heats of reactions. The effects of the scaling factors on the activation energies and the heats of reactions are also evaluated.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 277-289, 1998

Notes: 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 α-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 Ooi et al. [Proc. Natl. Acad. Sci. USA, 84, 3086 (1987)] and from Wesson and Eisenberg [Prot. Sci., 1, 227 (1992)] showed a good correlation between solvation energies of the conformations and their root-mean-square deviations from the X-ray crystal structure of APP.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 319-333, 1998

Notes: 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.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 377-395, 1998

Notes: Possible crystal structures of acetic acid were generated, considering eight space groups and assuming one molecule in the asymmetric unit. Our grid-search method was compared with a Monte Carlo approach as implemented in the Biosym/MSI Polymorph Predictor. This revealed no sampling deficiencies. A large number of possible crystal structures were found (∼100 within only 5 kJ/mol), including the experimental structure. Energy minimizations were done with a united-atoms force field (GROMOS), an all-atoms force field (AMBER), and a potential that describes the electrostatic interactions with distributed multipoles (DMA). In all cases, the experimental structure had a low lattice energy. The number of hypothetical crystal structures was reduced considerably by removing space-group symmetry constraints, or by a primitive molecular dynamics shake-up. Nevertheless, sufficient structures of equal or lower energy compared with the experimental structure remained to suggest that other factors need to be considered for genuine structure prediction.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 459-474, 1998

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

Notes: The potential energy surface (PES) for the cyclooctane molecule was comprehensively investigated at the Hartree-Fock (HF) level of theory employing the 3-21G, 6-31G, and 6-31G* basis sets. Six distinct true minimum energy structures (named B, BB, BC, CROWN, TBC, and TCC1), characterized through harmonic frequency analysis, were located on the multidimensional PES. Two transition state structures were also located on the PES for the cyclooctane molecule. Electron correlation effects were accounted for using the Møller-Plesset second-order perturbation theory (MP2) approach. The predicted global minimum energy structure on the ab initio PES for the cyclooctane molecule is the BC conformer. A gas phase electron diffraction study at 300 K suggested a conformational mixture while an NMR study in solution at 161.5 K predicted the BC conformer as the predominant form. The equilibrium constants reported in the present study, which were evaluated from the ab initio calculated total Gibbs free energy change values, were in good agreement with both experimental investigations. The ab initio results showed that the low temperature condition significantly favored the BC conformer while above room temperature both BC and CROWN structures can coexist.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 524-534, 1998

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

Notes: We propose a stochastic optimization technique based on a generalized simulated annealing (GSA) method for mapping minima points of molecular conformational energy surfaces. The energy maps are obtained by coupling a classical molecular force field (THOR package) with a GSA procedure. Unlike the usual molecular dynamics (MD) method, the method proposed in this study is force independent; that is, we obtain the optimized conformation without calculating the force, and only potential energy is involved. Therefore, we do not need to know the conformational energy gradient to arrive at equilibrium conformations. Its utility in molecular mechanics is illustrated by applying it to examples of simple molecules (H2O and H2O3) and to polypeptides. The results obtained for H2O and H2O3 using Tsallis thermostatistics suggest that the GSA approach is faster than the other two conventional methods (Boltzmann and Cauchy machines). The results for polypeptides show that pentalanine does not form a stable α-helix structure, probably because the number of hydrogen bonds is insufficient to maintain the helical array. On the contrary, the icoalanine molecule forms an α-helix structure. We obtain this structure simulating all Φ, Ψ pairs using only a few steps, as compared with conventional methods.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 647-657, 1998

Notes: A method to construct the equivalent of multidimensional Ramachandran plots for nucleic acids on the basis of singular value decomposition (SVD) is presented. For this purpose, a data matrix containing 244 DNA dinucleoside monophosphate steps, represented by nine torsion angles, was decomposed into a score and loading matrix. It is shown that biplots, containing both score points and loading vectors, provide a simple tool to interpret the principles of DNA class separation. Scores separate the data matrix into one A-DNA class, two different B-DNA classes, and one so-called crankshaft class. Loading vectors correlate torsion angles. The projections of scores on loading vectors indicate which torsion angles play a dominant role in DNA class separation. The results of the biplots are supported by (simple) physical interpretations. From a three-dimensional score space the nine original torsion angles can be reconstructed. Hence, the potential to create the multidimensional equivalent of a Ramachandran plot is available; that is, forbidden and accessible regions in the reduced space reflect these same regions in the nine-dimensional original space.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 695-715, 1998

Notes: Interaction with the ligand binding domain of receptors for natural chemicals present one potential mechanism for the biological effects of environmental chemicals. Evidence suggests that the electrostatic interaction between the ligand and the receptor is an important component for binding to some of the relevant receptors. The presence of charged residues near the binding site suggests that the charge distribution of the free ligand may be different from the charge distribution of the ligand as it approaches the binding domain of the protein. In this study a new type of potential is computed for a series of dibenzo-p-dioxin (dioxin) ligands. This quantum mechanically computed potential results from interaction between the ligand and a trimethyl ammonium probe at a set of grid points. This interaction potential is compared with the molecular electrostatic potential computed from the wave function of the isolated ligands. Three types of local minima are found: (1) above the oxygen; (2) above the conjugated ring; and (3) above the chlorine(s). The molecular electrostatic potential emphasizes the minima associated with the chlorine atoms and, in that potential, the minima associated with the oxygen atoms disappear with chlorination. In the new potential, the minima over the oxygen atoms are maintained even in tetrachlorodioxin. As chlorination is increased the differences between the two potentials increases. The new potential shows the influence of the π-cation interaction, which is largest when there is little substitution on the ring. The presence of the probe induces a dipole component of 1 debye perpendicular to the plane of the ligand. Local minima in the interaction potential are then used as starting structures for the determination of the most stable ligand-probe complexes. The most stable structures are obtained from the minima associated with the oxygen atoms. These structures are stabilized by a hydrogen bond formation between the probe and the oxygen and the molecule is bent by 30° about the O(SINGLE BOND)O axis. For this series of molecules, the new potential retains some of the features that determine the hydrogen bond whereas the molecular electrostatic potential does not.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 673-684, 1998

Notes: The sensitivity of aqueous solvation free energies (SFEs), estimated using the GB/SA continuum solvent model, on charge sets, protocols, and force fields, was studied. Simple energy calculations using the GB/SA solvent model were performed on 11 monofunctional organic compounds. Results indicate that calculated SFEs are strongly dependent on the charge sets. Charges derived from electrostatic potential fitting to high level ab initio wave functions using the CHELPG procedure and “class IV” charges from AM1/CM1a or PM3/CM1p calculations yielded better results than the corresponding Mulliken charges. Calculated SFEs were similar to MC/FEP energies obtained in the presence of explicit TIP4P water. Further improvements were obtained by using GVB/6-31G** and MP2/6-31+G** (CHELPG) charge sets that included correlation effects. SFEs calculated using charge sets assigned by the OPLSA* force field gave the best results of all standard force fields (MM2*, MM3*, MMFF, AMBER*, and OPLSA*) implemented in MacroModel. Comparison of relative and absolute SFEs computed using either the GB/SA continuum model or MC/FEP calculations in the presence of explicit TIP4P water showed that, in general, relative SFEs can be estimated with greater accuracy. A second set of 20 mono- and difunctional molecules was also studied and relative SFEs estimated using energy minimization and thermodynamic cycle perturbation (TCP) protocols. SFEs calculated from TCP calculations using the GB/SA model were sensitive to bond lengths of dummy bonds (i.e., bonds involving dummy atoms). In such cases, keeping the bond lengths of dummy bonds close to the corresponding bond lengths of the starting structures improved the agreement of TCP-calculated SFEs with energy minimization results. Overall, these results indicate that GB/SA solvation free energy estimates from simple energy minimization calculations are of similar accuracy and value to those obtained using more elaborate TCP protocols.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 769-780, 1998

Notes: A general, fast, and exact optimization, called neighbor-list reduction (NLR), is presented, which can be used to accelerate the computation of hard-sphere molecular surface areas. NLR allows selected neighbors of a central atom to be removed from the computation in a preprocessing step, thus allowing the calculation of the atom's surface area to proceed with a shorter list of neighbors. The atoms removed are those having intersections with the central atom falling entirely within unions of other atoms' intersections with the central atom. We describe explicit methods for two levels of neighbor-list reduction: 3NLR considers three hard spheres at a time - the central atom, the candidate for removal, and one other neighbor; whereas 4NLR considers two other neighbors. We demonstrate the correctness and efficiency of this optimization by means of a modified version of the NACCESS program, which computes atomic and molecular surface areas numerically. As test cases we used compounds of different size and class, with and without explicit hydrogens. When van der Waals surface (vdWSA) is computed, the NLR methods reduce the length of the neighbor list by as much as 41%; when solvent-accessible surface area (SASA) is computed, the reduction is as great as 74%. The overall speed improvement due to these reductions is a factor of only about 1.2 for vdWSA, but is about 2.0 for the computation of SASA, in the context of this particular program. All 39,554 calculated atomic surface areas (vdWSA and SASA) were found to be identical to within 0.001 Å2 to those obtained without NLR.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 797-808, 1998

Notes: Quantitative assessment of the “best fit” between neutral molecules and the cavity of a “rigid” neutral receptor is a challenging task in supramolecular chemistry, drug design, and biology. We investigate this question by molecular dynamics and free-energy perturbation simulations performed on the macrocyclic ligand cryptophane-E (L) and its L·S complexes with three tetrahedral guests (S=CH2Cl2, CHCl3, and CCl4) in the gas phase and in chloroform solution. The van der Waals interactions are shown to play a crucial role in the calculated complexation selectivity. Calculations using Lennard-Jones 6-12 potentials and “standard” OPLS R*Cl and εCl parameters for the Cl atoms of S lead to a preference for CCl4, in contrast to the selectivities observed experimentally in solution (CHCl3 〉 CH2Cl2 〉 CCl4). Based on systematic investigations of the relative free energies of binding of CHCl3/S, we derive a set of R*Cl and εCl van der Waals parameters that account for experimental binding data. Although the complexes are of the van der Waals type, their electrostatic representation is also crucial for correct calculation of relative stabilities. Thus, the recognition of the “best guest” stems from a subtle balance of distance and time-dependent, cumulative noncovalent interactions between atoms of S and of L, which require an accurate representation. In addition, even in a weakly polar solvent, like chloroform, solvation effects are shown to modulate the recognition of the neutral substrates.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 820-832, 1998

Notes: The evaluation of the first-order scalar relativistic corrections to MP2 energy based on either direct perturbation theory or the mass-velocity and Darwin terms is discussed. In a basis set of Lévy-Leblond spinors the one- and two-electron matrix elements of the relativistic Hamiltonian can be decomposed into a nonrelativistic part and a relativistic perturbation. Thus, a program capable of calculating nonrelativistic energy gradients can be used to calculate the cross-term between relativity and correlation. The method has been applied to selected closed-shell atoms (He, Be, Ne, and Ar) and molecules (CuH, AgH, and AuH). The calculated equilibrium distances and harmonic frequencies were compared with results from first-order relativistic density functional calculations. It was found that the cross-term is not the origin of the nonadditivity of relativistic and correlation effects.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1596-1603, 1998

Notes: 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-Å-radius sphere and a 2.5-Å-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

Notes: Insertion of ethylene into the Ti-methyl bond in TiH2CH+3 is chosen as a model reaction for investigating the performance of a range of contemporary quantum chemical models in polymerization studies. Basis set effects are investigated at the self-consistent-field level, covering Hartree-Fock, pure DFT, and hybrid DFT. In agreement with findings in part I of this study, the basis set sensitivity of ethylene is shown to introduce a bias in computed energetics, amounting to 2-3 kcal/mol when DZP bases are used to compute the overall heat of monomer insertion. The geometry of stationary points relevant to the insertion reaction is determined using hybrid density functional theory. Based on these structures, the energy profile of the insertion reaction is computed using a range of popular quantum chemical approximations. The methods include Hartree-Fock and Møller-Plesset (MP) perturbation theory up through the fourth order in spin-restricted, spin-unrestricted, and spin-projected formalisms. Furthermore, configuration-interaction-based methods are included, of which the top level method is singly and doubly excited coupled clusters with a perturbative estimate of the contribution from triply excited configurations added [CCSD(T)]. The performance of the methods just mentioned, as well as three pure density functional and three hybrid density functional methods, are compared with respect to “best” relative energies, defined through extrapolation of CCSD(T) correlation energies according to the PCI scheme of Siegbahn and coworkers. Even though the MP series show poor convergence, spin-projected MP2, as well as two pure DFT methods (BPW91, BP86) and PCI-78 based on the MCPF method, show similar and very good agreement with best relative energies for the insertion reaction.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 947-960, 1998

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

Notes: We have developed a parallel version of our pseudospectral localized Møller-Plesset electronic structure code. We present timings for molecules up to 1010 basis functions and parallel speedup for molecules in the range of 260-658 basis functions. We demonstrate that the code is scalable; that is, a larger number of nodes can be efficiently utilized as the size of the molecule increases. By taking advantage of the available distributed memory and disk space of a scalable parallel computer, the parallel code can calculate LMP2 energies of molecules too large to be done on workstations.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1030-1038, 1998

Notes: We present a method to define van der Waals, solvent-accessible, and solvent-excluding molecular surfaces with their partition in nonoverlapping surface portions (tesserae). The procedure is more efficient than those available in the literature to describe solvent effects on molecular systems of large size, and it can also be applied to solutes of small size without reducing the accuracy of the output and without increasing computational times. All the tesserae are expressed in terms of spherical triangles, having all the characterizing elements (vertices, centers, etc.) analytically defined. The method was tested by comparing the results for the surface area and the solvation free energy (decomposed in electrostatic, dispersion, and steric contributions) obtained using the GEPOL procedure within the framework of the polarizable continuum model solvation method. These comparisons regard 87 molecules at the molecular mechanics level and 28 molecules at the ab initio Hartree-Fock level: the results are quite satisfactory.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1758-1776, 1998

Notes: A new algorithm is proposed for approximation to the molecular surface. It starts with a triangular mesh built on an ellipsoid embracing the whole molecular surface. The triangular mesh is obtained from an icosahedron subdivision sphere with highly uniform vertex distribution, and the embracing surface is deflated stepwise to the best adherence of its triangles onto the surface of the molecule. The deflating direction of each vertex of a triangle is defined by the vector normal at this point to the previous deflated embracing surface. Our results show that the speed of the triangulation embracing ellipsoid method and the quality of the surface obtained by the method are faster and better than the method that starts with a quadrilateral mesh built from meridian and parallel representations on an embracing sphere to get the molecular surface. Furthermore, the surface obtained by the method can be used directly to approximate the molecular surface by spherical harmonic expansions.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1805-1815, 1998

Notes: Molecular mechanics calculations were applied to the conformational analysis of two diasteroisomers, the pyrrolizidine alkaloids (PAs) retronecine and heliotridine. The application of reoptimized parameters for H bonding corrected the tendency of MM3(92) calculations to give unrealistic H(DOTTED BOND)O distances for intramolecular OH interactions occurring in both diasterisomers. Inversions in the H-bond direction of exo-retronecine and in the relative stability of heliotridine endo-exo conformers were also observed with the application of the new parameters. A set of probable conformers was obtained for each diasterisomer, based on conformational and Boltzmann population analysis. Only exo-puckered conformers were found in the retronecine set, whereas both exo- and endo-puckered conformers were obtained for heliotridine. Transition state conformations supplied arguments supporting the design of models for H-bond interconversion in the case of exo-retronecine and for the exo-endo interconversion of heliotridine. Reactivity behaviors and 1H-NMR data of both diasterisomers were elucidated in light of the theoretical results.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1853-1861, 1998

Notes: Ab initio geometry optimizations were performed on gaseous protonated glycine using the second-order Møller-Plesset perturbation theory with the 6-31G*, 6-31G**, 6-31+G**, and 6-311+G** basis sets. Eight energy minima and 12 saddle points in the low-energy region of the electronic potential energy surface were characterized. The global minimum was an amino N-protonated conformer containing an ionic H bond between the (SINGLE BOND)NH3+ and O(DOUBLE BOND)C(DIAGONAL BOND)(DIAGONAL BOND) groups. The lowest energy O-protonated conformer was stabilized by a conjugative attraction between the nitrogen lone-pair electrons and the positively charged planar fragment (SINGLE BOND)C(OH)2+. Relative electronic energies of the nine N- and 11 O-protonated species fall in the ranges of 0-10 and 30-40 kcal mol-1. At room temperature the equilibrium distribution contained the most stable N-protonated conformer almost exclusively. Additional subjects for investigation include the effects of basis set and electron correlation on the predicted structures, nonbonded interactions that influence the relative stability of protonated conformers, conformational interconversions based on intrinsic reaction coordinate calculations, and kinetic pathways for protonation and associated changes in Gibbs free energy. The work provides geometric, energetic, and thermodynamic data pertinent to the study of gas-phase ion chemistry of amino acids and peptides.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1862-1876, 1998

Notes: Self-consistent field calculations are done using two-body density functionals for the correlation energy. The corresponding functional derivatives are obtained and used in pseudo-eigenvalue equations analogous to the Kohn-Sham ones. The examples studied include atomic systems from He to Ar. The values obtained for ionization potentials, electron affinities, dipole polarizabilities, and virial ratios from these calculations are given, and the effect of exchange is addressed. The results obtained are in good agreement with experimental values, and are of the same quality as those given by accurate exchange-correlation functionals.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1887-1898, 1998

Notes: We have applied a discretized version of the generator coordinate Hartree-Fock method to generate adapted Gaussian basis sets for atoms Cs (Z=55) to Lr (Z=103). Our Hartree-Fock total energy results, for all atoms studied, are better than the corresponding Hartree-Fock energy results attained with previous Gaussian basis sets. For the atoms Cs to Lr we have obtained an energy value within the accuracy of 10-4 to 10-3 hartree when compared with the corresponding numerical Hartree-Fock total energy results.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 858-865, 1998

Notes: The traditional method of analyzing solution structuring properties of solutes using atom-atom radial distribution functions (rdfs) can give rise to misleading interpretations when the volume occupied by the solute is ignored. It is shown by using the examples of O(4) in α- and β-D-allose that a more reliable interpretation of rdfs can be obtained by normalising the rdf using the available volume, rather than the traditional volume of a spherical shell.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 363-367, 1998

Notes: We comment on the convergence of the general coupling operator for all types of one-configuration or multiconfigurational wave functions that still preserve the one-configuration structure for the energy expression. The choice on the best arbitrary real and antisymmetric parameters inherent in the coupling operator methodology is discussed, giving a theoretical reason. Another type of coupling operator is defined, presented, and analyzed. Finally, we give some numerical examples related to the low-lying electronic states of a cluster model for K2NiF4 solid.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 368-376, 1998

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

Notes: Electron correlation at the Møller-Plesset second-order level was incorporated into the π-system portion of MM3 calculations for several conformers of [10]annulene, [18]annulene, bicyclo[5.3.1]undecapentaene, and bicyclo[4.4.1]undecapentaene. The conformers with “localized” C(SINGLE BOND)C π bonds (strongly alternating bond lengths) were found to be of lower energy than their counterparts with “delocalized” C(SINGLE BOND)C π bonds (similar bond lengths) before correlation energy was included. Correlation always lowered the energies of the delocalized conformation more than it did that of the localized conformation, such that often the latter was found to be more stable after correlation energy was included in the calculation. When a delocalized structure was not at a stationary point on the MM3 energy surface, such comparison could not be made. An example is the porphin molecule.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 475-487, 1998

Notes: A simple, customizable connectivity scheme is rigorously defined in which pairs of atoms are classified into three categories. The tools of graph theory are used to analyze the molecular graph and to efficiently find rings and ring assemblies through a combination of pruning and homeomorphic reduction. The definition of natural internal coordinates is extended in a nonredundant fashion for the various cases of weakly interacting components and for fused ring systems. The ring system coordinates were tested and found to be superior to Z-matrix coordinates.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 504-511, 1998

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

Notes: Semiempirical (MNDO and PM3) molecular orbital calculations have been undertaken to study the structures of the ground and excited states of 2,5-distrylpyrazine dye to assess its activity as a laser dye. In the ground and first excited singlet states, the trans-trans structure of C2h symmetry is the most stable structure in the gas phase and in DMSO, which agrees with the experimental findings. Upon excitation, the flexibility of the molecule decreases, leading to a subsequent decrease in the radiationless deactivation pathway and this increases the fluorescence efficiency of DSP. The absorption, excitation, and emission spectra have been calculated at the MNDO level using the PM3 optimized geometries in DMSO. At this level the agreement between theory and experiment is quite good. An estimated absorption band at 377 nm (expt 380 nm) is assigned to the S0→S1 transition. The excited state absorption band at 457 nm (expt 460 nm) is assigned to the S1→S12 transition. The emission band at 458 nm (expt 460 nm) is assigned to the S′1→S′0 transition. The overlap between the emission and the excited-state absorption spectra is presumably the main reason behind the reduced laser activity of the investigated dye.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 585-592, 1998

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

Notes: Computational studies of the minimum energy pathway for internal rotation of a methyl group are often made by constraining one dihedral angle at a sequence of values and optimizing all other parameters. When this is done, the methyl group adopts an asymmetric configuration at intermediate values of the torsion angle, with unequal bond lengths, bond angles, and torsion angles, even though the moiety against which it is rotating is another methyl group. The potential surface leading to this phenomenon is investigated using Hartree-Fock SCF calculations at the 6-31G* and 6-311G** levels and the detailed structural behavior of the methyl group during the course of the internal rotation is examined. It is shown that the nature of the constraint governs the resulting deformation of the methyl group symmetry.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1141-1145, 1998

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

Notes: A parallel implementation of the internally contracted (IC) multireference configuration (MRCI) module of the MOLPRO quantum chemistry program is described. The global array (GA) toolkit has been used in order to map an existing disk-paging small-memory algorithm onto a massively parallel supercomputer, where disk storage is replaced by the combined memory of all processors. This model has enabled a rather complicated code to be ported to the parallel environment without the need for the wholesale redesign of algorithms and data structures. Examples show that the parallel ICMRCI program can deliver results in a fraction of the time needed for equivalent uncontracted MRCI computations. Further examples demonstrate that ICMRCI computations with up to 107 variational parameters, and equivalent to uncontracted MRCI with 109 configurations, are feasible. The largest calculation demonstrates a parallel efficiency of about 80% on 128 nodes of a Cray T3E-300.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1215-1228, 1998

Notes: A method is presented to generate and triangulate molecular surfaces rapidly. It is based on the ‘marching tetrahedra’ approach. The method is fast, simple and easy to implement. Our approach is not analytical in nature. Hence no special treatment is required for complications with singularity, degeneracy, or with self-intersecting re-entrant surfaces. A quick test for determining the solvent accessibility of a point in space forms an important part of the method. This test has potential use outside of the surface generation algorithm such as in molecular field analysis where the solvent accessibility of a point needs to be determined. The triangulated surface generated is suitable for molecular graphics display as well as boundary element continuum dielectric calculations.   © 1998, Government of Canada. Exclusive worldwide publication rights in the article have been transferred to John Wiley & Sons, Inc. in perpetuity.   J Comput Chem 19: 1268-1277, 1998

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

Notes: The structure of lithium and sodium cyanates and isocyanates and their related ion pair SN2 reactions were investigated using molecular quantum mechanics at the Hartree-Fock (HF)/6-31G**//HF/6-31G** level. Extensive further calculations of some of the lithium systems at higher levels gave no significant changes in relative properties. The isocyanate ion pairs are the most stable monomeric forms. The lowest energy dimers are planar eight-membered rings. For the ionic SN2 reaction of cyanate ion and methyl fluoride and chloride, methyl isocyanate is the predicted major product. The monomer ion pair inversion mechanism is predicted to preferentially form methyl cyanate but the calculations also indicate that reaction with dimeric ion pairs have no clear preference for forming methyl cyanate or isocyanate.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1325-1336, 1998

Notes: Four alternatives are compared for estimating vibrational anharmonicity constants without explicitly calculating the expensive fourth derivatives of the potential curves. In the first, semiempirical approach, fourth derivatives for 53 diatomic molecules are estimated from ab initio second and third derivatives by using the Morse model potential. Vibrational anharmonicities ωexe are then computed from the third and fourth derivatives. The second approach invokes a purely empirical linear correlation between ωexe and the harmonic frequencies ωe. The third and fourth empirical approaches suppose that the effective harmonic and anharmonic force constants are proportional (with an additive constant in the fourth approach). Experimental values for ωexe are compared with empirical predictions and with semiempirical estimates based upon Hartree-Fock (HF), Møller-Plesset (MP2), and local, nonlocal, and hybrid density-functional theories (DFT), using the small 6-31G* basis set. Ab initio values of ωe and bond lengths re are also compared against experiment. The (U)MP2 results are the worst and include several anomalies. The other semiempirical methods yield results of comparable accuracy for ωexe of hydrides, although the DFT methods are markedly better for ωe and re and for ωexe of nonhydrides. The empirical estimates are nearly as good as the semiempirical ones. We conclude that: (1) both empirical and semiempirical approximations are useful for predicting stretching anharmonicity constants ωexe to precisions of σ≈5 cm-1 for hydrides and σ≈1.5 cm-1 for nonhydrides; and (2) MP2 theory is relatively unreliable for such calculations. In addition, we find the following tests to be useful when evaluating the reliability of vibrational constants calculated at the UMP2 level: (a) the calculated values of ωe and ωexe should not deviate substantially from the empirical relations; (b) harmonic frequencies and intensities calculated at the MP2 level should be smaller than those calculated at the corresponding HF level; (c) a large distance-dependence of the spin contamination, d〈S2〉/dR≳0.05 Å-1, suggests that calculated constants are too large.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1315-1324, 1998

Notes: A newly developed acyl halide molecular mechanics (MM3) force field can accurately calculate molecular geometry, usually to within experimental error. The new force field can also calculate vibrational frequencies. The acyl halides studied were formyl halides, acetyl halides, propionyl halides, n-butyryl halides, 2-methylpropionyl halides, and 2,2-dimethylpropionyl halides. The rms deviation for vibrational frequencies was 28 cm-1.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1387-1401, 1998

Notes: A completely general two-dimensional (2D) methodology for the classical simulation of reactive and nonreactive events on ab initio potential energy surfaces is introduced and tested. The methodology requires the minimum amount of information given a priori - geometries and energies at these geometries. From a list of ab initio geometries and energies, simulations may be executed and a distribution of outcomes obtained. The method introduced attempts a local approach at simulating the dynamics of the system, rather than a global analytic fit to the potential energy surface.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1431-1444, 1998

Notes: Atomic radii used to define the solute cavity in continuum-based methods are determined by reproducing the solvent-accessible surface defined as the loci of minima in a potential (solvent interaction potential) between the solute and a probe. This potential includes electrostatic interaction (ion-dipole, ion-quadrupole, and ion-induced dipole) terms as well as a Lennard-Jones energy term. The method alleviates the need to distinguish solute atoms in different chemical environments. These radii, when used in the calculation of solvation free energies, are shown to be superior to fixed atom-specific radii or to radii obtained from the electron isodensity surface from quantum-mechanical calculations.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1482-1493, 1998

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

Notes: We evaluate several multiple time step (MTS) molecular dynamics (MD) methods with respect to their suitability for protein dynamics simulations. In contrast to the usual check of conservation of total energy or comparisons of trajectory details, we chose a problem-oriented approach and selected a set of relevant observables computed from extended test simulations. We define relevance of observables with respect to their role in the description of protein function. Accordingly, the use of quantities that exhibit chaotic behavior, like trajectory details, is shown to be inappropriate for the sake of the evaluation of methods. The accuracy of a cutoff method and of six MTS methods is evaluated, which differ in their treatment of the computationally crucial long-ranged Coulomb interactions. For each of the observables considered, the size of purely statistical fluctuations is determined to allow identification of algorithmic artifacts. The obtained ranking of the considered MD methods differs significantly from that obtained by the usual measures of algorithmic accuracy. One particular distance class method, DC-1d, is shown to be clearly superior in that no algorithmic artifacts were detected.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1534-1552, 1998

Notes: A new molecular descriptor of hydrophobicity based on molecular quantum similarity measures (MQSM), which can be used to replace the log P parameter in QSAR studies, is proposed. Unlike the majority of existing approaches for calculation of log P, the present methodology does not rely on the use of fragment additive contributions, but rather it is based on the comparison of quantum chemically calculated electron density distributions of a given molecule in water and in 1-octanol, using MQSM. The method has been tested on a broad series of 58 molecules including such structural types as aliphatic hydrocarbons, alcohols, amines, halides, carboxylic acids, esters, amides, and ketones, as well as more complex systems with two functional groups. In all cases investigated, an excellent linear relationship between calculated MQSM and log P values was found. Additionally, an example of QSAR analysis is presented using MQSM instead of log P values, corresponding to predict the narcosis of tadpoles.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1575-1583, 1998

Notes: The equilibrium geometries and relative stabilities of several structural isomers of tungsten hexahydride, WH6, have been obtained at different levels of quantum chemical calculations. The performance of various strategies to (i) include electron correlation, viz. density functional theory based approaches, Møller/Plesset perturbation and coupled cluster theory, and to (ii) account for scalar relativistic effects, viz. various relativistic effective core potentials, first order perturbation theory, a quasi-relativistic treatment employing a Pauli Hamiltonian, and use of the Douglas/Kroll operator, are compared to the best theoretical data available. It is shown that relativistic and electron correlation effects are most important for the high-symmetry species, that these effects give rise to opposite trends in relative energies, and that overall the relativistic effects dominate. The most efficient way to incorporate relativistic effects appears to be via the use of relativistic effective core potentials, while the correlation energies are best taken account of using a conventional method such as CCSD(T).   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1604-1611, 1998

Notes: We introduce error weighting functions into the perturbative Monte Carlo method for use with a hybrid ab initio quantum mechanics/molecular mechanics (QM/MM) potential. The perturbative Monte Carlo approach introduced earlier provides a means to reduce the number of full SCF calculations in simulations using a QM/MM potential by evoking perturbation theory to calculate energy changes due to displacements of an MM molecule. The use of weighting functions, introduced here, allows an optimal number of MM molecule displacements to occur between the performance of the full self-consistent field calculations. This will allow the ab initio QM/MM approach to be applied to systems that require more accurate treatment of the QM and/or MM regions.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1632-1638, 1998

Notes: A new method for quantitatively comparing calculated vibrational modes is described that relies on projecting the vectors describing the normal modes of one molecule onto those of a basis molecule. The procedure virtually automates the assignment of vibrational modes from one molecule to a second, structurally similar one. We illustrate the concept by using the classical Wilson modes of benzene as a basis for describing normal modes of the monosubstituted benzene derivatives phenol, phenol-d5, and phenol radical cation. These examples demonstrate the method's power for accurately assigning and comparing the normal modes of molecules perturbed by chemical substitution, isotopic substitution, or oxidation state. The vibrational projection analysis method - a special case of vector projection analysis - is compared and contrasted with total energy distribution analysis, perhaps the most commonly used technique for quantitatively comparing vibrational modes, and is found superior in each case when comparing normal modes. Vibrational projection analysis need not be limited to single molecules and quantum calculations, because normal modes may be obtained for much larger systems using molecular mechanics or molecular dynamics techniques. The method should therefore prove useful for interpreting the normal modes of ordered periodic solids and structures perturbed by noncovalent contacts, including proteins and polymers. The method may also prove useful in evaluating new computational methods by allowing direct, quantitative comparison of the vibrational modes obtained from different techniques. The power of this technique will make vibrational projection analysis a valuable tool for computational chemists as well as those using calculations to complement vibrational spectroscopic measurements.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1663-1674, 1998

Notes: A parallel algorithm for efficient calculation of the second derivatives (Hessian) of the conformational energy in internal coordinates is proposed. This parallel algorithm is based on the master/slave model. A master processor distributes the calculations of components of the Hessian to one or more slave processors that, after finishing their calculations, send the results to the master processor that assembles all the components of the Hessian. Our previously developed molecular analysis system for conformational energy optimization, normal mode analysis, and Monte Carlo simulation for internal coordinates is extended to use this parallel algorithm for Hessian calculation on a massively parallel computer. The implementation of our algorithm uses the message passing interface and works effectively on both distributed-memory parallel computers and shared-memory parallel computers. We applied this system to the Newton-Raphson energy optimization of the structures of glutaminyl transfer RNA (Gln-tRNA) with 74 nucleotides and glutaminyl-tRNA synthetase (GlnRS) with 540 residues to analyze the performance of our system. The parallel speedups for the Hessian calculation were 6.8 for Gln-tRNA with 24 processors and 11.2 for GlnRS with 54 processors. The parallel speedups for the Newton-Raphson optimization were 6.3 for Gln-tRNA with 30 processors and 12.0 for GlnRS with 62 processors.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1716-1723, 1998

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

Notes: A comparative study was carried out to test the efficiency with which Metropolis Monte Carlo (MC) and stochastic dynamics (SD) sample the potential energy surface of the N-acetyl glycyl glycine methylamide peptide as defined by the united atom AMBER* force field. Boltzmann-weighted ensembles were generated with variations of all internal degrees of freedom (i.e., stretch, bend, and torsion) for a single N-acetyl glycyl glycine methylamide molecule at 300 K by 108-step MC and 100-ns SD simulations. As expected, both methods gave the same final energetic results. However, convergence was found to be ∼10 times faster with MC than with SD as measured by comparisons of the populations of all symmetrically equivalent conformers.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1294-1299, 1998

Notes: The concepts of averaged electron densities and averaged nuclear arrangements are discussed with reference to small nuclear distortions and an electron density transformation based on the inverse Löwdin transform.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1337-1344, 1998

Notes: High-level ab initio and DFT molecular orbital calculations have been used to investigate the physical properties of a model low-barrier hydrogen bond (LBHB) system: formic acid-formate anion. In the gas phase, it is found that the hydrogen bond formed is extraordinarily short and strong [ca. 27 kcal/mol at B3LYP/6-31++G(d, p)], with a calculated enthalpy of activation for proton transfer from donor to acceptor that is less than the zero-point vibrational energy available to the system. Several perturbations to this system were studied. Forcing a mismatch of pKas between donor and acceptor, via the use of substituents, causes the strength of the hydrogen bond to decrease. Microsolvation of the hydrogen-bonded complex does not affect the strength of the low-barrier hydrogen bond very much. Small variations in the structure of the LBHB results in a decrease in hydrogen-bond strength. Increasing the effective polarity of the cavity surrounding the LBHB was found to have a significant impact on the strength of the hydrogen bond. Implications for enzyme catalysis are discussed.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1345-1352, 1998

Notes: The relative reactivity of a series of nucleophiles that includes ethylene, sulfides, sulfoxides, amines, and phosphines toward dioxirane, dimethyldioxirane, carbonyloxide and dimethylcarbonyloxide has been examined at the MP4/6-31G*//MP2/6-31G*, QCISD(T)/6-31G*//MP2/6-31G*, and B3-LYP/6-31G* levels of theory. The barriers for the oxidations with dimethyldioxirane are higher (up to 2.5 kcal/mol for the oxidation of H2S) than those for the oxidations with the parent dioxirane. The oxidation barriers for dioxirane are larger than those for the oxidations with peroxyformic acid, except the barriers for the oxidation of sulfoxides. The reactivity of dimethylsulfide toward dimethyldioxirane was found to be comparable to that of dimethylsulfoxide both in the gas phase and in solution (chloroform). The classical gas phase barrier for the oxidation of trimethylamine to trimethylamine oxide was higher (6.3 kcal/mol at the MP4//MP2/6-31G* level) than that for oxygen atom transfer to trimethylphosphine. When the transition states were examined by self-consistent reaction field (SCRF) methods, the predicted barriers for the oxidation of amines and phosphines were found to be in good agreement with experiment. The general trend in reactivity for oxidation by dioxirane was R2S≈R2SO, R3P〉R3N in the gas phase, and R2S≈R2SO, R3N≈R3P (R=Me) in solution. The oxidation barriers calculated using the B3-LYP functional were lower than those computed at the MP4 and QCISD(T) levels.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1353-1369, 1998

Notes: Heats of formation of alkanes have been calculated with an accuracy of better than 0.36 kcal/mol by using the total energy calculated by density functional theory, plus bond and group equivalents and statistical mechanical corrections. The necessary equivalents were assigned to bonds and groups in molecules. Once such equivalents have been derived from the fit to available experimental values for a large and diverse set of compounds, they can be used to predict heats of formation for compounds of the same class for which these quantities are not experimentally available. Expanding the method to a new class of compounds requires that only new groups of equivalents for that class be added to the scheme. This provides a path for the systematic expansion of the model to new classes of compounds, and gives us a computational method for getting around the lack of experimental information about systems of interest.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1421-1430, 1998

Notes: The PESP (Parameterized ElectroStatic Potential) method for calculating molecular electrostatic potentials, previously parameterized for H, C, N, O, F, P, S, Cl, and Br, is extended to molecules containing Li+, Na+, Mg2+, K+, Ca2+, Zn2+, and I. For a collection of 166 molecules containing 1668 atoms with at least one metal or iodine atom, PESP achieves an average absolute deviation in electrostatic potential-derived atomic charges of 0.042e- compared with ab initio MP2/6-31G** calculations, with a correlation coefficient of 0.996. For a larger data set, consisting of 311 molecules encompassing all of the 16 elements just listed (2488 total atoms), PESP achieves an average absolute deviation of 0.040e- and a correlation coefficient of 0.995. PESP calculations are an order of magnitude faster than the simplest ab initio method (STO-3G) on large molecules, while achieving a level of accuracy that rivals much more elaborate ab initio methods.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1456-1469, 1998

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

Notes: No abstract.

Notes: The potential energy hypersurfaces for the triple inversion, from chair to boat and α to β conformations, are explored theoretically in 3-azabicyclo[3.3.1]nonan-9-one and its N-methyl derivative, by using ab initio quantum-mechanical calculations. Both compounds are precursors of rigid analogs of the potential GABAA and GABAB receptor antagonists. In contrast to results from semiempirical calculations, the chair-chair β conformers are found to be, by far, the most stable structures for both the nonmethylated and N-methylated compounds. The inversion barriers are found to be relatively low, so that the conformers could be expected to exist in thermodynamic equilibrium at room temperature. A population analysis reveals, however, that, in the ab initio approach, the molecules seem to exist practically only in the chair-chair-β conformation. The theoretical results compare well with the available experimental data.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1567-1574, 1998

Notes: When canonical molecular orbitals are expanded in terms of a set of localized molecular orbital building blocks, called bond orbitals, the character of the canonical molecular orbitals can be characterized according to the component bond orbitals resembling the core, lone pair, and localized bond building blocks in an intuitive Lewis structure. Weinhold's natural bond orbital method can produce a unique Lewis structure with total occupancy of its occupied bond orbitals exceeding 99.9% of the total electron density for simple molecules. Two useful indices, Lewis bond order and weight of lone pair orbitals, can be defined according to the weights of the bonding and lone pair components of this unique Lewis structure. Calculation results for molecules N2, CO, CS, NO, HCN, C2H2, H2O, and H2S show that the former index can account for the vibrational structures of photoelectron spectroscopy, whereas the latter index can account for the band intensity enhancement of Penning ionization electron spectroscopy.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 882-892, 1998

Notes: A detailed investigation of the reaction path for the thermal rearrangement of 3,4-dihydro-1aH-azirine[2,3-c]pyrrol-2-one to yield a cyanoketene-formaldimine complex is carried out at the MP2/6-31G* and B3LYP/6-31G* levels of theory. The ring opening of the five-membered pyrrolinone ring and the formation of the nitrile group takes place in a concerted manner, presenting a significant strain energy release and allowing for an electronic stabilization by coarctate conjugation of the transition structure (TS). These two factors make possible a moderate energy barrier. Although the structural features B3LYP/6-31G* theoretical levels, it is found that the MP2 energy barrier (28.8) CCSD(T)/6-31G*//MP2/6-31G* value (17.1 kcal/mol). The complex electronic rearrangement can be rationalized using the theory of coarctate transition structures developed by Herges as the evolution of an azirine structure without referring to a hypothetical vinyl nitrene intermediate.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 912-922, 1998

Notes: We have carried out theoretical calculations to analyze molecular interactions and proton transfer mechanisms in the formate-imidazole-water system, which may be considered the simplest model of catalytic triads in serine proteases. Computations were carried out at the density functional theory level. The effect of a dielectric environment on energy surfaces is considered using a polarizable continuum model and the self-consistent reaction field approach. The role played by inertial and noninertial polarization of this environment is emphasized. Nonequilibrium solvation effects have been estimated. The results show that there are different reaction mechanisms, concerted or stepwise, that may be competitive, depending on the nature of the molecular environment.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1675-1688, 1998

Notes: Asymmetric donor-acceptor-substituted π-conjugated systems with low dipole moments and structural components that favor parallel alignment of neighboring molecules are potential molecular organic materials with nonlinear optical properties for which near prefect dipole parallel alignment is possible. The asymmetrical 4-methoxysubstituted acetophenone azines with the substituents fluorine (1), chlorine (2), bromine (3), cyano (4), and nitro (5) in the 4′-position have been studied in this context, and for 2 and 3 the dipole parallel alignment has indeed been accomplished in crystals of the pure material. In the present study, the effects of asymmetrization on the structure and the electronic structures of the ground states of 1-5 have been explored at the RHF/6-31G* level. The properties of the optimized structures of the azines 1-5, a comparative analysis of asymmetrical and symmetrical azines, and natural population analyses, all show no significant evidence for conjugation over the azine bridge. The concept of azine spacers as “conjugation stoppers” is useful and the dipole moments of asymmetrical azines are relatively small. The analysis suggests that the charge transfer between the donor and acceptor contributes to the dipole moment much less than is generally assumed.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1130-1140, 1998

Notes: Transition metal complexes have been treated by a molecular mechanics approach using an MM2-like valence force field for the ligands and a Urey-Bradley force field about the metal. Metal-ligand atom bonds are parameterized with dependence on the metal covalent radius. The current metal-ligand distances and covalent radii were determined by a Simplex minimization and visual optimization of the deviations in bond distance between 230 structures from the Cambridge Crystallographic Structure Database and the force field calculated values. Moderately good reproducibility of structures is obtained with a variety of metal types with different oxidation states, and degree of unsaturation and with square-planar compounds. The current model also includes coordination of π ligands to the metal center with a 1/r6 attractive term whose minimum is 25 kcal/mol per p orbital at a distance dependent on the covalent radius of the metal and an angular attentuation factor.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1167-1178, 1998

Notes: The spectroscopic properties of XSiO (X=F, Cl, or Br) have been predicted using the B3-LYP/6-311+G(2d) level of theory. It has been shown that the halogen atom is Si bonded in a bent structure, with ∠(XSiO) bond angles close to 126°. The binding energy of the halogen with the SiO subunit was calculated to be -80.1, -40.9, and -29.3 kcal/mol for FSiO, ClSiO, and BrSiO, respectively. The harmonic frequencies and isotopic shifts have been calculated. A comparison between XSiO and X2SiO has also been made. For the X2SiO (X=F or Cl) compounds, the calculated frequencies are in fair agreement with the available experimental data. Characterization of bonding has been investigated with different approaches (natural bond orbital approach, topological analysis of the charge density, and of the electron localization function ELF).   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1205-1214, 1998

Notes: An algorithm is presented for the four-index transformation of electron repulsion integrals to a localized molecular orbital (MO) basis. Unlike in most programs, the first two indices are transformed in a single step. This and the localization of the orbitals allows the efficient neglect of small contributions at several points in the algorithm, leading to significant time savings. Thresholds are applied to the following quantities: distant orbital pairs, the virtual space before and after the orthogonalizing projection to the occupied space, and small contributions in the transformation. A series of calculations on medium-sized molecules has been used to determine appropriate thresholds that keep the truncation errors small (below 0.01% of the correlation energy in most cases). Benchmarks for local second-order Møller-Plesset perturbation theory (MP2; i.e., MP2 with a localized MO basis in the occupied subspace) are presented for several large molecules with no symmetry, up to 975 contracted basis functions, and 60 atoms. These are among the largest MP2 calculations performed on a single processor. The computational time (with constant basis set) scales with a somewhat lower than cubic power of the molecular size, and the memory demand is moderate even for large molecules, making calculations that require a supercomputer for the traditional MP2 feasible on workstations.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1241-1254, 1998

Notes: When calculating free energy differences between two molecular systems by means of molecular dynamics simulation, accessory potential functions can help eliminate uninteresting configurational entropy contributions, improve convergence, and facilitate reversibility. In this work, we demonstrate that the use of a harmonic potential function to restrain key portions of a molecular system in a free energy perturbation dual-topology molecular dynamics approach dramatically improves convergence and precision of the calculation. Limitations of this technique are illustrated, and its use in conjunction with a fixed bond-length constraint is developed.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1278-1283, 1998

Notes: We present a method of fitting arbitrary functions to linear combinations of Gaussians. In particular, we discuss an adaptation of Prony's method, or separation of exponentials, which allows us to automatically select appropriate exponents for these Gaussians. We then apply this technique to the selection of dealiasing sets for pseudospectral electron correlation methods. We show that it can successfully choose functions that generally improve the accuracy of pseudospectral correlation energies while reducing the size of the dealiasing set chosen.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1300-1314, 1998

Notes: A new molecular mechanics (MM3) force field has been developed based on various experimental data as well as ab initio calculations. Computer-generated molecular structures and energy values were compared with experimentally determined data. The acyl halides studied were formyl halides, acetyl halides, propionyl halides, n-butyryl halides, 2-methylpropionyl halides, and 2,2-dimethylpropionyl halides. The rms deviations were 0.005 Å and 1.06° for bond lengths and bond angles, respectively. MM3 was in good overall agreement with the available structural, conformational, and thermodynamic data.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1370-1386, 1998

Notes: We report ab initio [MP2(full)/6-31G*] and force field (MMP2 extended to carbocations) results on the parent and alkylated singlet cyclopentadienyl cations. The ab initio results are complemented by force field calculations of 1,3 π-interactions and antiaromatic destabilization energies. The antiaromatic destabilization of the singlet cyclopentadienyl cation is larger than that of cyclobutadiene. The MMP2 heats of formation of various cyclopentadienyl cations agree with experimental and ab initio data. Two Jahn-Teller distorted isomeric structures with almost identical energies are found for the parent singlet cyclopentadienyl cation. Dynamic and nondynamic correlation calculations employing the CASSCF(4,5)/6-31G*//MP2(full)/6-31G* and CASPT2(4,5)/6-31G*//MP2(full)/6-31G* levels of theory clearly show the energetic favorability of an allylic structure for the parent singlet cyclopentadienyl cation. Aklylated cyclopentadienyl cations prefer one of the two possible isomeric cyclopentadienyl cation forms. The substitution pattern determines the preferred cyclopentadienyl structure. As a consequence of its strong antiaromatic destabilization the singlet t-butylcyclopentadienyl cation is predicted to rearrange to a homoallylic 2-(1-methyl-2,4-cyclopentadienyl)-2-propyl cation.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1402-1420, 1998

Notes: A general method designed to isolate the global minimum of a multidimensional objective function with multiple minima is presented. The algorithm exploits an integral “coarse-graining” transformation of the objective function, U, into a smoothed function with few minima. When the coarse-graining is defined over a cubic neighborhood of length scale ∊, the exact gradient of the smoothed function, U∊, is a simple three-point finite difference of U. When ∊ is very large, the gradient of U∊ appears to be a “bad derivative” of U. Because the gradient of U∊ is a simple function of U, minimization on the smoothed surface requires no explicit calculation or differentiation of U∊. The minimization method is “derivative-free” and may be applied to optimization problems involving functions that are not smooth or differentiable. Generalization to functions in high-dimensional space is straightforward. In the context of molecular conformational optimization, the method may be used to minimize the potential energy or, preferably, to maximize the Boltzmann probability function. The algorithm is applied to conformational optimization of a model potential, Lennard-Jones atomic clusters, and a tetrapeptide.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1445-1455, 1998

Notes: A new and efficient method for overcoming the multiple minima problem of polypeptides, the systematic stepsize variation (SSV) method, is presented. The SSV is based on the assumption that energy barriers can be passed over by sufficiently large rotations about rotatable bonds: randomly chosen dihedral angles are updated starting with a small stepsize (i.e., magnitude of rotation). A new structure is accepted only if it possesses a lower energy than the precedent one. Local minima are passed over by increasing the stepsize systematically. When no new structures are found any longer, the simulation is continued with the starting structure, but other trajectories will be followed due to the random order in updating the torsional angles. First, the method is tested with Met-enkephalin, a peptide with a known global minimum structure; in all runs the latter is found at least once. The global minimum conformations obtained in the simulations show deviations of ±0.0004 kcal/mol from the reference structure and, consequently, are perfectly superposable. For comparison, Metropolis Monte Carlo simulated annealing (MMC-SA) is performed. To estimate the efficiency of the algorithm depending on the complexity of the optimization problem, homopolymers of Ala and Gly of different lengths are simulated, with both the SSV and the MMC-SA method. The comparative simulations clearly reveal the higher efficiency of SSV compared with MMC-SA.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1470-1481, 1998

Notes: The modeling of voltage-gated ion-channel proteins is a continuing challenge for force-field calculations because of the diverse range of interactions involved. In particular, current force fields are not parameterized for either ion-amino acid or amino acid-electric field interactions. To address the parameterization of ion-amino acid interactions, we have tested the use of empirical correction terms, derived from ab initio calculations of single amino acids (representing the peptide backbone) interacting with K+ ions. Having demonstrated the utility of such an approach, we then extended the application to the amino acid side chains. The calculation of the interaction of K+ with serine, cysteine, methionine, lysine, arginine, aspartate, histidine (uncharged), tyrosine, tryptophan, and phenylalanine, both completes the parameterization of the molecular environments contained in the amino acids, and allows specific comment on these ion-functional group interactions. The cation-π interactions were of particular interest, given recent proposals in the literature and the fear that force fields would not be able to treat such interactions. We present a comprehensive comparison of the ab initio (DFT [BLYP], 6-311 G**) and force field (CHARMm22.0) assessments of these interactions.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1515-1525, 1998

Notes: We have implemented geometry optimization using an analytic gradient to a two-component Kramers' restricted Hartree-Fock (KRHF) method for polyatomic molecules with closed-shell configurations. The KRHF method is a Hartree-Fock method based on relativistic effective core potentials with effective spin-orbit operators. The derivatives of spin-orbit integrals are obtained by numerical differentiation. Geometries for the various forms of polyatomic hydrides containing row 6 p-block elements are optimized with and without spin-orbit interactions. The structural changes due to spin-orbit interactions are small, but show definite trends, which correlate well with the p1/2 spinor population. Atomization energies are reduced significantly by incorporating spin-orbit interactions for all molecules considered. The KRHF calculations of several methylhalides demonstrate that the spinor energies from the KRHF method can be useful for the interpretation of experimental photoelectron spectra of molecules exhibiting spin-orbit splittings.   © 1998 John Wiley & Sons, Inc.   J Comput Chem 19: 1526-1533, 1998