Intramolecular Interaction Energies in Model Alanine and Glycine Tetrapeptides. Evaluation of Anisotropy, Polarization, and Correlation Effects. A Parallel Ab Initio HF/MP2, DFT, and Polarizable Molecular Mechanics Study NOHAD GRESH,1 SHERIF A. KAFAFI,2 JEAN-FRANC ¸ OIS TRUCHON,3 DENNIS R. SALAHUB3,4 1 Laboratoire de Pharmacochimie Mole´culaire et Cellulaire, FRE 2718 CNRS, Universite´ Rene´-Descartes, IFR Biome´dicale, 45, Rue des Saints-Pe`res, 75006 Paris, France 2 Density Functional Technologies Inc., 3803 Gateway T., Burtonsville, Maryland 20866 3 De´partement de Chimie, Universite´ de Montre´al, C.P. 6128, Montre´al, Que´bec H3X 2H9 Canada 4 University of Calgary, 2500 University Drive NW, Calgary T2N 1N4, Canada
Received 22 June 2003; Accepted 17 November 2003
Abstract: An extension of the SIBFA polarizable molecular mechanics procedure to flexible oligopeptides is reported. The procedure is evaluated by computing the relative conformational energies, ␦Econf, of the alanine tetrapeptide in 10 representative conformations, which were originally derived by Beachy et al. (J Am Chem Soc 1997, 119, 5908) to benchmark molecular mechanics procedures with respect to ab initio computations. In the present study, a particular emphasis is on the separable nature of the components of the energy and the particular impact of the polarization energy component on ␦Econf. We perform comparisons with respect to single-point HF, DFT, LMP2, and MP2 computations done at the SIBFA-derived energy minima. Such comparisons are made first for the 10 conformers derived from / torsional angle energy-minimization (the rigid rotor approach), and, in a second step, after allowing additional relaxation of the C␣ centered valence angles. In both series of energy-minimization, the SIBFA ␦Econf compared best with the LMP2 results using the 6-311G** basis set, the rms being 1.3 kcal/mol. In the absence of the polarization component, the rms is 3.5 kcal/mol. In both series of minimizations, the magnitudes of ␦Econf, computed as differences with respect to the most stable conformer taken as energy zero, decrease along the series: HF ⬎ DFT ⬎ LMP2 ⬎ SIBFA ⬎ MP2, indicative of increasing stabilization of the most highly folded conformers. © 2004 Wiley Periodicals, Inc.
J Comput Chem 25: 823– 834, 2004
Key words: model tetrapeptides; polarizable molecular mechanics study; ab initio HF/MP2
Introduction Two essential refinements to molecular mechanics are inclusion of anisotropy and of nonadditive effects, and the need for both is recognized. An essential determinant of anisotropy relates to the electrostatic contribution. Distributed multipoles (charges, dipoles, quadrupoles, and possibly octupoles) can be derived from the ab initio molecular wave function of a given molecule or molecular fragment and distributed on its atoms and chemical bonds.1 For peptides, this approach was shown to afford an accurate representation of electrostatics.2 In intermolecular interactions, the electrostatic contribution, Eel, is then computed as sums of multipole– multipole interaction energy terms. We have recently shown that,
in agreement with the results from ab initio Localized Møller– Plesset (LMP2) computations,3 such a representation of Eel was essential to account for the more favorable dimerization energy of N-methylacetamide than that of the -conformer of the alanyl dipeptide4 while conventional atom-centered point charges pre-
Correspondence to: N. Gresh; e-mail:
[email protected] Contract/grant sponsors: IDRIS, CINES, NSERC (Canada), and the Ligue Nationale Contre le Cancer (Comite´ Ile-de-France) This article includes Supplementary Material available from the authors upon request or via the Internet at http://www.interscience.wiley.com/ jpages/0192-8651/suppmat.
© 2004 Wiley Periodicals, Inc.
824
Gresh et al.
•
Vol. 25, No. 6
•
dicted the reverse preference. In fact, all the components of the ab initio binding energy are anisotropic. In the case of the SIBFA (Sum of Interactions Between Fragments Ab Initio computed) procedure,5 this can be accounted for, concerning the repulsion and charge-transfer short-range components, by an explicit introduction of heteroatom lone pairs, and concerning the polarization component, by the use of distributed polarizability tensors, along with the actual nonisotropic character of the polarizing field. The directional features of the HF components of the intermolecular interaction energy, ⌬E, were analyzed previously in the case of several Zn2⫹-ligand complexes using in parallel ab initio energy decomposition analyses.5e Nonadditivity is an essential feature of both multiply hydrogenbonded complexes and polycoordinated complexes of metal cations. Energy-decomposition analyses at the Hartree–Fock level with the Restricted Variational Space Analysis (RVS) procedure6 quantify the individual weights of polarization and charge-transfer to nonadditivity and evaluate the ability of polarizable molecular mechanics to account for them.7 In complexes involving charged entities, the correlation contribution was also found to be nonadditive.7b– c Very strong cooperative effects were computed for models of peptide H-bonded networks.8 As demonstrated by parallel ab initio HF/MP2, DFT, and polarizable molecular mechanics, these effects are strongly enhanced in the presence of a formate anion, the side chain of Asp or Glu residues, whose H-bonding interaction with the NH group of an N-methylformamide (NMF) can stabilize the formation of a linear array of NMF molecules modeling the ␣-helix.9 A molecular orbital study has shown cooperativity to be essential to ensure the preference for extended rather than cyclic complexation modes in crystals of 1,3-diones.10 The enhancement of intramolecular amide–amide hydrogen bonding through cooperativity was demonstrated by infrared and variable-temperature NMR spectroscopy.11 The importance of cooperativity in the formation of helices of types 310 and ␣, stabilized by hydrogen bonds between a carbonyl oxygen of a residue at position i and an amide hydrogen of the residue at position i⫹3 and i⫹4 respectively,12 and in the formation of sheets by -peptides13 was recently put forth in a DFT study of polyglycine models of up to 14 amino acid residues and underlined the need for semiempirical and molecular mechanics to account for such effects. Thus, inclusion of cooperativity appears as essential in view of reliable molecular mechanics/dynamics simulations of oligopeptide and protein folding. A promising step in this direction consists of the use of the fluctuating charge method within classical molecular mechanics/dynamics.14 Nevertheless, the need for models going beyond point-charges was mentioned.14d,e Ab initio and DFT methods are essential to validate molecular mechanics in conformational analysis, namely for the computation of actual intramolecular interactions. This was shown in the study by Beachy et al.,3 in which the energies of 10 conformers of the alanine tetrapeptide were computed by several different molecular mechanics force fields and compared to the ab initio LMP2 results. The present study evaluates the corresponding performances of the polarizable molecular mechanics procedure SIBFA in view of future extensions to larger oligopeptides and proteins. Our approach builds on the separable character of the SIBFA procedure, in which each of the components of the total energy has a representation of its own, grounded in ab initio quantum chemistry.
Journal of Computational Chemistry
Although energy decomposition analyses are an essential asset for the validation of polarizable molecular mechanics, such analyses can only be used for complexes between separate and distinct molecular entities. Thus they cannot be used to evaluate the variations of the relative weights of the individual components of the energy upon a conformational change or in competing conformations of a flexible molecule. On the other hand, within the SIBFA molecular mechanics procedure, a large molecule is assembled from a library of its rigid constitutive fragments, and the multipoles along the junctional bonds are relocated on the two bond extremities as well as on their midpoints, and the intramolecular energy is computed as the sum of interfragment interaction energies.5a,b The separable nature of the potential energy into its individual components can then lend itself to such an analysis, which was reported by us in the case of the Ala dipeptide.15 However, this analysis was done with an older representation of the intramolecular polarization, which is strongly damped at short distances. To fully account for cooperativity effects, and to ensure full consistency with the intermolecular interaction computations, it is necessary to use an identical formulation and calibration for the inter- and the intramolecular polarization energy contributions. In such a case, there are two mutually exclusive requirements. On the one hand, it is necessary to redistribute the multipoles and polarizabilities along the interfragment junctional bonds. On the other hand, it is necessary to retain nonfractional charges (0 for neutral entities, ⫾1 for charged ones) on the constitutive fragments, because otherwise this could result in divergences of the polarization component owing to the very short distances between two connecting fragments. We have recently reported studies on Zn2⫹ binding by flexible ligands having two or more Zn2⫹ coordinating groups, such as neutral and zwitterionic glycine,16a mercaptocarboxamides,16b,c and triphosphate.16d For such complexes, we resorted to an approach in which the Zn–ligand and ligand–ligand interactions are computed simultaneously as actual “intermolecular” interactions. To avoid overlaps between two connecting fragments, the H atoms and the bond barycenters of the junctional bonds were carried back on the C atoms from where they originate. This procedure was shown to afford a reproduction of the conformation-dependent HF and MP2 Zn2⫹ interaction energies of mercaptocarboxamides16b,c and triphosphate16d with a relative error of ⬍3%. Not accounting for intramolecular polarization effects simultaneously with the intermolecular ones resulted, on the other hand, in much higher relative errors, amounting to 10 –15%. As an extension of the previous approach to neutral molecules, the present study bears on model Ala and Gly tetrapeptides, built out of their N-methylformamide (NMF) and methane constitutive fragments. We now evaluate a procedure in which, similar to our previous treatment,5a,b the electrostatic contribution EMTP uses multipoles redistributed along the interfragment junctions. To compute the polarization component, on the other hand, another set of multipoles is used, in which, as in ref. 16, no redistribution is done along the junctions, and the hydrogens and bond barycenters are carried back on the covalently bound C or N atoms. In such a way, net nonfractional charges are preserved on the constitutive NMF and methane fragments for the computation of the Epol polarization contribution.
Intramolecular Interaction Energies in Model Alanine and Glycine Tetrapeptides
The present tests are done by computing the conformational energies of 10 conformers of the Ala tetrapeptide, which were originally derived by Beachy et al., and whose relative energetical ordering was used to assess the performances of empirical force fields3 as well as of semiempirical methods.17 These conformers are represented in Figure 1.
Procedure
825
On the other hand, for hydrogen-bonded systems such as the water dimer, the K2-BVWN method was found to parallel the Becke 88 three-parameter method.27 The fourth functional is the hybrid B3LYP exchange-correlation functional,34 that has recently been used by us, along with HF and MP2 computations, in tests of SIBFA on large molecular complexes.16c,d,35 The B3LYP computations used the 6-311G** and the cc-pVTZ(-f) basis sets. SIBFA Computations
Ab initio computations
The Hartree–Fock (HF) and MP2 computations used two basis sets: the CEP 4-31G(2d) developed by Stevens et al.18 and the 6-311G** basis set. These computations were done on the Gamess19 and Gaussian 9820 packages. We also performed HF and localized Møller–Plesset (LMP2) computations, using both 6-311G** and cc-pVTZ(-f) basis sets.21 These computations used the Jaguar code.22 The LMP2 computations are based on the approach developed by Pulay and Saebo.23 As indicated in the original articles, the occupied orbitals are localized, and the correlation basis is limited to the Atomic Orbitals (AO) in the spatial vicinity of the correlated pairs. Distant pairs are either neglected or treated at a lower level pair correlation. In addition to significant savings in CPU time, this approach reduces or eliminates the Basis Set Superposition Error (BSSE)24 and further enables to unravel the different contributions to the correlation energy in intermolecular interactions.25 The HF and LMP2 methods in Jaguar use a pseudospectral method to compute the integrals.26
DFT Computations
The DFT computations used four functionals: the first two ones are the Becke 8827 and the PLAP328 functionals. We used the TZVP basis set29 and the deMon code.30 The third functional, denoted by K2-BVWN, is the method recently developed by Kafafi et al.31 implemented in the Gaussian 96 suite of programs.32 All of the computations reported in this paper using this functional were done at the K2-BVWN/6-311G** level of theory. The advantages of this methodology over other DFT functionals have been reported earlier.31b– d Briefly, the K2-BVWN code completely eliminates the use of gradient corrections (GC) to the Vosko-Wilk-Nusair free electron gas correlation functional,33 while retaining the corresponding GC to exchange. In addition, along similar line as the popular Becke 88 method, the K2-BVWN approach empirically mixes 25% of the exact HF exchange to improve the asymptotic behavior of the functional at large distances.31b– d The above corrections have been shown to improve the performance of the K2-BVWN approach in handling weak, noncovalent, van der Waals systems such as the noble gas dimers, as well as charge transfer complexes. The latter systems are predicted by other DFT methods as noninteracting, or even strongly repulsive.
In the SIBFA procedure, a large, flexible, molecule is built from its constitutive fragments, and the intramolecular (conformational) energy is computed as a sum of intermolecular interfragment energies: ⌬E ⫽ E MTP ⫹ Erep ⫹ Epol ⫹ Ect ⫹ Edisp which denote, respectively, the first-order electrostatic (multipolar) and repulsion components, the polarization, charge-transfer, and dispersion terms. Upon a conformational change, the multipoles and polarizabilities undergo the appropriate matrix rotational changes. The multipoles (charges, dipoles, and quadrupoles) are derived from the ab initio SCF molecular orbitals on the fragment and are distributed on its atoms and bond barycenters using the procedure developed by Vigne´–Maeder and Claverie.1c The polarizabilities are distributed on the fragment’s localized molecular orbitals (bond barycenters and heteroatom lone pairs) using the procedure developed by Garmer and Stevens.36 The ab initio computations use the CEP 4-31G(2d) basis set and the GAMESS package.19 The constitutive fragments of an Ala oligopeptide are N-methylformamide (NMF) for the backbone and methane for the side chain. NMF is split into two fragments, formamide and methane, to allow for rotation along the N—C bond (the torsional angle), while retaining the multipolar distribution of NMF. We have resorted to such a representation for the backbone rather than resorting to two independent formamide and methane fragments in keeping with our results from our initial extensions of SIBFA to peptides5b to account for the inductive effects exerted on the formamide moiety due to N-methylation. In this process, two superimposing fictitious bonds are created, N—H and H—C. The multipoles on their H atoms are set to zero, while the multipoles on the two bond barycenters have each half the intensity of the multipoles of the original N—C bond. Similarly, each of the two junctional N—H and H—C bonds carries half the polarizability tensor of the original N—C bond. In this fashion, NMF as a unit retains a net charge of zero. The fact that the two separate NMF fragments do not retain a nonfractional charge is not detrimental for ␦Econf because the intra-NMF polarization energy is not counted, as specified below. To compute EMTP, the multipoles along the two C—H bonds connecting two successive NMF fragments are redistributed on the two C atoms of the newly created C—C bond as well as on its middle, following the procedure published in ref. 5a. The fragments used for the central dipeptide unit of the Ala tetrapeptide are shown in Figure 2 (the two backbone methyl groups at its N- and C-termini belong themselves
826
Gresh et al.
•
Vol. 25, No. 6
•
Chart 1. Value of the Valence Angles Used for the SIBFA Computations on Oligopeptides and of the Corresponding X-ray Values.
Journal of Computational Chemistry
Results and Discussion We have proceeded along the following steps:
C␣—C—N O—C—N C—N—C␣ C␣—C—O N—C␣—C N—C␣—C N—C␣—C
SIBFA
X-ray
114.0 125.0 123.0 121.0 109.5 109.5 109.5
115.8 124.2 120.1 119.7 109.2 107.9 107.9
to NMF fragments). Another set of multipoles is used to compute Epol, in which junctional bonds are shrunk by carrying back their H atoms and barycenters on their C or N atoms. This retains a net null charge on the NMF and methane fragments while simultaneously preventing divergences of Epol due to overlaps along the junctions. This procedure should be generalized to neutral as well as to charged peptide residues, because multipole redistribution to compute EMTP is confined to the junctional bonds, and there is no flow of charge beyond the junctions. The Ala tetrapeptide is therefore built out of 12 fragments: four NMF fragments each further split into a formamide-like and a methane-like fragment, one methane to cap the N-terminal residue, and three methane side chains. Upon computing Epol, the interactions between two successive fragments belonging to the same initial NMF entity are not taken into account. This choice was made to avoid double counting, because the mutual polarizations of the formamide-like and methane-like fragment should be reflected in the ab initio SCF multipolar distribution of the NMF entity. For consistency with our previous studies,5b,15 we have used the same bond lengths and valence angles as in these, as reported in chart 1. Overlay of the five heavy atoms (C, O, C␣, N, C␣) of a peptide unit built with this internal geometry on the corresponding unit built with X-ray values gives rise to a rms value of ⬍0.03 Å. Upon performing variations on the valence angles, we used the same values of 0 as in AMBER.37 The value of the valence angle bending force constant K0 was fit as follows. The Ala tetrapeptide was held in its extended conformation. Starting from 100°, 12 step-wise 2° variations were done on the N—C␣–C angle centered around the second C␣ carbon. The variations of the SIBFA intramolecular energies (without Edisp) were compared to those from corresponding 6-311G** HF calculations. We preferred to fit K0 based on uncorrelated rather than on correlated calculations due to some uncertainty regarding the most adequate correlated calculation for the fitting (see below). The value of K0 that best fit the HF intramolecular energy variations was 100 kcal/(mol 䡠 rad2). Such a value is about twice as large as in conventional force fields. As discussed below, this reflects the different formulation of the potential energy in SIBFA than in conventional molecular mechanics, and the inclusion of one to three interactions in the energy computations whenever atoms or atom centers 1 and 3 belong to different fragments. The Merlin package38 was used for the energy-minimizations.
1. Starting from the geometries made available from Ref. 3, we performed local SIBFA energy minimization of the C—N— C␣—C () and N—C␣—C—N () torsional angles, setting the bond lengths and valence angles at standard values. Singlepoint ab initio HF, MP2, as well as LMP2 and DFT computations were performed on the converged minima. Several distinct basis sets were used for the ab initio and DFT computations. This enabled us to evaluate in a consistent fashion, for each of the 10 geometries: (a) the validity of the introduction of Epol using the same formulation as in refs. 4, 5d–f, 7, 9, and 16, for intermolecular computations, and the improvement in the relative energetical ordering compared to the corresponding ordering in its absence; (b) the impact of using distributed multipoles rather than monopoles; (c) the weight of the contribution of correlation to the relative ordering, as obtained from the different quantum-chemical approaches and by molecular mechanics; (d) the influence of the functional used in the DFT computations; (e) the effect of the basis set. 2. Single-point SIBFA and ab initio computations were then done on the 10 conformers now removing the side-chain methyl and N- and C-terminal groups and replacing them by H atoms, thus passing from the Ala tetrapeptides to the “unprotected” Gly ones. The multipoles of the C-terminal formamide were those of an actual formamide fragment, and not those of the pseudoformamide fragment from NMF, as the latter fragment does not have a net null charge in the absence of its N-methyl substituent. Such removals of the methyl groups were done to evaluate the effects of the methyl groups on the energy ordering. These effects can translate into increases of: (a) short-range repulsion; (b) induction and/or polarizability; (c) correlation/dispersion. It was important to assess how well these effects translate in polarizable molecular mechanics as compared to ab initio results. 3. In present uses of the SIBFA procedure, conformational variations are done in the rigid rotor approximation. As a first step towards relaxing the valence angles in this procedure, we have allowed for variations of the C␣ valence angles along with the and torsional angles, and compared the resulting energy ordering following energy minimization to the ab initio ones on the so-derived converged minima. Energy minimizations on the and Torsional Angles
Figure 1 represents the 10 conformers investigated by Beachy et al., at the outcome of and energy minimization using the above-described procedure to compute EMTP and Epol. These energy-minimized conformations will be those used throughout to evaluate the conformational energies in the Ala and the Gly tetrapeptides prior to valence-angle minimization. Table 1 reports the conformational energy differences with respect to conformer I, in the absence of dispersion/correlation effects. It lists the ab initio HF values using the CEP 4-31G(2d), 6-311G**, and cc-pVTZ(-f) basis sets, and the SIBFA values without the Edisp term. The
Intramolecular Interaction Energies in Model Alanine and Glycine Tetrapeptides
Figure 1. Representation of the 10 energy minimized conformers of the alanine tetrapeptide.
SIBFA conformational energies are denoted as ␦E in the presence of intramolecular polarization Epol and as ␦E* in the absence of Epol. We also compute the corresponding value of ␦E, ␦Emon, when EMTP is limited to the monopolar contribution. Table 2 provides the values of the individual components of the SIBFA intramolecular energy. The essential conclusions from Table 1 are: (a) The CEP 4-31G(2d) and cc-pVTZ(-f) basis sets give the same ordering of
827
relative conformational energies: 1 ⬎ 2 ⬎ 4 ⬎ 5 ⫽ 6 ⬎ 3 ⬎ 7 ⬎ 8 ⬎ 10 ⬎ 9. The same ordering is found with the 6-311G** basis with the sole exception that 6 is found to have a slightly smaller (by 0.8 kcal/mol) conformational energy than 5. All three basis sets give very close values of the conformational energy differences ␦E computed with respect to conformer 1 taken as the energy zero. This could to some extent reflect the fact that the three basis sets have in common that their most diffuse s and p functions, as well as their polarization functions, are uncontracted Gaussians. It does not prejudge on the results that would obtain if different classes of polarization functions were used. (b) The ordering of conformers by ␦E(SIBFA) values is identical to the ab initio cc-pVTZ(-f) and CEP 4 31G(2d) ordering, and the values of ␦E are extremely close. The maximal error with respect to the cc-pVTZ(-f) basis set is 1.4 kcal/mol, with an rms of 0.7 kcal/mol. In the absence of the polarization component, the agreement is downgraded. Thus, the difference with respect to the cc-pVTZ(-f) basis set increases to up to 4.4 kcal/mol for folded conformer 3 and to 6 kcal/mol for conformers 9 and 10. Although conformer 5 has a lowered ␦E, conformer 6 has an increased ␦E, resulting in a 3.4 kcal/mol difference favoring 5, as contrasted to 0.4 in the presence of Epol and 0.1 with the cc-pVTZ(-f) basis set. This, along with the conclusions from Banks et al.14d underlines the need for an explicit Epol component, even concerning small uncharged molecules such as the alanine tetrapeptide. The present results also suggest that the introduction of Epol in the intramolecular interaction energy, with the same calibration as for intermolecular interactions, could be realized in a balanced fashion, without overestimating the stability of the most folded conformers. In addition, as a consequence of the separable nature of the SIBFA potential, this occurred without having to rescale the other energy components, which themselves retain the same calibration as in intermolecular interactions. It is also important to indicate that the impact of Epol on the relative E values is different from other polarizable molecular mechanics procedures. Thus, the
Table 1. Ala Tetrapeptide.
Ab initio HF
SIBFA
Conformer
4-31G(2d)
6-311G**
cc-pVTZ(-f)
⌬E
⌬E*
␦ E mono
1 2 3 4 5 6 7 8 9 10
0.0 1.0 10.2 3.2 7.3 7.5 13.4 17.6 30.0 28.2
0.0 1.1 9.3 3.1 7.3 6.1 12.2 16.2 28.4 26.6
0.0 1.3 10.5 3.4 7.6 7.7 13.7 17.8 29.8 28.5
0.0 0.7 10.8 3.2 7.6 7.8 11.9 18.0 29.8 28.9
0.0 1.0 14.9 3.2 6.5 9.8 13.4 21.8 35.7 34.6
0.0 ⫺0.7 6.1 1.9 6.3 9.2 12.4 22.0 27.4 38.1
Values of the HF and SIBFA (without the dispersion component) conformational energy differences ␦ E. Single-point ab initio computations are performed on the SIBFA minima. The ␦ E values (kcal/mol) are computed with respect to the energy of the most stable conformer taken as energy zero. To ease the evaluation, the most significant conformation energies are in bold type. ␦ E*: SIBFA energy value in the absence of E pol.
828
Gresh et al.
•
Vol. 25, No. 6
•
Journal of Computational Chemistry
Table 2. Values (kcal/mol) of the Intramolecular Interaction Energies in the 10 SIBFA and Energy-
Minimized Tetrapeptides and of Their Individual Components. Conformer
E MTP
E rep
E1
1 2 3 4 5 6 7 8 9 10
⫺3.3 ⫺2.3 ⫺8.7 ⫺0.3 3.0 0.1 ⫺1.3 ⫺1.8 ⫺17.4 1.9
114.7 114.6 135.0 114.9 114.9 121.1 126.0 135.0 164.5 144.1
111.4 112.3 126.3 114.6 117.9 121.2 124.7 133.2 147.1 146.0
E ct
E2
E1 ⫹ E2
⫺1.6 ⫺1.4 ⫺2.8 ⫺1.4 ⫺1.5 ⫺1.7 ⫺2.0 ⫺2.4 ⫺4.6 ⫺2.8
⫺8.2 ⫺8.4 ⫺12.3 ⫺8.2 ⫺7.1 ⫺10.1 ⫺9.6 ⫺12.0 ⫺14.1 ⫺13.9
103.2 103.9 114.0 106.4 110.8 111.1 115.1 121.2 133.0 132.1
E pol ⫺6.6 ⫺7.0 ⫺9.5 ⫺6.8 ⫺5.6 ⫺8.4 ⫺7.6 ⫺9.6 ⫺9.5 ⫺11.1
E disp
E tot
⫺92.5 ⫺92.7 ⫺98.3 ⫺93.1 ⫺93.5 ⫺96.0 ⫺97.9 ⫺98.5 ⫺102.2 ⫺101.6
10.7 11.2 15.7 13.3 17.3 15.1 17.2 22.7 30.8 30.5
A multiplicative factor of 0.8 is used for the E disp component.
energy difference E separating folded conformer 7 from extended conformer 1 is lowered by 1 kcal/mol upon including polarization. The corresponding difference actually increases by 2 and 1.2 kcal/mol in the OPLS/FQ14d and nonadditive AMBER 9939 forcefields respectively. (c) Limiting the electrostatic energy component to the monopolar terms, even in the presence of Epol, downgrades the agreement with respect to the ab initio conformational energies. Thus, the relative stabilities of 2 and 1 are reversed, conformer 3 has a significantly smaller ␦E value (6.1 kcal/mol compared to 10.5 from ab initio cc-pVTZ(-f) while conformer 10 has a ␦E that is increased by 10 kcal/mol. Not allowing for multipole redistribution, even in the presence of Epol, gives rise to even larger errors (unpublished). We have represented on Figure 3 the evolutions as a function of the conformation number, of the following relative conformational energies: ␦E(HF) with the CEP 4-31G(2d) basis
set; ␦E(SIBFA) without dispersion; ␦E(SIBFA*) without dispersion and without charge-transfer; and ␦(E1), the variations of the sole first-order SIBFA contribution (EMTP ⫹ Erep) with respect to the E1 value taken as energy zero. Although the least quantitative agreement with ␦E(HF) is that of ␦(E1), the explicit presence of Ect also appears necessary for the match of ␦E(SIBFA) with ␦E(HF). Table 3 lists the corresponding values when correlation effects are taken into account. These are: DFT results with the following functionals: Becke 88 /Perdew 86 (a), PLAP3 (b), K2-BVWN (c), and B3LYP with either the 6-311G** (d) or the cc pVTZ(-f) basis sets (e); LMP2 results with these respective basis sets (f) and (g) respectively; MP2 results with the 6-311G** basis set (h). These are compared with the SIBFA results (i) and (j) upon using a multiplicative factor of 1.0 (default option) or of 0.8 for the dispersion energy. The latter factor was previously found4 to enable the reproduction by SIBFA of the ⫺5.1 kcal/mol water–
Table 3. Ala Tetrapeptide.
DFT
LMP2
MP2
SIBFA
Conformer
B88/PD86 TZVP
PLAP3 TZVP
K2-BVWN 6-311G**
B3LYP/ 6-311G**
B3LYP/ cc-pVTZ(-f)
6-311G**
cc-pVTZ(-f)
6-311G**
a
b
1 2 3 4 5 6 7 8 9 10
0.0 0.7 6.9 3.5 7.7 7.2 11.5 15.3 22.3 24.0
0.0 1.0 11.4 3.7 8.1 9.5 14.7 19.9 33.0 32.5
0.0 0.7 7.7 2.8 7.0 6.8 11.2 15.0 23.9 24.0
0.0 0.7 6.6 2.7 7.2 5.4 9.8 13.6 22.5 21.4
0.0 1.0 7.7 3.0 7.3 6.6 11.1 15.0 23.7 23.6
0.0 0.6 5.8 1.9 5.8 3.9 7.5 11.7 21.5 17.6
0.0 3.4 10.8 2.1 5.8 4.2 9.7 16.3 28.3 20.2
0.0 ⫺0.1 1.5 1.1 4.4 0.8 3.1 6.6 15.4 10.7
0.0 0.3 3.4 2.3 6.5 3.6 5.2 10.7 17.5 17.5
0.0 0.5 5.0 2.6 6.6 4.4 6.5 12.0 20.1 19.8
Values of the DFT, LMP2, MP2 quantum-chemical, and SIBFA conformational energy differences ␦ E. Single-point ab initio computations are performed on the SIBFA minima. The ␦ E values (kcal/mol) are computed with respect to the energy of the most stable conformer taken as energy zero. To ease the evaluation, the most significant conformation energies are in bold type. a A multiplicative factor of 1 is used for the E disp component. b A multiplicative factor of 0.8 is used for the E disp component.
Intramolecular Interaction Energies in Model Alanine and Glycine Tetrapeptides
Figure 2. Representation of the fragments and their association to construct the central dipeptide unit. The dotted lines indicate the junctions along which multipole redistribution occurs.
829
water dimerization energy that resulted from a large basis set MP2 study of this dimer.40 The main results are as follows. With respect to the HF computations, the values of ␦E are reduced in magnitude, that is, the folded conformations have their relative stabilities improved relative to the unfolded conformations. However, with the present set of bond lengths and valence angles, the extended conformer 1 remains the most stable. The amount of ␦E reduction is the least pronounced with the DFT procedure, and the most pronounced with MP2. The DFT values with the B88/PW86, K2-BVWN, and B3LYP functionals give very close ␦E values. The ␦E/B3LYP values are not very sensitive to the basis set, whether cc-pVTZ(-f) or 6-311G**. Both the B88/PD86 and the ␦E/K2 values are closer to their B3LYP counterparts with the cc-pVTZ(-f) basis set than to
Figure 3. Evolution as a function of the conformation number of ␦E(HF), ␦E(SIBFA), ␦E*(SIBFA), and ␦(E1) (see text for definition). These curves are marked by rhombi, squares, triangles, and crosses, respectively.
830
Gresh et al.
•
Vol. 25, No. 6
•
the 6-311G* set. The PLAP3 functional, on the other hand, gives somewhat larger (app. 20%) ␦E values. With respect to their B3LYP counterparts, the LMP2 computations bring about a further reduction in ␦E values. However, these values are much more sensitive to the basis set used, cc-pVTZ(-f) or 6-311G**, than at the DFT level. This is particularly noteworthy for conformers 2 (3.4 kcal/mol vs. 0.6 for these respective basis sets), 3 (10.8 kcal/mol vs. 5.8), and 8 (16.3 kcal/mol vs. 11.7). This indicates a greater sensitivity of the cc-pVTZ(-f) basis set at the LMP2 level to the detailed features of the considered conformer, because in the geometries derived by Beachy et al. following full geometry optimization, both basis sets were in fact observed to give closely similar ␦E values. Such a sensitivity of the cc-pVTZ(-f) basis set is found even though the energy-minimized conformations at this stage (i.e., on the sole / dihedral angles) have small superposition rms with respect to their corresponding fully minimized counterparts whose geometries are available from Beachy et al.3 The average rms on heavy atoms is 0.3 Å over the 10 conformers, being smallest for the extended conformer 1 (0.1 Å) and largest for the highly folded conformer 10 (0.6 Å). An overlay of the SIBFA torsion angle energy-minimized conformations on the corresponding ones derived in ref. 1 is given as Supporting Information. The MP2 computations give rise to significantly larger reductions of ␦E than those from LMP2 computations. This is illustrated by the 6-311G** results for conformers 3 (1.5 vs. 5.8 kcal/mol), 6 (0.8 vs. 3.4 kcal/mol), 7 (3.1 vs. 6.9 kcal/mol), and 8 (6.6 vs. 11.3 kcal/mol). The magnitude of these differences is somewhat surprising. It is likely that the amount of intramolecular BSSE in MP2, a short-range effect, is larger the more folded the conformer. On the other hand, LMP2 while mostly devoid of BSSE effects,24 may not be capturing all the correlations between “distant” pairs, that is, the long-range stabilization energy from MP2, denoted as “dispersion.” To gain insight into the comparative MP2 and LMP2 interaction energy values, we have computed the formamide– formamide interactions in two extreme cases, namely a linear H-bonded dimer and a stacked complex at the interplanar separation of 3.3 Å. In the first, which privileges the C⫽O—HN interaction, the values of ⌬Eint with the 6-311G** basis set are ⫺6.2, ⫺6.2, and ⫺6.5 kcal/mol, at the DFT (B3LYP), LMP2 and MP2 levels, respectively. In the second, with no privileged interformamide interactions but maximized overlaps of the formamide planes, they amount to ⫺2.8, ⫺3.1, and ⫺4.1 kcal/mol, respectively. The LMP2 procedure in this extreme case would perform similar to the DFT one, which with the present functionals was found to underestimate the long-range stabilization energy in stacked complexes.41 In the case of tetrapeptides, the situation is less clear cut, because all interacting fragments being connected through chemical bonds, the distinction between “close” and “remote” entities is not clear. Due to the constraints imposed by the molecular structure, there can occur significant overlaps involving the four constitutive NMF entities. These give rise to interactions that are neither pure NH–CO hydrogen bonds nor pure stacking. Additional studies on model peptides and related molecules to further compare the results from correlated approaches in intramolecular interactions may appear necessary to clarify these points. The SIBFA results agree best with the 6-311G** LMP2 calculations with which they give a 1.3 kcal/mol rms.
Journal of Computational Chemistry
Table 4. Gly Tetrapeptide.
Conformer
4-31G(2d)
6-311G**
cc-pVTZ(-f)
SIBFA
1 2 3 4 5 6 7 8 9 10
0.0 2.3 3.3 3.9 7.7 8.0 8.8 10.5 7.0 19.9
0.0 2.1 2.2 3.8 7.4 7.5 7.5 9.1 5.6 18.2
0.0 2.3 3.1 4.0 7.8 8.1 8.8 10.4 6.5 19.8
0.0 2.5 3.1 3.4 7.4 7.0 8.1 11.1 7.7 21.7
Values of the HF and SIBFA (without the dispersion component) conformational energy differences ⌬E. Single-point ab initio computations are performed on the SIBFA minima. The ␦ E values (kcal/mol) are computed with respect to the energy of the most stable conformer taken as energy zero. To ease the evaluation, the most significant conformation energies are in bold type.
Overall, the order of increasing ␦E values is HF ⬎ DFT ⬎ LMP2 ⱖ SIBFA ⬎ MP2. Gly Tetrapeptides
We deemed it important to evaluate whether the presence of the methyl groups is correctly translated in terms of their effects on the conformational energy differences. Indeed, our initial calibration of the van der Waals radii of C atoms, which referred to the Ala dipeptide15 could have been biased by some possible improper representation of methyl/aliphatic C atoms. As mentioned in the Introduction, their effects could be threefold: repulsive, inductive, and van der Waals attractive. They should be translated at the levels of Erep, Epol/Ect, and Edisp respectively. We have redone single-point computations on the 10 tetrapeptide conformers, in which the three —CH3 side chains as well as the N- and Cterminal —CH3 groups were replaced by H atoms. The results are reported in Tables 4 and 5 along with HF DFT, LMP2, and MP2 calculations with different basis sets. The DFT calculations were done with the B3-LYP, K2-BVWN, and PLAP3 functionals. Tables 4 and 5 report the results obtained at the uncorrelated and correlated levels, respectively. At the uncorrelated level (Table 4), the rms between the ␦E(SIBFA) and those using the CEP 4-31G(2d), cc-pVTZ(-f), and 6-311G** basis sets amount to 0.8, 0.9, and 1.5 kcal/mol, respectively. As in the case of the Ala tetrapeptide, all three basis sets give closely similar ␦E values, and the ␦E values obtained with the CEP 4-31G(2d) basis set are closer to their cc-pVTZ(-f) than to their 6-311G** counterparts. The same trends are found as with the Ala tetrapeptide conformers concerning the relative ␦E values, whose magnitudes from the HF to the MP2 values are along the sequence HF ⬎ DFT ⬎ LMP2 ⫽ SIBFA ⬎ MP2. It is noted that owing to the absence of —CH3 side chains, conformer 3 is now more stable than conformer 1 at the LMP2 (6-311G** basis set), SIBFA, and MP2 levels. Again, as in the
Intramolecular Interaction Energies in Model Alanine and Glycine Tetrapeptides
831
Table 5. Gly Tetrapeptide.
DFT LMP2
MP2
Conformer
PLAP3 TZVP
K2-BVWN 6-311G**
B3LYP/ 6-311G**
6-311G**
cc-pVTZ(-f)
CEP 4-31G(2d)
6-311G**
SIBFA
1 2 3 4 5 6 7 8 9 10
0.0 1.8 3.6 3.7 7.8 9.1 9.1 10.9 7.4 21.0
0.0 1.8 1.2 3.4 6.9 7.1 6.9 8.1 3.2 16.4
0.0 1.7 0.5 3.3 7.2 6.7 6.1 7.3 1.9 14.3
0.0 1.2 ⴚ0.7 2.2 5.5 4.4 3.7 4.7 1.8 10.3
0.0 3.5 1.5 2.3 5.7 4.6 4.1 7.4 5.9 10.4
0.0 0.4 ⫺5.7 0.9 4.2 0.5 ⫺0.5 0.2 ⫺3.5 2.8
0.0 0.6 ⫺3.8 1.4 4.6 2.2 0.9 1.6 ⫺2.5 5.4
0.0 1.9 ⴚ1.2 2.5 6.5 4.2 4.4 6.8 1.8 14.8
Values of the DFT, LMP2, MP2 quantum-chemical, and SIBFA conformational energy differences ␦ E. Single-point ab initio computations are performed on the SIBFA minima. The ␦ E values (kcal/mol) are computed with respect to the energy of the most stable conformer taken as energy zero. To ease the evaluation, the most significant conformation energies are in bold type.
case of the Ala tetrapeptides, the best agreement of SIBFA is with the 6-311G** LMP2 calculations, with an rms of 1.6 kcal/mol. Concerning the LMP2 computations, we observe similar effects as with the Ala tetrapeptides, namely, the fact that some conformers have significantly larger ␦E values with the cc-pVTZ(-f) basis set than with the 6-311G**. This is the case for conformers 3 and 8 (1.5 vs. ⫺0.7 kcal/mol in the case of 3, and 7.4 vs. 4.7 kcal/mol in the case of 8). Such differences seem rather large considering the small size of the systems considered, even in the absence of geometry reoptimization. At the DFT level, and similar to the Ala tetrapeptides, the K2-BVWN functional gives a closer agreement with the B3LYP(6-311G**) results than PLAP3. Energy Minimization with Valence-Angle Relaxation
The SIBFA ␦E values were found to give a good agreement with the 6-311G** basis set results, at both HF and LMP2 levels. The results obtained at this stage, however, cannot be compared to those published by Beachy et al.,3 because energy minimization has so far involved only the and dihedral angles, while it was perfomed by Beachy et al. on the whole set of cartesian coordinates. The most obvious difference is in the fact that folded conformer 3 is less stable energetically than unfolded conformers 1, 2, and 4, while 3 is the global minimum from the Beachy et al. unrestrained energy minimizations. Furthermore, the compact conformers 9 and 10 are much higher energetically (␦E ⬎ 15 kcal/ mol) than from the unconstrained energy minimizations of Beachy et al., in which both 9 and 10 have ␦E values of 6.9 kcal/mol.3 The SIBFA procedure operating on the set of internal coordinates, inclusion of further relaxation in the procedure would be done through valence angle and bond length relaxation. This would require parameterization of the appropriate angle bending and stretching bond constants. Initial efforts along this direction (Creuzet and Langlet, unpublished) used the harmonic approxima-
tion, and we have adopted this approach here. We have limited this investigation to the variations, in the process of energy minimization, of the valence angles that are centered around the sp3 C␣ carbons along with the and dihedral angles. As noted above, the K0 parameters for C sp3 centered carbons were twice as large as in conventional force-field bending. This could reflect a greater sensitivity of SIBFA than conventional potentials to valence angle changes close to the interfragment junctions. There are two possible reasons for such a sensitivity. The first is the fact that in conventional atom-centered molecular mechanics, the interactions between atoms that are covalently linked to a common atom (denoted as one to three interactions) are not computed. Such interactions are computed in SIBFA, however, whenever the pair of corresponding atoms or centers belong to two distinct fragments. The second reason is the different formulation of the potential energy, which uses distributed multipoles and polarizabilities rather than atom-centered point charges, and the explicit introduction of short-range terms, Erep and Ect, which are represented by overlap-dependent formulas. At the uncorrelated level (Table 6), the SIBFA values are seen to remain in good agreement with both cc-pVTZ(-f) and 6-311G** basis set results, with rms values limited to 1.3 and 1.2 kcal/mol, respectively. This implies that the valence angle-bending constants parameterized for sp3 C are reasonable. The same observations can be made as in the absence of sp3 valence angle bending concerning the correlated computations. The SIBFA ␦E values compared best to their LMP2 counterparts with the 6-311G** basis set, with an rms of 1.2 kcal/mol, while the rms increases to 4.1 kcal/mol with the cc-pVTZ(-f) basis set. With respect to the DFT computations, the rms values amount to 3.3 and 2.3 kcal/mol, with respect to the cc-pVTZ(-f) and 6-311 G** computations. It is finally noted that, similar to the rigid rotor results from Table 1, the lowering of SE is more pronounced with the LMP2 than the DFT computations. Although ␦E for conformer 3 is smaller than with the rigid rotor
832
Gresh et al.
•
Vol. 25, No. 6
•
Journal of Computational Chemistry
Table 6. Ala Tetrapeptide.
Uncorrelated HF or SIBFA without E disp
DFT
LMP2
Conformer
6-311G**
cc-pVTZ(-f)
SIBFA
6-311G**
cc-pVTZ(-f)
6-311G**
cc-pVTZ(-f)
SIBFA
1 2 3 4 5 6 7 8 9 10
0.0 1.9 7.5 3.8 5.7 6.9 13.8 15.7 23.8 22.0
0.0 1.8 9.1 3.9 5.8 8.3 15.5 17.4 25.4 23.7
0.0 1.2 9.3 3.7 7.2 7.4 12.3 16.7 26.2 22.6
0.0 1.3 5.7 3.1 5.7 5.7 11.3 13.2 19.5 17.6
0.0 1.5 7.3 3.5 5.7 7.0 12.6 14.9 21.0 19.6
0.0 0.9 4.5 2.5 4.6 4.6 8.8 11.6 19.5 14.7
0.0 3.9 7.6 2.5 4.2 7.0 11.1 16.5 26.2 17.5
0.0 0.8 3.9 2.8 6.0 3.6 6.1 10.6 18.2 14.7
Comparisons between HF, LMP2, DFT, and SIBFA conformational energy differences ␦ E values following SIBFA energy-minimization of , , and sp3 C-centered valence angles. The ␦ E values (kcal/mol) are computed with respect to the energy of the most stable conformer taken as energy zero. To ease the evaluation, the most significant conformation energies are in bold type. a A multiplicative factor of 1 is used for the E disp component. b A multiplicative factor of 0.8 is used for the E disp component.
approach, it still remains less stable than conformers 1, 2, and 4, warranting a more complete relaxation. As a further step towards total relaxation of the tetrapeptides, energy minimization should be carried out on the three C—N— C␣, C—N—H, and H—N—C␣ N-centered valence angles, and on the three C␣—C⫽O, C␣—C—N, and O⫽C—N C carbonyl-centered valence angles. For the reasons commented above, the K0 bending constants from conventional molecular mechanics could not be used as such in the SIBFA potential. Instead, as reported in this study for the C␣-centered valence angles, their calibration should be undertaken on the individual valence angles of an N-methylformamide moiety embedded in a small peptide. Such a calibration is beyond the scope of this work, and will be reported elsewhere.
Conclusions and Perspectives The main objective of this work was to extend in consistent fashion the formulation of the SIBFA potential, designed for intermolecular interactions, to intramolecular (conformational) interactions within flexible molecules. It was essential to retain both anisotropy and nonadditivity features of the intermolecular potential. Regarding the Eel term, a key asset is the use of multicentered multipoles derived from quantum-chemical computations. The derivation of such multipoles for biomolecular modelling was pioneered by Dreyfus and Claverie42 and Rein et al.43 The evolution of Eel, as computed with cumulative atomic multipole moments, involving the interactions between rotated fragments upon conformational changes were found, in model heteroatom-containing molecules, to have good parallelism with respect to the corresponding ab initio HF curves upon inclusion of one to two interactions.44 (In the approach followed in ref. 44, the molecules were not built as done here using constitutive fragments). To attain more
quantitative agreement with these ab initio results, nevertheless, the need to account for induction effects was also mentioned.44 In the present study, we have elaborated on an approach in which, consistent with the original formulation of SIBFA5a,b the multipoles belonging to junctional H atoms or barycenters are redistributed along the interfragment junctions to compute the EMTP component, but carried back on the C or N atoms whence they originate to compute Epol. In this fashion, nonfractional charges (0 for neutral fragments, ⫾1 for charged ones) are retained for the interacting fragments. This should prevent divergences of Epol due to too short distances between two consecutive fragments as well as to nonfractional charges on the polarizing fragments. Such a procedure to compute Epol within a “pseudointermolecular” approach (treating all fragments within a flexible molecule as separate molecules connected together by shrunken C—H or N—H bonds) was justified by several tests on complexes of Zn2⫹ with flexible molecules having several Zn-binding ligands and their comparisons with ab initio values.16 Therefore, while permanent multipoles are used to compute the first-order electrostatic energy, the second-order polarization and charge-transfer contributions should enable one to account for the variations of multipole intensities that result from the construction of a large molecule from rigid fragments, as well as for their conformational dependencies.2b,c,f,g The issue of multipole transferability was raised for the first time in 1990.2b However, to our knowledge, aside from cation-flexible ligand complexes,16a– e analyses of its impact on actual interaction energies concerning neutral peptides were undertaken in only a few recent studies.2f,g Thus, in view of future simulations of oligopeptides, an important test of the present procedure related to the values of ␦E, the conformational energy differences, in 10 conformers of the alanine tetrapeptide, initially derived by Beachy et al.3 to benchmark molecular mechanics potentials with respect to ab initio computations. As a first step towards this evaluation, energy minimiza-
Intramolecular Interaction Energies in Model Alanine and Glycine Tetrapeptides
tion was done on the / angles and the conformational energies of the 10 conformers were compared to those from single-point computations from HF, DFT, LMP2, and MP2 computations. ␦E(SIBFA) in the absence of Edisp was found to predict the corresponding ␦E(HF) values with 4-31G(2d), 6-311G**, and cc-pVTZ(-f) basis sets with rms values ⬍1 kcal/mol. All three basis sets gave closely similar ␦E values, indicative of basis set independence at the HF level. In the absence of Epol, the agreement of ␦E(SIBFA) with HF values is significantly less, and the rms increased from 0.7 to 3.5 kcal/mol. Having verified that Epol was thus correctly introduced for intramolecular interactions, and that Ect is well behaved, we next evaluated the effect of Edisp on relative ␦Es. Presently, the effect of correlation is included in SIBFA through the Edisp component, because the multipoles and polarizabilities are those derived from the uncorrelated molecular orbitals of the constitutive fragments. The validity of this approximation for intermolecular interactions was shown previously in several computations devoted to polycoordinated complexes7a– c as well as multiply H-bonded complexes.7d,e,9 In contrast with the HF results, an unexpected result from single-point correlated quantum chemical computations was in the differing amplitudes of ␦E values, and the loss of basis set independence. Overall, the ␦E values were found to decrease in magnitude along the sequence HF ⬎ DFT ⬎ LMP2 ⫽ SIBFA ⬎ MP2. Such decreases reflect the increasing stabilization of folded conformers. It is likely that MP2 computations overestimate such a stabilization, because of intramolecular BSSE effects, but on the other hand, LMP2 and DFT may be providing a lower bound for the attractive effect of correlation. For the SIBFA computations, the best agreement was with the LMP2 computations with the 6-311G** basis set, with an rms of 1.3 kcal/mol. Single-point comparisons with the 10 corresponding conformers of the glycine tetramer, in which all —CH3 groups were removed, gave the same overall features as for the alanine tetramer. The closest agreement of ␦E(SIBFA) was again with the LMP2 calculations using the 6-311G** basis set. The retention by SIBFA of the same features as such calculations indicate that the effects of the methyl groups could be translated correctly. As a first step beyond the rigid rotor approach, additional relaxation of the angles centered around the sp3–C atoms was done. The same features obtained as in the rigid rotor approach. Thus, at the uncorrelated level, a close agreement obtains with both 6-311G** and cc-pVTZ(-f) basis sets, while at the correlated level the best agreement of ␦E(SIBFA) values (rms 1.2 kcal/mol) obtains with the LMP2 6-311G** basis set. Calibration of sp2 Cand N-centered valence angles will constitute the next step towards full relaxation of peptide conformations. The results of ␦E calculations using correlated multipoles and polarizabilities instead of uncorrelated ones will also be evaluated. Presently, we observed that the use of the Edisp component with a calibration factor of 0.8 allowed to translate correctly in SIBFA the evolution of ␦E (ab initio) upon passing from the HF level to the LMP2 one using the 6-311G** basis set. Additional tests of the approach will be done on larger oligopeptides, to assess if nonadditivity and cooperative/anticooperative effects could be correctly accounted for. These tests are now possible for 100 –200 atom complexes, which have become amenable to ab initio/DFT computations. Such effects come into play
833
in the preferential stabilization of an ␣-helix by a Glu or an Asp side-chain binding one main-chain NH hydrogen,9 and in the stabilization of double- or triple-stranded -sheets.12 They are also involved in the formation of complexes between divalent cations and metallooligopeptides such as Zn fingers or calcium-binding loops. The availability of ab initio-based polarizable molecular mechanics, in conjunction with that of reliable Continuum reaction field approaches for solvation effects,45,46 should constitute important steps towards predictions of the folding of short oligopeptides and possibly small proteins and nucleic acids. Three applications were recently published in which inter- and intramolecular interaction energies are computed simultaneously as a single integrated total energy.35a– c Supporting Information
S1 overlay of the SIBFA torsion angle energy-minimized conformations on the corresponding ones derived in ref. 1.
Acknowledgments The DFT and ab initio computations reported in this work using the Jaguar and Gaussian codes have been performed on the computers of I.D.R.I.S (Orsay, France) and C.I.N.E.S (Montpellier, France).
References 1. (a) Stone, A. J.; Alderton, M. Mol Phys 1985, 56, 1047; (b) Sokalski, W. A. ; Sawaryn, A. J Chem Phys 1987, 87, 526; (c) Vigne´–Maeder, F.; Claverie, P. J Chem Phys 1988, 88, 4934. 2. (a) Etchebest, C.; Lavery; R.; Pullman, A. Theoret Chim Acta 1982, 62, 17; (b) Faerman, C. H.; Price, S. L. J Am Chem Soc 1990, 112, 4915; (c) Price, S. L.; Stone, A. J. J Chem Soc Faraday Soc 1992, 88, 1755; (d) Dudek, M. J.; Ponder, J. W. J Comput Chem 1995, 16, 791; (e) Dardenne, L. E.; Werneck, A. S.; de Oliveira Neto, M.; Bisch, P. M. J Comput Chem 2001, 22, 689; (f) Kedzierski, P.; Sokalski, W. A. J Comput Chem 2001, 22, 1082; (g) Hermida–Ramon, J. Brdarski, S.; Karlstrom, G.; Berg, U. J Comput Chem 2003, 24, 161; (h) Whitehead, C. E.; Breneman, C. .M.; Sukumar, N.; Ryan, M. D. J Comput Chem 2003, 24, 512. 3. Beachy, M. D.; Chapman, D.; Murphy, R. B.; Halgren T. A.; Friesner, R. A. J Am Chem Soc 1997, 119, 5908. 4. Gresh, N.; Guo, H.; Salahub, D. R.; Roques, B. P.; Kafafi, S. A. J Am Chem Soc 1999, 121, 7885. 5. (a) Gresh, N.; Claverie, P.; Pullman, A. Theoret Chim Acta 1984, 66, 1; (b) Gresh, N.; Pullman, A.; Claverie, P.; Pullman, A. Theoret Chim Acta 1985, 67, 11; (c) Gresh, N.; Claverie, P.; Pullman, A. Int J Quantum Chem 1986, 29, 101; (d) Gresh, N.; Leboeuf, M.; Salahub, D. R. In Modelling the Hydrogen Bond; Smith, D. A., Ed.; American Chemical Society: Washington, DC, 1994, 82; (e) Gresh, N. J Comput Chem 1995, 16, 856; (f) Gresh, N.; Garmer, G. J Comput Chem 1996, 17, 1481. 6. Stevens, W. J.; Fink, W. Chem Phys Lett 1987, 139, 15. 7. (a) Garmer, D.; Gresh, N.; Roques, B. P. Proteins 1998, 31, 42; (b) Tiraboschi, G.; Roques, B. P.; Gresh, N. J Comput Chem 1999, 20, 1379; (c) Tiraboschi, G.; Gresh, N.; Giessner–Prettre, C.; Pedersen,
834
8.
9. 10. 11. 12. 13. 14.
15. 16.
17. 18. 19.
20.
21. 22. 23. 24. 25.
Gresh et al.
•
Vol. 25, No. 6
•
L. G.; Deerfield, D. W. J Comput Chem 2000, 21, 1011; (d) Gresh, N. J Phys Chem A 1997, 101, 8680; (e) Masella, M.; Gresh, N.; Flament, J. P. J Chem Soc Faraday Trans 1998, 94, 2745. (a) Guo, H.; Karplus, M. J Phys Chem 1992, 96, 7273; (b) Guo, H.; Karplus, M. J Phys Chem 1994, 98, 7104; (c) Suhai, S. J Chem Phys 1994, 104, 9766; (d) Edison, A. S.; Weinhold, F.; Markley, J. L. J Am Chem Soc 1995, 117, 12603; (e) Guo, H.; Sirois, S.; Proynov, E. I.; Salahub, D. R. In Theoretical Treatment of Hydrogen Bonding; Hadzi, D., Ed.; John Wiley & Sons Ltd: New York, 1997, Chapt 3; (f) Guo, H.; Salahub, D. R. Angew Chem Int Ed Engl 1998, 37, 2985; (g) Ludwig, R.; Reis, O.; Winter, R.; Weinhold, F.; Farrar, T. C. J Phys Chem B 1998, 102, 9312; (h) Kobko, N.; Paraskevas, L.; del Rio, E.; Dannenberg, J. J. J Am Chem Soc 2001, 123, 4348. Guo, H.; Gresh, N.; Roques, B. P.; Salahub, D. R. J Phys Chem B 2000, 105, 9746. Turi, L.; Dannenberg, J. J. J Phys Chem 1992, 96, 5819. (a) Gung, B. W.; Zhu, Z. Tetrahedron Lett 1996, 37, 2189; (b) Gung, B. W.; Zhu, Z.; Everingham, B. J Am Chem Soc 1997, 62, 3436. Wu, Y. D.; Zhao, Y. L. J Am Chem Soc 2001, 123, 5313. Lin, J.-Q.; Luo, S.-W., Wu, Y.-D. J Comput Chem 2002, 23, 1551. (a) Rappe´, A. Z.; Goddard, W. A. J Phys Chem 1991, 95, 3358; (b) Rick, S. W.; Stuart, S. J.; Berne, B. J. J Chem Phys 1994, 101, 6141; (c) Sternberg, U.; Koch, F.-T.; Molhoff, M. J Comput Chem 1994, 15, 524; (d) Banks, L.; Kaminski, G. A.; Zhou, R.; Mainz, D. T.; Berne, B. J.; Friesner, R. A. J Chem Phys 1999, 111, 741; (e) Kaminski, G. A.; Stern, H. A.; Berne, B. J.; Friesner, R. A.; Cao, Y. X.; Murphy, R. B.; Zhou, R.; Halgren, T. A. J Comput Chem 2002, 23, 1515. Gresh, N.; Tiraboschi, G.; Salahub, D. R Biopolymers 1998, 45, 405. (a) Rogalewicz, F.; Ohanessian, G.; Gresh, N. J Comput Chem 2000, 21, 963; (b) Tiraboschi, G.; Fournie´–Zaluski, M. C.; Roques, B. P.; Gresh, N. J Comput Chem 2001, 22, 1038; (c) Antony, J.; Gresh, N.; Hemmingsen, L.; Olsen, L.; Schofield, C.; Bauer, R. J Comput Chem 2002, 23, 1281; (d) Gresh, N.; Shi, G. B. J Comput Chem 2004, 25, 160. Mohle, K.; Hofmann, H.-J.; Thiel, W. J Comput Chem 2001, 5, 509. Stevens, W. J.; Basch, H.; Krauss, M. J Chem Phys 1984, 81, 6026. Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A., Jr. J Comput Chem 1993, 14, 1347. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andres, J. L.; Gonzalez, C.; Head–Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, Revision A.6; Gaussian, Inc.: Pittsburgh, PA, 1998. Dunning, T. H. J Chem Phys 1989, 90, 1007. Jaguar 4.1. Schrodinger Inc.: Portland, OR, 2000. (a) Saebo, S.; Pulay, P. J Chem Phys 1987, 86, 914; (b) Saebo, S.; Pulay, P. J Chem Phys 1988, 88, 1884. Saebo, S.; Tong, W.; Pulay, P. J Chem Phys 1993, 98, 2170. (a) Schutz, M.; Rauhut, G., Werner, H.-J. J Chem Phys A 1998, 102, 5997; (b) Reinhardt, P. Chem Phys Lett 2003, 370, 338.
Journal of Computational Chemistry
26. Murphy, R. B.; Beachy, M. D.; Friesner, R. A. J Chem Phys 1995, 103, 1481. 27. Becke, A. D. J Chem Phys 1988, 88, 1053. 28. (a) Proynov, E. I.; Sirois, S.; Salahub, D. R. Int J Quantum Chem 1997, 64, 427; (b) Perdew, J. P. ; Wang, Y. Phys Rev B 1986, 33, 8800; (c) Perdew, J. P. Phys Rev B 1986, 34, 7406E. 29. Godbout, N.; Salahub, D. R.; Andzelm, J.; Wimmer, E. Can J Chem 1992, 70, 560. 30. Andzelm, J. W.; Casida, M. E.; Koester, A.; Proynov, A., St-Amant, A., Salahub, D. R.; Duarte, H.; Godbout, N.; Guan, J.; Jamorski, C.; Leboeuf, M.; Malkin, V.; Malkina, O.; Sim, F. Vela, A demon software; University of Montre´al: Montreal, 1995. 31. (a) Kafafi, S. A.; El-Gharkawy, E. R. H. J Phys Chem A 1998, 102, 3202; (b) Kafafi, S. A. J Phys Chem A 102, 1998, 10404; (c) Kafafi, S. A.; Krauss, M. Int J Quantum Chem 1999, 75, 289; (d) Gregurick, S. K.; Kafafi, S. A. J Carbohydr Chem 1999, 18, 867. 32. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B.; Robb, M. A.; Cheeseman, J. R.; Keith, T.; Peterson, G. A; Montgomery, J. A., Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V ; Foresman, J. B.; Cioslowski, J.; Stefanov, B. B.; Nanayakkara, A.; Challacombe, M.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, J., Baker, J. P.; Stewart, J. P.; Head–Gordon, M.; Gonzalez, C.; Pople, J. A. Gaussian 96, Revision B.3; Gaussian Inc.: Pittsburgh, PA, 1998. 33. Vosko, S.; Wilk, L.; Nusair, M. Can J Phys 1980, 58, 1200. 34. (a) Lee, C.; Yang, W.; Parr, R. G. Phys Rev 1988, B37, 785; (b) Becke, A. J Chem Phys 1993, 98, 5648. 35. (a) Gresh, N.; Derreumaux, P. J Phys Chem B 2003, 108, 4862; (b) Gresh, N.; Sponer, J. E.; Spackova, N., Leszczynski, J.; Sponer, J. J Phys Chem B 2003, 107, 8669; (c) Ledecq, M.; Lebon, F.; Durant, F.; Marquez, A.; Giessner–Prettre, C.; Gresh, N. J Phys Chem B 2003, 107, 10640. 36. Garmer, D. R.; Stevens, W. J. J Phys Chem 1989, 93, 8263. 37. Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Fergusoin, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J Am Chem Soc 1995, 117, 5179. 38. Evangelakis, G. A.; Rizos, J. P.; Lagaris, I. E.; Demetropoulos, I. N. Comput Phys Commun 1987, 46, 401. 39. Wang, J.; Cieplak, P.; Kollman, P. A. J Comput Chem 2000, 21, 1049. 40. Feyereisen, M. W.; Feller, D.; Dixon, D. A. J Phys Chem 1996, 100, 2993. 41. Hobza, P.; Sponer, J. Chem Rev 1988, 99, 3247. 42. (a) Dreyfus, M. PhD Thesis, 1970, University of Paris; (b) Claverie, P. Thesis, Paris, 1973, CNRS library number A.O. 8214; (c) Claverie, P. In Localization and Delocalization in Quantum Chemistry; Chalvet, O.; Daudel, R.; Diner, S.; Malrieu, J. P., Eds.; Reidel: Dordrecht, 1976, p. 127, vol. II. 43. (a) Rein, R.; Rabinowitz, J. R.; Swissler, T. J. J Theoret Biol 1972, 34, 215; (b) Rein, R.; Swissler, T. J.; Renugopalakrishnan, V.; Pack, G. Jerusalem Symposia on Quantum Chemistry and Biochemistry, Jerusalem, 1973; (c) Rein, R. Adv Quantum Chem 1973, 7, 335. 44. Sokalski, W. A.; Lai, J.; Luo, N.; Shibata, M.; Ornstein, R.; Rein, R. Int J Quantum Chem Quantum Biol Symp 1991, 18, 61; (b) Sokalski, W. A.; Shibata, M.; Ornstein, R. L.; Rein, R. Theoret Chim Acta 1993, 85, 209; (c) Strasburger, K.; Sokalski, W. A. Chem Phys Lett 1994, 221, 129. 45. Langlet, J.; Claverie, P.; Caillet, J.; Pullman, A. J Phys Chem 1988, 92, 1631. 46. (a) Roux, B.; Simonson, T. Biophys Chem 1999, 78, 1; (b) Simonson, T. Curr Opin Struct Biol 2001, 11, 243.
