A Mechanistic Study of the Spontaneous Hydrolysis of Glycylserine as the Simplest Model for Protein Self‐Cleavage

A common feature of several classes of intrinsically reactive proteins with diverse biological functions is that they undergo self‐catalyzed reactions initiated by an N→O or N→S acyl shift of a peptide bond adjacent to a serine, threonine, or cysteine residue. In this study, we examine the N→O acyl shift initiated peptide‐bond hydrolysis at the serine residue on a model compound, glycylserine (GlySer), by means of DFT and ab initio methods. In the most favorable rate‐determining transition state, the serine ?COO^? group acts as a general base to accept a proton from the attacking ?OH function, which results in oxyoxazolidine ring closure. The calculated activation energy (29.4 kcal mol^?1) is in excellent agreement with the experimental value, 29.4 kcal mol^?1, determined by ¹H NMR measurements. A reaction mechanism for the entire process of GlySer dipeptide hydrolysis is also proposed. In the case of proteins, we found that when no other groups that may act as a general base are available, the N→O acyl shift mechanism might instead involve a water‐assisted proton transfer from the attacking serine ?OH group to the amide oxygen. However, the calculated energy barrier for this process is relatively high (33.6 kcal mol^?1), thus indicating that in absence of catalytic factors the peptide bond adjacent to serine is no longer a weak point in the protein backbone. An analogous rearrangement involving the amide N‐protonated form, rather than the principle zwitterion form of GlySer, was also considered as a model for the previously proposed mechanism of sea‐urchin sperm protein, enterokinase, and agrin (SEA) domain autoproteolysis. The calculated activation energy (14.3 kcal mol^?1) is significantly lower than the experimental value reported for SEA (≈21 kcal mol^?1), but is still in better agreement as compared to earlier theoretical attempts.     

Charge reversal of the conjugate base of formamide

Charge reversal collisional activation mass spectremetry of negative ions has been used in conjunction with positive ion collisional activation to investigate several isomeric [H₂, C, N, O]⁺ ions. Generation of m/z 44 ions from formamide, acetamide, JV-methylformamide, acetaldoxime and by charge reversal of the [M–1]^? ion formed from formamide yields several different isomeric structures. Charge reversal of the conjugate base of formamide appears to yield a mixture of singlet and triplet H₂NC?O⁺ ions; experiments with deuterium-labeled compounds have been used to support this. Ab initio molecular orbital calculations indicate that the triplet ion is a stable structure, existing in a potential minimum 390.6 kJ mol^?1 above the ground state singlet.     

Vapor-phase Beckmann rearrangement of oxime molecules over H-Faujasite zeolite.

The Beckmann rearrangement (BR) plays an important role in a variety of industries. The mechanism of this reaction rearrangement of oximes with different molecular sizes, specifically, the oximes of formaldehyde (H₂C?NOH), Z‐acetaldehyde (CH₃HC?NOH), E‐acetaldehyde (CH₃HC?NOH) and acetone (CH₃)₂C?NOH, catalyzed by the Faujasite zeolite is investigated by both the quantum cluster and embedded cluster approaches at the B3LYP level of theory using the 6‐31G (d,p) basis set. To enhance the energetic properties, single point calculations are undertaken at MP2/6‐311G(d,p). The rearrangement step, using the bare cluster model, is the rate determining step of the entire reaction of these oxime molecules of which the energy barrier is between 50–70 kcal mol^?1. The more accurate embedded cluster model, in which the effect of the zeolitic framework is included, yields as the rate determining step, the formaldehyde oxime reaction rearrangement with an energy barrier of 50.4 kcal mol^?1. With the inclusion of the methyl substitution at the carbon‐end of formaldehyde oxime, the rate determining step of the reaction becomes the 1,2 H‐shift step for Z‐acetaldehyde oxime (30.5 kcal mol^?1) and acetone oxime (31.2 kcal mol^?1), while, in the E‐acetaldehyde oxime, the rate determining step is either the 1,2 H‐shift (26.2 kcal mol^?1) or the rearrangement step (26.6 kcal mol^?1). These results signify the important role that the effect of the zeolite framework plays in lowering the activation energy by stabilizing all of the ionic species in the process. It should, however, be noted that the sizeable turnover of a reaction catalyzed by the Brønsted acid site might be delayed by the quantitatively high desorption energy of the product and readsorption of the reactant at the active center.     

A Highly Oxidizing and Isolable Oxoruthenium(V) Complex [RuV(N4O)(O)]2+: Electronic Structure,Redox Properties,and Oxidation Reactions Investigated by DFT Calculations

The electronic structure and redox properties of the highly oxidizing, isolable Ru^V?O complex [Ru^V(N₄O)(O)]²⁺, its oxidation reactions with saturated alkanes (cyclohexane and methane) and inorganic substrates (hydrochloric acid and water), and its intermolecular coupling reaction have been examined by DFT calculations. The oxidation reactions with cyclohexane and methane proceed through hydrogen atom transfer in a transition state with a calculated free energy barrier of 10.8 and 23.8 kcal mol^?1, respectively. The overall free energy activation barrier (ΔG^≠=25.5 kcal mol^?1) of oxidation of hydrochloric acid can be decomposed into two parts: the formation of [Ru^III(N₄O)(HOCl)]²⁺ (ΔG=15.0 kcal mol^?1) and the substitution of HOCl by a water molecule (ΔG^≠=10.5 kcal mol^?1). For water oxidation, nucleophilic attack on Ru^V?O by water, leading to O? O bond formation, has a free energy barrier of 24.0 kcal mol^?1, the major component of which comes from the cleavage of the H? OH bond of water. Intermolecular self‐coupling of two molecules of [Ru^V(N₄O)(O)]²⁺ leads to the [(N₄O)Ru^IV? O₂? Ru^III(N₄O)]⁴⁺ complex with a calculated free energy barrier of 12.0 kcal mol^?1.     

A theoretical study of C4H5N isomers: Comparison of cyano and isocyano as substituents

Ab initio molecular orbital calculations using a 3-21G basis set have been used to optimize geometries for pyrrole, CH₃(X)CCH₂, CH₃(H)CCHX (both cis and trans), c-C₃H₅X, and CH₂CHCH₂X, where X is CN and NC. In all the alkenyl derivatives methyl groups are found to adopt the conformation in which the methyl hydrogen eclipses the double bond. 6-31G*∥3-21G level calculations show the alkenyl cyanides to be of similar energy to pyrrole, but the isocyanides are ～20 kcal mol^?1 higher in energy. For both substituents the cyclopropyl derivatives are higher in energy by ～10 kcal mol^?1. At the 6-31G* level ring strain is 27.7 kcal mol^?1 for the cyanide and 30.6 kcal mol^?1 for the isocyanide. Data on the relative energies of RCN and RNC are compared when R is (i) a saturated hydrocarbon, (ii) an unsaturated hydrocarbon, (iii) an α-carbenium ion, (iv) an allyl cation, and (v) an α-carbanion.     

First Principles (DFT) Characterization of RhI/dppp‐Catalyzed CH Activation by Tandem 1,2‐Addition/1,4‐Rh Shift Reactions of Norbornene to Phenylboronic Acid

The C?H activation in the tandem, “merry‐go‐round”, [(dppp)Rh]‐catalyzed (dppp=1,3‐bis(diphenylphosphino)propane), four‐fold addition of norborene to PhB(OH)₂ has been postulated to occur by a C(alkyl)?H oxidative addition to square‐pyramidal Rh^III?H species, which in turn undergoes a C(aryl)?H reductive elimination. Our DFT calculations confirm the Rh^I/Rh^III mechanism. At the IEFPCM(toluene, 373.15 K)/PBE0/DGDZVP level of theory, the oxidative addition barrier was calculated to be 12.9 kcal mol^?1, and that of reductive elimination was 5.0 kcal mol^?1. The observed selectivity of the reaction correlates well with the relative energy barriers of the cycle steps. The higher barrier (20.9 kcal mol^?1) for norbornyl–Rh protonation ensures that the reaction is steered towards the 1,4‐shift (total barrier of 16.3 kcal mol^?1), acting as an equilibration shuttle. The carborhodation (13.2 kcal mol^?1) proceeds through a lower barrier than the protonation (16.7 kcal mol^?1) of the rearranged aryl–Rh species in the absence of o‐ or m‐substituents, ensuring multiple carborhodations take place. However, for 2,5‐dimethylphenyl, which was used as a model substrate, the barrier for carborhodation is increased to 19.4 kcal mol^?1, explaining the observed termination of the reaction at 1,2,3,4‐tetra(exo‐norborn‐2‐yl)benzene. Finally, calculations with (Z)‐2‐butene gave a carborhodation barrier of 20.2 kcal mol^?1, suggesting that carborhodation of non‐strained, open‐chain substrates would be disfavored relative to protonation.     

Evaluation of DFT Methods and Implicit Solvation Models for Anion‐Binding Host‐Guest Systems

Although supramolecular chemistry is traditionally an experimental discipline, computations have emerged as important tools for the understanding of supramolecules. We have explored how well commonly used density functional theory quantum mechanics and polarizable continuum solvation models can calculate binding affinities of host‐guest systems. We report the calculation of binding affinities for eight host–guest complexes and compare our results to experimentally measured binding free energies that span the range from ?2.3 to ?6.1 kcal mol^?1. These systems consist of four hosts (biotin[6]uril, triphenoxymethane, cryptand, and bis‐thiourea) with different halide ions (F^?, Cl^?, Br^?) in various media including organic and aqueous. The mean average deviation (MAD) of calculated from measured ΔG_a is 2.5 kcal mol^?1 when using B3LYP‐D3 with either CPCM or PCM. This MAD value lowers even more by eliminating two outliers: 1.1 kcal mol^?1 for CPCM and 1.2 kcal mol^?1 for PCM. The best DFT and implicit solvation model combination that we have studied is B3LYP?D3 with either CPCM or PCM.     

Addition of ethyl radicals to carbon monoxide. Kinetic and thermochemical properties of the propionyl radical

Azoethane was irradiated in the presence of carbon monoxide in the temperature range of 238 to 378 K. Kinetic parameters for the addition of ethyl radicals to carbon monoxide and for the decomposition of propionyl radicals were determined. The rate constants were found to be log k(cm³ mol^?1 sec^?1) = 11.19 - 4.8/θ and log k(sec^?1) = 12.77 - 14.4/θ, respectively. Estimated thermochemical properties of the propionyl radical are ΔH_f⁰ = -10.6 ± 1.0 kcal mol^?1, S⁰ = 77.3 ± 1.0 cal K^?1 mol^?1, and D(C₂H₅CO? H) = 87.4 kcal mol^?1.     

Ab initio calculations on phosphorus compounds. II. Effects of disubstitution on ligand apicophilicity in phosphoranes

Geometry optimizations at the HF/3-21G(*) and HF/6-31G* levels of ab initio theory have been carried out for various isomers of model disubstituted phosphoranes PH₃XY(X, Y?OH, CH₃, NH₂, and SH). Reasonable agreement was obtained between the optimized geometries and available crystal structure data for analogous compounds. The isomers were further characterized by frequency calculations. The MP2/6-31G*//6-31G* + ZPE energy data reveal that the interactions between the ligands are relatively small (0–4 kcal mol^?1) for the most stable conformations of the isomers. Hence, for these conformations the apicophilicities (based upon monosubstituted phosphoranes) are approximately additive. The less stable PH₃XY conformations are in general transition states or higher-order saddle points, and their interligand interactions are larger in magnitude (up to 10 kcal mol^?1); the results with these conformations suggest that apicophilicities may not be as additive for some highly substituted phosphoranes. © 1993 John Wiley & Sons, Inc.     

Mechanism of the Tyrosine Ammonia Lyase Reaction—Tandem Nucleophilic and Electrophilic Enhancement by a Proton Transfer

Quantum mechanics/molecular mechanics calculations in tyrosine ammonia lyase (TAL) ruled out the hypothetical Friedel–Crafts (FC) route for ammonia elimination from L ‐tyrosine due to the high energy of FC intermediates. The calculated pathway from the zwitterionic L ‐tyrosine‐binding state (0.0 kcal mol^?1) to the product‐binding state ((E)‐coumarate+H₂N? MIO; ?24.0 kcal mol^?1; MIO=3,5‐dihydro‐5‐methylidene‐4H‐imidazol‐4‐one) involves an intermediate (IS, ?19.9 kcal mol^?1), which has a covalent bond between the N atom of the substrate and MIO, as well as two transition states (TS1 and TS2). TS1 (14.4 kcal mol^?1) corresponds to a proton transfer from the substrate to the N1 atom of MIO by Tyr300? OH. Thus, a tandem nucleophilic activation of the substrate and electrophilic activation of MIO happens. TS2 (5.2 kcal mol^?1) indicates a concerted C? N bond breaking of the N‐MIO intermediate and deprotonation of the pro‐S β position by Tyr60. Calculations elucidate the role of enzymic bases (Tyr60 and Tyr300) and other catalytically relevant residues (Asn203, Arg303, and Asn333, Asn435), which are fully conserved in the amino acid sequences and in 3D structures of all known MIO‐containing ammonia lyases and 2,3‐aminomutases.     

Dynamic proton magnetic resonance studies on complex spin systems. Non-mutual three-spin exchange in N-(trideuteriomethyl)-2-cyanoaziridine

At room temperature and below, the proton NMR spectrum of N-(trideuteriomethyl)-2-cyanoaziridine consists of two superimposed ABC patterns assignable to two N-invertomers; a single time-averaged ABC pattern is observed at 158.9°C. The static parameters extracted from the spectra in the temperature range from –40.3 to 23.2°C and from the high-temperature spectrum permit the calculation of the thermodynamic quantities ΔH⁰ = ?475±20 cal mol^?1 (?1.987 ± 0.084 kJ mol^?1) and ΔS⁰ = 0.43±0.08 cal mol^?1 K^?1 (1.80±0.33 J mol^?1 K^?1) for the cis ? trans equilibrium. Bandshape analysis of the spectra broadened by non-mutual three-spin exchange in the temperature range from 39.4–137.8°C yields the activation parameters ΔH_tc^≠ = 17.52±0.18 kcal mol^?1 (73.30±0.75 kJ mol^?1), ΔS_tc^≠ = ?2.08±0.50 cal mol^?1 K^?1 (?8.70±2.09 J mol^?1 K^?1) and ΔG_tc^≠ (300 K) = 18.14±0.03 kcal mol^?1 (75.90±0.13 kJ mol^?1) for the trans → cis isomerization. An attempt is made to rationalize the observed entropy data in terms of the principles of statistical thermodynamics.     

Heat of formation of Benzophenone Oxide

The heat of formation of benzophenone oxide, Ph₂CO₂, was measured using photoacoustic calorimetry. The enthalpy of the reaction Ph₂CN₂ + O₂ → Ph₂CO₂ + N₂ was found to be ?48.0 ±0.8 kcal mol^?1 and ΔH_f(Ph₂CN₂) was determined by measuring the reaction enthalpy for Ph₂CN₂ + EtOH → Ph₂CHOEt + N₂ (?53.6 ±1.0 kcal mol^?1). Taking ΔH_f(PhCHOEt) = ?10.6 kcal mol^?1 led to ΔH_f(Ph₂CN₂) = 99.2 ± 1.5 kcal mol^?1 and hence to ΔH_f(Ph₂CO₂) = 51.1 ± 2.0 kcal mol^?1. The results imply that the self-reaction of benzophenone oxide i.e., 2Ph₂CO₂ → 2Ph₂CO + O₂ is exothermic by ?76.0 ±4.0 kcal mol^?1.     

Stabilization of HHeF by Complexation: Is it a Really Viable Strategy?

Ab initio calculations at the MP2 and CCSD(T) levels of theory have disclosed the conceivable existence of fluorine‐coordinated complexes of HHeF with alkali‐metal ions and molecules M⁺ (M⁺=Li⁺–Cs⁺), M⁺–OH₂, M⁺–NH₃ (M⁺=Li⁺, Na⁺), and MX (M=Li, Na; X=F, Cl, Br). All these ligands L induce a shortening of the H? He distance and a lengthening of the He? F distance accompanied by consistent blue‐ and redshifts, respectively, of the H? He and He? F stretching modes. These structural effects are qualitatively similar to those predicted for other investigated complexes of the noble gas hydrides HNgY, but are quantitatively more pronounced. For example, the blueshifts of the H? He stretching mode are exceptionally large, ranging between around 750 and 1000 cm^?1. The interactions of HHeF with the ligands investigated herein also enhance the (HHe)⁺F^? dipole character and produce large complexation energies of around 20–60 kcal mol^?1. Most of the HHeF–L complexes are indeed so stable that the three‐body dissociation of HHeF into H+He+F, exothermic by around 25–30 kcal mol^?1, becomes endothermic. This effect is, however, accompanied by a strong decrease in the H? He? F bending barrier. The complexation energies, ΔE, and the bending barriers, E*, are, in particular, related by the inverse relationship E*(kcal mol^?1)=6.9exp[?0.041ΔE(kcal mol^?1)]. Therefore the HHeF? L complexes, which are definitely stable with respect to H+He+F+L (ΔE≈25–30 kcal mol^?1), are predicted to have bending barriers of only 0.5–2 kcal mol^?1. Overall, our calculations cast doubt on the conceivable stabilization of HHeF by complexation.     

N.m.r. studies of a rate process in a bridged biphenyl: Resolution of a discrepancy between n.m.r. and polarimetric kinetic data

Total bandshape analysis of the temperature-dependent n.m.r. spectrum of the bridge methylene protons in 6,6-dicarbethoxy-1,2,3,4-dibenzcyclohepta-1,3-diene has been carried out to yield the following activation parameters: ΔG = 60.96 ± 0.17 kJmol^?1 (14.57 ± 0.04 kcal mol^?1), ΔH^? = 47.3 ± 0.8 kJ mol^?1 (11.3 ± 0.2 kcal mol^?1) and ΔS^? = ?45.6 ± 3.3 J mol^?1 K^?1 (?10.9 ± 0.8 cal mol^?1 K^?1). A value of ΔG^? in agreement with the above may be obtained from a bandshape analysis (at a suitable temperature) of the spectrum of the methylene protons in the carbethoxy groups. In our opinion, the rate process described by these activation parameters must be the configurational inversion of the biphenyl system, and thus the polarimetric data on this compound reported by Iffland and Siegel [J. Am. Chem. Soc. 80, 1947 (1958); half-life of 80 min at room temperature] are seriously called to question, as previously noted by Sutherland and Ramsay [Tetrahedron 21, 3401 (1965)]. We have attempted unsuccessfully to repeat the polarimetric work, and have been able to show by means of thin-layer chromatography in the cold, mass spectral analysis and polarimetry that an optically active impurity (as yet unidentified) is the most probable cause of the discrepancy.     

Anionic CH⋅⋅⋅X− Hydrogen Bonds: Origin of Their Strength,Geometry, and Other Properties

CF₃H as a proton donor was paired with a variety of anions, and its properties were assessed by MP2/aug‐cc‐pVDZ calculations. The binding energy of monoanions halide, NO₃^?, formate, acetate, HSO₄^?, and H₂PO₄^? lie in the 12–17 kcal mol^?1 range, although F^? is more strongly bound, by 26 kcal mol^?1. Dianions SO₄^2? and HPO₄^2? are bound by 27 kcal mol^?1, and trianion PO₄^3? by 45 kcal mol^?1. When two O atoms are available on the anion, the CH???O^? H‐bond (HB) is usually bifurcated, although asymmetrically. The CH bond is elongated and its stretching frequency redshifted in these ionic HBs, but the shift is reduced in the bifurcated structures. Slightly more than half of the binding energy is attributed to Coulombic attraction, with smaller contributions from induction and dispersion. The amount of charge transfer from the anions to the σ*(CH) orbital correlates with many of the other indicators of bond strength, such as binding energy, CH bond stretch, CH redshift, downfield NMR spectroscopic chemical shift of the bridging proton, and density at bond critical points.     

Data on the thermochemistry of azoalkanes

The rate constant of the primary decomposition step was determined for four symmetrical and four unsymmetrical azoalkanes. From the experimental activation energies and some literature enthalpy data, the following enthalpies of formation of radicals and group contributions were calculated: ΔH_? (CH₃N₂) = 51.5 ± 1.8 kcal mol^?1, ΔH_? (C₂H₅N₂) = 44.8 ± 2.5 kcal mol^?1, ΔH_? (2?C₃H₇N₂) = 37.9 ± 2.2 kcal mol^?1, [N_A-(C)] = 27.6 ± 3.7 kcal mol^?1, [N_A-(?_A) (C)] = 61.2 ± 3.1 kcal mol^?1.     

The fragmentations of [M-H]- anions derived from underivatised peptides. The side-chain loss of H2S from Cys. A joint experimental and theoretical study

Loss of H₂S is the characteristic Cys side‐chain fragmentation of the [M? H]^? anions of Cys‐containing peptides. A combination of experiment and theory suggests that this reaction is initiated from the Cys enolate anion as follows: RNH‐^?C(CH₂SH)CONHR′ Ø [RNHC(?CH₂)CONHR′ (HS^?)] Ø [RNHC(?CH₂)CO‐HNR′‐H]^?+H₂S. This process is facile. Calculations at the HF/6‐31G(d)//AM1 level of theory indicate that the initial anion needs only ≥20.1 kcal mol^?1 of excess energy to effect loss of H₂S. Loss of CH₂S is a minor process, RNHCH(CH₂SH)CON^?‐R′ Ø RNHCH(CH₂S^?)CONHR′ Ø RNH ^?CHCONHR+CH₂S, requiring an excess energy of ≥50.2 kcal mol^?1. When Cys occupies the C‐terminal end of a peptide, the major fragmentation from the [M–H]^? species involves loss of (H₂S+CO₂). A deuterium‐labelling study suggests that this could either be a charge‐remote reaction (a process which occurs remote from and uninfluenced by the charged centre in the molecule), or an anionic reaction initiated from the C‐terminal CO₂^? group. These processes have barriers requiring the starting material to have an excess energy of ≥79.6 (charge‐remote) or ≥67.1 (anion‐directed) kcal mol^?1, respectively, at the HF/6‐31G(d)//AM1 level of theory. The corresponding losses of CH₂O and H₂O from the [M? H]^? anions of Ser‐containing peptides require ≥35.6 and ≥44.4 kcal mol^?1 of excess energy (calculated at the AM1 level of theory), explaining why loss of CH₂O is the characteristic side‐chain loss of Ser in the negative ion mode. Copyright © 2003 John Wiley & Sons, Ltd.     

Formation of the iron-oxo hydroxylating species in the catalytic cycle of aromatic amino acid hydroxylases

The first part of the catalytic cycle of the pterin‐dependent, dioxygen‐using nonheme‐iron aromatic amino acid hydroxylases, leading to the Fe^IV?O hydroxylating intermediate, has been investigated by means of density functional theory. The starting structure in the present investigation is the water‐free Fe? O₂ complex cluster model that represents the catalytically competent form of the enzymes. A model for this structure was obtained in a previous study of water‐ligand dissociation from the hexacoordinate model complex of the X‐ray crystal structure of the catalytic domain of phenylalanine hydroxylase in complex with the cofactor (6R)‐L ‐erythro‐5,6,7,8‐tetrahydrobiopterin (BH₄) (PAH‐Fe^II‐BH₄). The O? O bond rupture and two‐electron oxidation of the cofactor are found to take place via a Fe‐O‐O‐BH₄ bridge structure that is formed in consecutive radical reactions involving a superoxide ion, O₂^?. The overall effective free‐energy barrier to formation of the Fe^IV?O species is calculated to be 13.9 kcal mol^?1, less than 2 kcal mol^?1 lower than that derived from experiment. The rate‐limiting step is associated with a one‐electron transfer from the cofactor to dioxygen, whereas the spin inversion needed to arrive at the quintet state in which the O? O bond cleavage is finalized, essentially proceeds without activation.     

Amide–π interactions between formamide and benzene

High‐level ab initio calculations have been carried out using a formamide–benzene model system to evaluate amide–π interactions. The interaction energies were estimated as a sum of the CCSD(T) correlation contribution and the HF energy at the complete basis set limit, for the geometries of the model structures at the energy minimum obtained by potential energy surface (PES) scans. NH/π geometry in a face‐on configuration was found to be the most attractive among the various geometries considered, with interaction energy of ?3.75 kcal/mol. An interaction energy of ?2.08 kcal/mol was calculated for the stacked N/Center type geometry, where the nitrogen atom of formamide points directly toward the center of the aromatic ring. The weakest C?O/π geometry, where a carbonyl oxygen atom points toward the plane of the aromatic ring, was found to have energy minimum at an intermolecular distance of 3.67 Å from the PES, with a repulsive interaction energy less than 1 kcal/mol. However, if there are simultaneous attractive interactions with other parts of the molecule besides the amide group, the weak repulsion could be easily overcome, to give a C?O/π geometry interaction. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009