- Email: [email protected]
Glycine and alanine: a theoretical study of solvent effects upon energetics and molecular response properties
Journal of Molecular Structure (Theochem) 500 (2000) 113–127 www.elsevier.nl/locate/theochem
Glycine and alanine: a theoretical study of solvent effects upon energetics and molecular response properties L. Gontrani, B. Mennucci, J. Tomasi* Dipartimento di Chimica e Chimica Industriale, Universita` di Pisa, Via del Risorgimento 35, 56126 Pisa, Italy
Abstract In the present work, we show that polarisable continuum model is able to reproduce the stability of zwitterionic forms of amino acids alanine and glycine in water solution at B3-LYP/6-31G (d) level of theory. The model is then extended to the calculation of vibrational frequencies, Vibrational circular dichroism spectra and nuclear magnetic resonance chemical shifts. The agreement with experimental data is good, except in the case of vibrations where specific hydrogen bond interactions are involved. In the latter case, a supermolecular approach may help in the predictions of some vibrational frequencies of the groups which form hydrogen bonds. 䉷 2000 Elsevier Science B.V. All rights reserved. Keywords: Polarisable continuum model; Solvation; Vibrational frequencies; Vibrational circular dichroism; Nuclear magnetic resonance chemical shifts
1. Introduction The polarisable continuum (PCM) solvation model [1] is widely used to accurately describe the solvent effect on the energy and on conformational barriers of various solutes, which do not undergo chemical transformation in solution. Recently, owing to the implementation of analytical first and second derivatives of the free energy [2–4], the calculation of various molecular properties has been made feasible, for instance vibrational/ Raman scattering frequencies and intensities (hyper) polarisabilities, nuclear shielding, ESR g-tensors and Vibrational circular dichroism (VCD) [5–7]. In this work, we present one of the first systematic applications of some of these new features to the compu* Corresponding author. Tel.: ⫹ 39-50-918244; fax: ⫹ 39-50502270. E-mail address: [email protected] (J. Tomasi).
tation of molecular properties in aqueous medium, selecting the two simplest molecules belonging to the class of amino acids, namely glycine and alanine, which show specific features in water solution, requiring a reappraisal of the solvation energy. The importance of amino acids in chemistry is well established [8]. As a matter of fact, they are among the most important biological molecules, since they constitute the “building blocks” of peptides and proteins. Most of their outstanding properties are due to their amphoteric nature, which derives from the presence, in their structure, of two functional groups of different polarity (the amino, NH2 and the carboxyl, COOH, groups). They are responsible, for example, of the existence of a tautomeric equilibrium between neutral and zwitterionic forms: in the gas phase (in case that the molecule vaporises without decomposing) the neutral form is energetically favoured, while in condensed phases, such as solutions in polar solvents and crystals, the ionic
0166-1280/00/$ - see front matter 䉷 2000 Elsevier Science B.V. All rights reserved. PII: S0166-128 0(00)00390-0
114
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
form is more stable [9,10]. These functional groups also play a fundamental role in the formation of hydrogen bonds, both with polar solvents and with other molecules in crystals; the presence of such interactions can be revealed through infrared (IR) or Raman spectroscopy and by nuclear magnetic resonance (NMR) experiments. Furthermore, all amino acids (except non-deuterated glycine) are chiral molecules and are naturally found in enantiomerically pure form (l). They exhibit, therefore, circular dichroism spectra. In particular, VCD has been extensively investigated recently, both experimentally and theoretically [11–16]. The choice of glycine and alanine, which, as stated before, are zwitterionic in water solutions, is rather challenging from the theoretical point of view because it is well known [17,18] that the definition of the cavity, which is at the basis of PCM model (see methodological section), turns out to be particularly crucial to get quantitative results on solvation energy in the case of anions and cations. Scope of this work is to verify and discuss how our model behaves in this particular case and to confront our findings with recent ab initio calculations that followed the more computationally expensive “supramolecular” approach [16].
2. Methodology and computational details The PCM [1,19] belongs to the category of solvation methods exploiting a continuum description of the solvent [20]. In brief, the physical system is represented by a charge distribution which describes the target molecule (or group of molecules) one has identified as solute. This charge occupies a given volume of space, called the molecular cavity, which is surrounded by an infinite continuum dielectric (the solvent) characterised by specific macroscopic properties (density, dielectric constant, refractive index, etc.). The interaction between the solute charge, or better the field it produces into the external volume occupied by the solvent, and the solvent itself, here described as a polarisable dielectric, is represented in terms of an apparent surface charge, s , spreading on the cavity surface. As a response, such charge produces a new field, called reaction field, which can modify the
previous state of the solute charge distribution (if the latter is not assumed to be composed of fixed point charges). The combined action of the solute field and the solvent reaction field can be analytically treated so as to obtain a final situation in which both solute and solvent distribution charges are mutually equilibrated. If this general picture is translated into a quantum mechanical language, the target solute is represented by its wavefunction and the interactions with the solvent are introduced in the Hamiltonian through specific perturbation, or reaction, operators. The Schro¨dinger equation to be solved becomes H 0 ⫹ V reac 兩C典 E兩C典
1
where H 0 is the Hamiltonian in the absence of the solvent. The solvent operator V reac acting on C is defined in terms of the surface apparent charge s , and thus, as said before, depends on the solute charge distribution (i.e. the solute wavefunction C ). This mutual interactions between C and V reac assures that the solution of Eq. (1) represents an equilibrated solute–solvent system. Here, we do not report any detail about the operators appearing in Eq. (1), and the procedures used to solve it, but we only recall that the solvent operator is defined as an integral over the cavity surface, to be solved analytically without any approximation. In the computational practice, the form of such integral, as well as of the whole apparent surface charge s , is simplified by discretising the cavity surface into small regions, called tesserae, of known area, ak. At each tessera k we can associate an apparent point charge, q k ak s , computed assuming that the surface charge s is constant in this small region. This discretisation allows to largely simplify all the operators involving the solvent. The main advantage of the PCM method is its almost immediate generalisability to different levels of the quantum mechanical description and to more complex physical systems involving, for example, external electromagnetic fields. In addition, the large freedom in the definition of the molecular cavity, that here has not to be limited to simple spherical or elliptical forms, but can be modelled according to the real shape of the solute geometry, allows a by far more realistic description of the solute–solvent system.
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
Thanks to all these peculiarities, in the last years the PCM model has been extended to the calculation of different molecular properties so as to make the solvated system a quantum mechanical problem, which can be studied with all the standard tools, originally defined for isolated molecules. This has first required the formulation of analytical gradients [2,3] and Hessians [4] taking into account the solvent effects, in order to allow the search, and the following characterisation, of all the interesting points on the potential energy surface of the system; geometry optimisations, definition of transition states, and vibrational analyses have then become available procedures exactly as in vacuo. Once this preliminary step is solved, the attention has been then shifted to more complex analyses, in particular, to the study of the combined action of the solvent and externally applied electric [21–24] and/or magnetic fields [25–28]. This research field is in fact characterised by numerous important applications, suffice it to cite the various spectroscopies, and thus it represents a big challenge for any theoretical solvation model. In the following analysis, besides to vibrational analysis, two specific external fields, and the related spectroscopies, will be considered; namely, an external magnetic field (to evaluate the nuclear shielding tensors determined in NMR measurements), and the combination of electric and magnetic fields (exploited in the VCD spectroscopy. To go a step further in the analysis of the molecular response properties of PCM solutes which are related to these processes, it becomes compulsory to define the free energy functional G. It is in fact this energetic quantity which has to be minimised with respect to the variational parameters defining the solute wavefunction so as to obtain the solution of Eq. (1); namely we can write G H0 ⫹
1 reac V 2
2
where the factor 1/2 takes into account the dependence of the solvent reaction operator on the solute wavefunction. For a PCM molecular solute G is the fundamental energetic quantity, as it univocally determines the behaviour of the molecular system in the presence of any internal and external perturbation. Thus, the
115
nuclear magnetic shielding which can be defined as the response of each nucleus X of the molecule to the contemporaneous effect of an external magnetic field and the nuclear magnetic moment, in the presence of the solvent becomes [24] X s ab
s 2G 2B a 2mXb
3
where Ba and mXb a; b x; y; z are the Cartesian components of the external magnetic field B, and of the nuclear magnetic moment of nucleus X. As in vacuo, also in condensed phases, the presence of a magnetic field introduces the problem of the definition of the origin of the related vector potential (historically called “gauge origin”). In ab initio calculations these problems can be dealt with introducing field dependent atomic orbitals (the so called GIAOs) [26–28]. The GIAO method has been recently applied to the PCM scheme [5], in order to obtain reliable nuclear shieldings for molecules in solutions. Still related to the GIAO approach, another recent extension of PCM is represented by the computation of VCD intensities for chiral molecules in solution [6]. In this case the external perturbation to consider is both magnetic and electric and the final effect is that the chiral sample responses to the left and right circularly polarised light in a different way, namely X Ri fi n ; n i 4 D1 n kn i
D1 (n ) is the differential absorption between the two differently polarised beams, n i and Ri are the excitation frequencies and the rotational strengths of the fundamental transition (0–1) in the ith normal mode, and the normalised band shape, fi, is assumed to be Lorentzian. Ab initio calculations of the frequencies and the rotational strengths are carried out within the harmonic approximation using the equation: Ri Imb具0兩mel 兩1典i ·具0兩mmag 兩1典c
5
where m el and m mag are the electric and the magnetic transition dipoles, respectively. These two quantities can be obtained by applying standard derivatives techniques (being the normal coordinates and the external magnetic field of the derivative parameters) that now, as we have shown
116
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
3. Results 3.1. Calculation of thermodynamic quantities for neutral and zwitterionic forms
Scheme 1.
above, are available also in the presence of a PCM solvent. In the present paper, we shall exploit a revised version of PCM called Integral Equation Formalism (IEF) [29–31] which represents a development of the original formulation [1] toward a more general and computationally efficient solvation model. In particular, in the following numerical study we shall report results obtained with the density functional extension of the PCM-IEF method, here implemented in a development version of Gaussian code [32]. All the features recalled above, geometry optimisation, vibrational analysis, nuclear shieldings and determination of VCD spectra, will be obtained in this framework by exploiting the Becke-3-LYP (B3LYP) hybrid functional with the 6-31G(d) basis set. As concerns the specificities of the solvation model, a cavity defined accordingly to the United Atom Topological Model (UATM) [33] will be exploited. In brief, this cavity is formed by interlocking spheres centred on selected nuclei (in our case C, N, and O) with radii defined according to the topological state of each nucleus (i.e. number and type of atoms bonded). The construction of the cavity and its tessellation are performed with the program gepol. [34] The solvent will be represented by a dielectric constant of 78.39 corresponding to the experimental value for liquid water at room temperature. To correctly assign normal modes of vibration, we performed a potential energy distribution (PED) calculation with the program gar2ped [35]. Table 1 Calculated free energy differences (DG) for selected amino acids processes (values in kcal/mol) Process
Glycine
Alanine
(1) NTaq ! ZWaq (2) NTg ! NTaq (3) NTaq ! ZWaq
⫺3.42 ⫺10.11 ⫺13.52
⫺4.48 ⫺7.42 ⫺11.90
Within the PCM framework, several thermochemical quantities can be calculated for amino acid systems. The processes we have selected may be summarised in the following Scheme 1: (the two symbols will identify the two tautomeric forms from now on). Core of any thermochemical calculation is the free energy of the system, since it is related (at constant pressure) to the equilibrium constant of a given chemical process and as it can be used to obtain other quantities through standard thermodynamics machinery. The free energy of a system in solution [20] may be decomposed into: electrostatic contribution Gel, cavitation term (related to the work necessary to create the cavity in the solvent) Gcav, dispersion free energy Gdis (due to dispersion interactions between solute and solvent), repulsion free energy Grep and a term due to molecular motions (thermal motions), GMm. The first term is provided by the PCM single point calculation (as explained in the methodological section); the Gcav, Gdis and Grep terms (often referred to as non-electrostatic contributions) are evaluated by exploiting semi-empirical relationships [36,37]; the last contribution (GMm) can be calculated applying the standard methods of statistical mechanics, once vibrational frequencies and moments of inertia for the system are known. The total free energy can therefore be expressed as follows: G Gel ⫹ Gcav ⫹ Gdis ⫹ Grep ⫺ RT ln qtr qrot qvib 6 In this formula, the subscripts identify the various components into which the molecular partition function can be factorised according to the Born–Oppenheimer approximation. Notice that the zero-point energy (ZPE) is included in the vibrational term. The computed results for all the processes indicated in Scheme 1 are collected in Table 1; the various contributions to G are listed in Table 2. The analysis of such data will be given in the following sections.
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127 Table 2 Calculated values of the various components of solvation free energy (kcal/mol)
Electrostatic (Gel)
Glycine NT
Glycine ZW
Alanine NT
Alanine ZW
⫺12.11
⫺14.15
⫺10.77
⫺13.94
(Values calculated with respect to neutral form optimised in the gas phase) Total non-electrostatic 1.41 ⫺0.14 2.13 0.63 Cavitation (Gcav) 11.83 10.45 14.31 12.99 Dispersion (Gdis) ⫺14.48 ⫺15.48 ⫺15.74 ⫺16.62 Repulsion (Grep) 4.06 4.89 3.56 4.25 Thermal motions at 32.02 32.19 48.38 48.58 298.15 K (GMm)
3.1.1. Calculation of D G for process 1 (NTaq ! ZWaq) As it was stated previously, amino acids are found to exist mainly in zwitterionic form in neutral water solutions. Unfortunately, the mole fraction of the neutral form in such media is so negligible that it cannot be revealed through experimental measurements. Therefore, the equilibrium constant of process 1 in solution is not known and thermodynamic quantities are not accessible. Nevertheless, they can be estimated. A possible approach was proposed by Haberfield [38], who suggested that, in most cases, thermodynamic quantities for a given compound (i.e. the free energy of formation, DfG) could be calculated from the values which the same property assumes in the “fragments” in which the molecule can be divided (“summation method”). Following this guideline, it is possible to devise a thermodynamic cycle which allows one to obtain DGZW-NT if gas phase acidities are known [39]. We report it briefly below; the reader is referred to the original article for further details. DGZW-NT, corresponding to step A, is equal to ⫺B ⫹ C ⫹ D ⫹ E ⫺ F (Scheme 2). Some of these quantities can be obtained by experimental measurements, others are estimated as
Scheme 2.
117
mentioned before. By following this procedure, Haberfield reported a value of ⫺7.67 kcal/mol for DGZW-NT of glycine. Another kind of estimate, which is performed through several equilibrium constant measurements, is based on the assumption that the esterification of the carboxylic group has limited effect on the basicity of the amino group [40]. In this case, DGZW-NT ⫺7:27 kcal=mol: Despite the results obtained by applying the “summation method” have been argued about [41] since it seems that the method induces significant errors (as large as 5–10 kcal/mol), we chose to refer our discussion to Haberfield’s findings because that procedure can be readily modified to obtain DGZW-NT values for every amino acid, provided that its gas phase acidity is known [39]. Moreover, in the case of glycine, the two approaches give very similar results and, as far as we know, the second method has not been applied to other amino acids yet. It should be clear, at this point, that estimation data in our possession are not very reliable as a test to check the validity of our computational results. We are aware that both “experimental” and theoretical results undergo large errors. Nonetheless, our aim is first of all to verify that zwitterions are stable molecules within the PCM computational scheme, so as to be able to use the geometries obtained theoretically for the evaluation of other properties, such as vibrational frequencies and chemical shifts. In the first step, the geometry of ZW and NT forms in water solution was optimised at the B3-LYP/631G(d) level of theory with the IEF version of PCM (see methodological section). The potential energy surface (PES) of neutral glycine and alanine has been studied in detail by several authors [42–46]. In spite of their apparent simplicity, these two molecules give rise in vacuo to quite complex potential surfaces which include several minima and transition states connecting them. In both cases, in the lowest energy structure, the two hydrogen atoms of the amino group lean towards a lone pair of carboxyl oxygen; the structure for glycine has Cs symmetry (see Fig. 1). The ZW case is rather more complex [47–49]. In fact, the zwitterionic form does not seem to exist in the gas phase as a stable structure; explicit water molecules [48], a crystal field or a continuum solvent distribution [36] is
118
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
Fig. 1. Optimised geometries of glycine and S-alanine B3-LYP/6-31G(d): left, neutral forms; right, Zwitterionic forms.
needed to stabilise it. Ding et al. [50] performed a detailed study of glycine ZW PES at the RHF/631G ⫹ (d) level of theory, in which they were able to characterise a local minimum with Cs symmetry, ˚ ) N–H···O hydrogen with a particularly short (1.55 A bond; this structure underwent proton transfer when the basis set was enlarged to Dunning’s DZP, thus moving to the lower energy neutral structure. Starting from that geometry, Jensen et al., in a more recent investigation, found that two water molecules are necessary to prevent the proton transfer at the MP2/ DZP⫹⫹//RHF/DZP level of theory [48]. As for alanine ZW, Sambrano et al. [49] recently dealt with intramolecular proton transfer in aqueous medium within the SCRF (self-consistent reaction field) approach by using an ellipsoidal cavity and multipolar expansion. The inclusion of both diffuse functions and electron correlation was deemed necessary to make the ZW a global minimum. Bearing in mind these apparently contradictory results, we chose in vacuo-type geometries (according
to Ding) as starting structures to perform our analysis, by applying the PCM method. For the neutral forms, we took the aforementioned lowest energy structures. The most significant optimised geometrical parameters for the two systems in water are reported in Table 3. Two features must be noted: the introduction of a methyl group in the side chain of the NT glycine induces a deviation from planarity of about 15⬚ in the O–C–C–N backbone; a similar distortion can be found passing to ZW too, owing to the formation of a ˚ long in glycine hydrogen bond. The latter is 1.834 A ˚ ZW and 1.800 A in alanine ZW. Such optimised geometries have been used for the vibrational analysis needed to calculate the thermal contributions to the free energies. As concerns this point some preliminary considerations have to be done. It is well known from statistical mechanics theory [51] that certain torsional low-frequency motions cannot be regarded as “vibrations” but should rather
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
119
Table 3 Relevant geometrical parameters for NT and ZW forms in water (see Fig. 1 for labels) a Parameter (NT) d N5–H6 d N5–H7 d C4–N5 d C1–C4 d C1–O2 d C1–O3 d O2–H8 a C1–C4–N5 d O3–C1–N4–N5 d H6–N5–C4–C1 d H7–N5–C4–C1 a
Glycine (Cs)
Alanine
1.019 1.019
1.020 1.021
1.452 1.525 1.355 1.211 0.976 112.62 0. 56.50 ⫺56.50
1.463 1.532 1.344 1.219 0.991 112.62 15.78 59.20 ⫺54.31
Parameter (ZW)
Glycine
Alanine
d N5–H6 d N5–H7 d N5–H8 d C4–N5 d C1–C4 d C1–O2 d C1–O3 d O3–H8 a C1–C4–N5 d O3–C1–N4–N5 d H6–N5–C4–C1 d H7–N5–C4–C1
1.030 1.029 1.053 1.503 1.553 1.247 1.265 1.834 102.59 8.53 ⫺12.67 103.31
1.030 1.030 1.054 1.515 1.561 1.247 1.266 1.800 105.60 ⫺9.13 12.52 ⫺103.51
˚ , angles in degrees. Bond lengths in A
be considered as hindered rotations. In polyatomic non-linear molecules, in fact, individual molecular groups may rotate relative to each other across the bond linking them. A barrier to internal rotation UIRMAX therefore exists. According to the magnitude of the barrier with respect to RT, the motion can be classified as free rotation UIR⫺MAX =RT 艑 1; hindered rotation UIR⫺MAX =RT 5–10 or ordinary normal mode of vibration UIR⫺MAX =RT ⬎ 10: The subject was studied by Pitzer and Gwinn [52], who compiled accurate tables containing the corrections to thermodynamic functions (the increase of the free energy from the free rotor value is reported) for the single hindered rotor. An extension of Pitzer procedure to multiple rotors, assuming negligible coupling between them, has recently been carried out by Ayala et al. [54] and it is now available in the latest version of the Gaussian package [32]. The procedure automatically identifies hindered rotation motions among the calculated normal modes and yields the necessary corrections to partition functions, by using simple expressions (of Lennard-Jones type) for the rotation potentials. Two hindered rotations were found for glycine, and three for alanine; in either case, they correspond to rotations around C–C and C–N bonds, whose harmonic frequencies occur between 70 and 320 cm ⫺1 (this result was predictable before running the calculation). Approximately the same result is obtained when the contribution for hindered rotation is evaluated through direct calculation of the various UIR-MAX by
performing a torsional energy scan at the same level of theory and then using Pitzer tables. The difference between the final computed DG is about 0.5 kcal/mol for alanine; we therefore decided to use the Ayala method. Still another contribution to free energy must be evaluated. In fact, even though the solvent is modelled as a continuum, it is made up of individual molecules, characterised by their internal and external motions. While moving, they can collide with the solute and exert a friction upon it; the solute–solvent interaction itself also plays a similar role. This kind of “medium hindrance” [53] can be described using Pitzer model as well. The increase in free energy was calculated to be about 2.3 kcal/mol by considering UIR⫺MAX =RT 5: By summing all the contributions to G, we found a value of ⫺3.42 kcal/mol for glycine DGZW-NT and ⫺4.48 kcal/mol for alanine (see Table 1). The DGZW-NT values reported in Table 1, although not in excellent agreement with “experimental” results (but see the former discussion upon the reliability of them), are able to correctly predict ZW stability without requiring explicit water molecules and to reproduce suitably the difference between alanine and glycine. In fact, the difference between the DGZW–NT of the molecules is equal to 0.68 kcal/mol from the estimates and 1.06 kcal/mol from our calculation. We intend to perform calculations on other amino acids to verify if the trend is respected and then calibrate our method against systematic errors that, according to us, may occur on both sides.
120
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
3.1.2. Calculation of D G for processes 2 (NTg ! ZWaq) and 3 (NTg ! NTaq) Gaffney et al. [55] performed an estimate of the enthalpy difference (DH) for the process (2) by calculating DHdiss ⫺ DHsub, the solution and sublimation heats of solid glycine, respectively, reporting a value of ⫺19.2 ^ 1 kcal/mol at 298 K (note that the uncertainty regards the sublimation heat only, while a relatively higher, but not quantifiable error, is to be associated with the solution heat). The calculation of this quantity within the PCM procedure was studied in the past by Bonaccorsi et al. [47,56], who appraised the purely electrostatic free energy (Gel) difference to be ⫺16.5 kcal/mol at the RHF/4-31G ⴱ level of theory; to obtain the required DH they evaluated the electrostatic entropic correction to G in solution with the relation 2 DG 8 DH DG ⫺ T 2T P Assuming that the temperature change of electrostatic DG is given only by modification of the dielectric constant 1 and of the cavity volume v, we obtain: 2 dv 2 d1 DG ⫺T DG DH DG ⫺ T 2v 21 P dT P dT 9 The values in brackets are calculated analytically, while the derivatives of v and 1 are taken from experimental data. They reported a value of 3.0 kcal/mol (see the original paper for details) for such correction, and ⫺19.5 kcal/mol for DH. The authors regarded the almost perfect numerical coincidence as due to error compensation, both from the theoretical point of view (neglect of thermal contributions to the free energy) and from the experimental one (the already cited large errors on solution enthalpy). We performed a free energy calculation for the neutral forms of glycine and alanine in gas phase, by optimising the geometry of the species and by subsequently calculating thermal contributions to G, as discussed previously for process (1). From the G values thus obtained for the gas phase species, we were able to appraise DG for processes (2) and (3). Similar calculations were also performed for alanine. A summary of all the quantities calculated is reported in Table 1.
A value of ⫺13.52 kcal/mol was calculated for DG of process 2 in the case of glycine. We also re-evaluated the electrostatic contribution (9), getting a value of ⫺4.07 kcal/mol; we computed, as well, the thermal motions entropy contribution from the calculated partition functions (⫺1.0 kcal/mol). Therefore, our estimate for DH of process (3) is ⫺18.57 kcal/mol. 3.2. Calculation of IR and VCD spectra The second subject we dealt with is the prediction of vibrational and VCD spectra of alanine and glycine in aqueous environment. As already cited in Section 1, the calculation of vibrational frequencies was performed at B3-LYP/6-31G (d) level on the geometries optimised with the same method. All the models were carried out following the harmonic approach and considering the solvent properties and the cavity unperturbed; moreover, we assume that the solvent is always equilibrated to the solute and that this equilibrium is preserved during vibrations. This hypothesis is reasonable for normal modes of intermediate frequency (bending region) while it may be unsatisfactory for high-frequency (stretching) and low-frequency (torsional) motions. In the former case, the solvent may not be able to attain equilibration during the solute motion; in the second case, for which large amplitudes are possible, there might be effects on the cavity size and the orientational motions of the solvent may perturb the solute displacements. In this study, we did not take into account any of these effects. Further details on the calculation of vibrational frequencies in solution can be found in a general review [20] and in a recent paper by our group [61]. 3.2.1. Vibrational spectra Among all the calculations we ran, we choose to report the aqueous ZW spectra only, because it is they that are worth discussing to assess the behaviour of the PCM method. The results are reported in Tables 4 and 5, together with the experimental data available and with a set of theoretical results obtained by Alper et al. [57,58] (glycine) and by Tajkhorshid et al. [16] for alanine. Experimental measurements go through several difficulties in water solution because of water absorptions interference (at least at the time the data we
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127 Table 4 Vibrational frequencies of Glycine ZW in water (cm ⫺1) PED HN–H2 ant. Str. HN–H2 sym. Str. C–H10 str. H2N–H str. ⫹ C–H9 str CH2 sym. str ⫹ H2N–H str. CO2⫺ ant. str. NH3 ant. bend. NH3 ant. bend. CH2 sciss. NH3 sym. bend. CH2 rock ⫹ NH3 sym. bend. CH2 wag. ⫹ NH7 bend. CH2 twist. ⫹ NH8 bend. CH2 wag. ⫹ CCN bend. HCC ⫹ HNC bend. C–N str. ⫹ CCN bend HN–H2 bend. ⫹ CHN bend. CH2 wag. ⫹ C–C str. ⫹ CO2 bend. CCN ⫹ CCH bend. Out of plane CO2 bend. ⫹ CCH bend. CCN ⫹ CCO bend. CCN bend. NH3 rotation NH3 and CH2 rotation a b c
Aa
Bb
Cc
3152 3058 2962 2968 1634 1615 1597 1533 1440 1410 1327 1315 1118 1119 1027 919 896
3397 3345 3126 3021 3011 1735 1659 1617 1454 1443 1388 1339 1289 1104 1086 983 906 868
3390 3366 3288 2924 2868 1831 1790 1700 1651 1547 1457 1437 1379 1257 1196 1032 1015 892
665 613
680 560
781 701
577 502 364 238
497 321 262 73
661 569 408 234
Takeda et al. [59], Ghazanfar et al. [60]. This work. Alper et al. [57,58]; scaled values.
possess [59,60] were collected); therefore, comparison with experiment is rather hampered. Whenever there are uncertainties in the datum, crystalline phase frequencies have been reported, as in other theoretical works [57,58]. The experimental spectrum of alanine was assigned by Diem et al. [11] by means of isotope substitution. To correctly assign the calculated frequencies, besides a graphical animation of vibrational modes with Gaussview package (included in Gaussian distribution [32]), a convenient way we encountered was to perform a PED calculation. We made it by employing the program gar2ped, a Gaussian output postprocessing utility written by Martin and Van Alsenoy [35] which uses a valence internal co-ordinate scheme proposed by Pulay and yields the fractional contribution of each defined co-ordinate within the calculated normal mode. Our assignment is reported in Tables 4
121
and 5 as well (only the major contributions to the mode are reported). As it can be seen, there is a fairly good agreement between calculated and experimental data (it is worth noticing that no scaling of the calculated values was applied). Yet, a couple of remarks are needed. We will limit our discussion to alanine spectrum, because only in this case an experimental assignment of normal modes is available (we consider it at all reliable, having been carried out by isotope substitution). First, the quality of the agreement is poorer for N– H stretching modes. In the experimental spectra, in fact, the three highest frequency modes (3080, 3060 and 3020 cm ⫺1) are attributed to the stretching vibrations of NH⫹ 3 group (two antisymmetrical and one symmetrical, respectively). The shift in frequency between these modes and the neutral gaseous form N– H stretching vibrations (occurring at about 3560 cm ⫺1[46]) is 500 cm ⫺1 on average. In our calculation, we obtained two N–H stretching frequencies in the upper part of the spectrum (3391 and 3338 cm ⫺1), which the PED calculation assigned to antisymmetrical and symmetrical stretching motions, respectively, of the two N–H bonds not involved in the intramolecular hydrogen bond (N5–H6 and N5–H7) and one N–H stretching frequency at 2992 cm ⫺1, which was assigned to a pure N5–H8 (that is, the H-bond acceptor fragment) stretching mode. On the contrary, the corresponding sixth normal mode (in decreasing order of frequency) of the experimental spectrum, which absorbs at 2962 cm ⫺1, was proved to be a methine stretching vibration. Therefore, it seems that the continuum alone is not able to account for such strong interactions with water (responsible of a shift of 500 wavenumbers, while only 150–200 are supplied by PCM). Moreover, in this case, the proximity of NH⫹ 3 and CO2 – groups is crucial. The resulting intramolecular H-bond (O3–H8) is much more effective than the interaction with the continuum solvent, even if the latter has a very high dielectric constant. This overestimate of a single N–H with respect to all the others is responsible for the permutation of the frequency order as well. We also performed a test calculation on a small cluster made up of an amino acid molecule and one water molecule in the vicinity of NH⫹ 3 and CO2 – groups; all the system was surrounded by a continuum
122
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
Table 5 A–C: Vibrational frequencies (cm ⫺1); D–E: VCD rotational strengths (esu 2 cm 2 × 10 44) for S-Alanine ZW in water PED a
Ab
Ba
Cc
Da
HN–H2 ant. str. HN–H2 sym. str. CH3 ant. str. CH3 ant. str. CH3 sym. str. H2N–H str. d Cm –H str. d CO2 ant. str. NH⫹ 3 bend. NH⫹ 3 bend. CH3 sym. bend. CH3 ant. bend. CH3 bend. ⫹ NH⫹ 3 bend. CH3 bend. ⫹ NH⫹ 3 bend. CH3 bend. ⫹ Cm –H bend. ⫹ CO2 s str. Cm –H bend. Cm –H bend. NH2 rock. ⫹ CH3 wag. C–N str. ⫹ CH3 bend. H2N–H bend. ⫹ Cm –H bend. NH⫹ 3 bend. ⫹ CH3 twist. HN–H2 bend. C–N str. ⫹ C–C str. CO2 bend. ⫹ OCC bend. Out of plane CO2 bend. C1–C–N bend. ⫹ C9–C–N bend. C–C–O bend. ⫹ C–C str. ⫹ C9–C–N bend. C–C–O bend. ⫹ C1–C–N bend. C–N tors. ⫹ C–CO2 tors. C–N tors. ⫹ C–CH3 tors. C–CH3 tors. CO2 –CH3 tors. ⫹ C–N tors.
3080 3060 3020 3003 2993 2962 2949 1645 1625 1607 1498 1459 1459 1410 1375
3391 3338 3149 3121 3051 2992 2952 1731 1664 1623 1520 1516 1445 1422 1388
3164 3142 3137 3135 3113 3093 3064 1759 1750 1695 1635 1536 1529 1441 1416
⫺5 19 0 0 ⫺3 34 ⫺20 ⫺36 42 ⫺100 ⫺11 ⫺2 ⫺19 5 104
⫺11 ⫺4 ⫺71 37 ⫺30 56 ⫺2 ⫺17 27 ⫺53 ⫺9 0 ⫺9 68 48
1351 1301 1220 1145 1110 1001 995 922 850 775 640
1364 1297 1213 1105 1084 1002 984 881 827 759 645
1390 1350 1290 1213 1133 1062 1028 918 840 775 635
⫺73 2 25 43 ⫺87 ⫺1 24 ⫺42 5 23 26
⫺171 38 50 1 ⫺11 3 ⫺2 ⫺22 20 22 ⫺8
527 477 399
520 381 342
625 528 432
⫺9 ⫺11 ⫺25
25 ⫺20 22
296 283 219 184
283 263 240 70
355 290 264 170
38 71 ⫺25 ⫺10
⫺33 50 ⫺30 ⫺9
a b c d
Ec
Present work. Diem et al. [11–13]. Tajkhorshid et al. [16]. See text.
distribution of water. The intramolecular · ··CO –H–bond was replaced by an interNH⫹ 2 3 molecular N–H· ··OH2 interaction, but the calculated frequency shift was only a little modified (2992– 2943 cm ⫺1). This result is compliant with recent findings by our group [61] about ROH···OR 0 interactions in water and methanol. In order to reproduce the observed frequency shifts, every H-bond site should
likely be allowed to interact directly with an explicit water molecule. In the rest of the spectrum, instead, the agreement of our frequencies is better than those calculated with the “super-molecular” approach [16]. Therefore, it seems that whenever the interaction is less specific, a continuum model is rather more efficient.
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
3.2.2. VCD spectra As it was already anticipated in the methodological section, the newly developed versions of PCM have recently made the calculation of several molecular response properties in solution feasible. Among them, VCD plays a fundamental role for the significant impact it has in the study of biomolecules. Amino acids, being the most easily available chiral compounds, have been the subject of extensive research for a couple of decades, both experimental and (more recently) theoretical [11–16]. Some amino acids VCD spectra show recurrent patterns (i.e. progression of positive and negative bands). One of the most important patterns can be found in the methine bending absorptions (between 1250 and 1450 cm ⫺1), which show a typical ( ⫹ , ⫺ , ⫹ ) pattern. It can be considered as a sort of “fingerprint” for the smallest amino acids (from deuterated glycine to proline) and for their polypeptides. Ab initio (or DFT) calculation of VCD spectra is a very delicate matter. In order to obtain good results, both for sign and for band intensities, optimisation and subsequent frequency calculation employing very large basis sets is recommended [62–64]. Such a design is hardly applicable in our case, since the stability of the zwitterionic form is largely dependent on the basis set and on the method used. Moreover, geometry convergence in these systems is already very difficult when a medium size basis is used, becoming almost prohibitive when the set is enlarged. Therefore, in the present work, we chose to calculate VCD spectra for l-alanine (corresponding to absolute S configuration) and S-Glycine-Cd1 at B3-LYP/6-31G (d) level of theory as well. Our findings are reported in Table 5 and Fig. 2a and b. In the former table, we confront our results with those of Tajkhorshid et al. [16], who employed a supermolecular model, while in Fig. 2a and b, the calculated VCD spectra of the two molecules in the methine bending region are drawn. As it can be seen (Table 5), the VCD spectra of lalanine yielded by the PCM calculation, is quite similar to that obtainable from a supermolecular model, with the exception of few vibrations (see also Section 3.2.1). Moreover, the spectra for both amino acids
123
(Fig. 2a and b) in the methine bending region reflects the experimental trend fairly well. 3.3. NMR chemical shift calculations Chemical shift calculations, in addition to being a theoretically interesting subject, may constitute nowadays a valid aid for complex NMR spectra interpretation. Recent works have shown that ab initio or density function calculations are able to give satisfactory results even for significantly large molecules, as taxol [65]. A great deal of experimental measurements is performed in solution; therefore, prevision of solvent effect is crucial. As it was mentioned in the methodological section, the newly implemented features of PCM allow one to calculate several properties, including nuclear shielding tensors. In this case, we decided to perform several calculations of shielding properties at the optimised B3-LYP/ 6-31G(d) geometry, employing different methods and basis sets. The reason for such a choice is that shielding properties are known to depend loosely on the geometry of the system and more tightly on the basis sets used; the quantum mechanical method applied also seems to play a role, but it has not been well established yet. We evaluated chemical shifts for B3-LYP/6-31G(d) optimised geometries by calculating the isotropic nuclear shielding of C and N atoms of glycine and alanine at the following levels of theory, namely HF/6-31G(d), HF/6-311⫹G(d,p), HF/6-311⫹G(2d,p), B3-LYP/6-31G(d), B3-LYP/6311⫹G(d,p), B3-LYP/6-311⫹G(2d,p) and then by subtracting these values from those calculated, at the same level of theory, for C atoms of neat tetramethylsilane and N atoms of neat nitromethane (data in our possession are referred to these standards). The values for neat liquids were obtained optimising the molecular geometry of TMS and CH3NO2 (at B3-LYP/6-31G(d) level), each embedded in the like liquid, described as a continuum according to the PCM scheme and then calculating the isotropic shielding with the required method. The calculated 13C and 15N chemical shift values for glycine and alanine in water solution are reported in Table 6, together with experimental values [66,67]. As it can be seen, in the case of 13C shifts, the best agreement was found using Hartree–Fock theory at the same basis set (smallest) used for geometry
124
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
Fig. 2. (a) Calculated VCD spectrum of S-Glycine-C-d1 in the methine bending region. (b) Calculated VCD spectrum of l-alanine in the methine bending region.
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
125
Table 6 13 C and 15N chemical shifts of glycine and alanine in water (in ppm with respect to TMS and nitromethane) Glycine
Expt. a
Calc.
Method
Alanine
Expt. a
Calc.
Method
C(yO)
173.2
171.1596 184.5958 184.1293 157.8906 178.1734 178.2840
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
C(yO)
176.5
173.6479 187.3207 186.7785 161.6340 182.1343 182.0875
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
C(–N)
42.2
42.8664 50.7837 51.4952 41.9151 47.1808 46.5477
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
C(–N)
51.3
50.4673 53.2342 52.5814 52.6163 59.0460 58.5760
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
350.2
468.6218 494.2359 488.4879 348.8527 380.1944 379.2182
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
Ca
17.0
16.9199 18.3043 17.7558 18.7769 20.6157 20.0401
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
458.2070 483.4378 477.7165 337.5204 367.7613 367.0281
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
N
N
a
338
Experimental data from [66,67].
optimisation at DFT level. The enlargement of the basis and the use of density functional did not lead to any improvement. When using DFT, the set augmentation seems to help only in the carbonyl case. The situation is reversed for 15N: Hartree–Fock methods yield poor values while B3-LYP behaves much better; even in this case, a larger basis set leads to worse results. The relative behaviour HF/ DFT with different nuclei has been already found to be unpredictable [65]. It is clear that our findings cannot easily be interpreted and we believe that the surprisingly good agreement of low-level calculations is due to a fortuitous cancellation of errors. Nevertheless, we can affirm that our calculated chemical shifts are, on average, in rather good agreement with the experiment, apart from 15N Hartree–Fock values.
4. Conclusions In this paper, we have described some of the capabilities of the PCM continuum solvation model, as well as some of its shortcomings. We have demonstrated that it is able to correctly reproduce solvent effects upon the energy of molecular systems, even in the “difficult” example of zwitterionic amino acids stability. A very good behaviour is found in the prediction of nuclear shielding properties as well, even though the results we obtained depend on the choice of the theoretical method and on the basis set. As regards vibrational properties and VCD, we have found that PCM yields a satisfactory description in all the cases in which the interactions between solute and solvent are not short-range and specific. In this case (i.e. the stretching modes of N–H fragments) a supermolecular model could give results in better agreement with experimental data.
126
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
Comprehensively, we may state that the model is a convenient and versatile method to account for solvent effects. Acknowledgements Partial support of CNR is here acknowledged. References [1] S. Miertus, E. Scrocco, J. Tomasi, J. Chem. Phys. 55 (1981) 117. [2] E. Cances, B. Mennucci, J. Chem. Phys. 109 (1998) 249. [3] E. Cances, B. Mennucci, J. Tomasi, J. Chem. Phys. 109 (1998) 260. [4] B. Mennucci, R. Cammi, J. Tomasi, J. Chem. Phys. 110 (1999) 6858. [5] R. Cammi, B. Mennucci, J. Tomasi, J. Chem. Phys. 110 (1999) 9877. [6] B. Mennucci, R. Cammi, J. Tomasi, P. Stephens, J.R. Cheeseman, M.J. Frisch, in preparation. [7] N. Rega, M. Cossi, V. Barone, J. Am. Chem. Soc. 119 (1997) 12 962. [8] P.J. Greenstein, M. Winitz, Chemistry of the Amino Acids, Wiley, New York, 1961, pp. 435–522. [9] K. Ijima, K. Tanaka, S. Onuma, J. Mol. Struct. 246 (1991) 257. [10] J. Almlo¨f, A. Kuick, O.J. Thomas, J. Chem. Phys. 59 (1973) 3901. [11] M. Diem, P.L. Polavarapu, M. Oboodi, L.A. Nafie, J. Am. Chem. Soc. 104 (1982) 3329. [12] B.B. Lal, M. Diem, P.L. Polavarapu, M. Oboodi, T.B. Freedman, L.A. Nafie, J. Am. Chem. Soc. 104 (1982) 3337. [13] T.B. Freedman, M. Diem, P.L. Polavarapu, L.A. Nafie, J. Am. Chem. Soc. 104 (1982) 3343. [14] M. Diem, J. Am. Chem. Soc. 110 (1988) 6967. [15] W.M. Zuk, T.B. Freedman, L.A. Nafie, J. Phys. Chem. 93 (1989) 1771. [16] E. Tajkhorshid, K.J. Jalkanen, S. Suhai, J. Phys. Chem. B 102 (1998) 5899. [17] J. Mestres, M. Solar, M. Orozco, J. Phys. Chem. 100 (1996) 606. [18] F.J. Luque, M. Orozco, Chem. Phys. 182 (1994) 237. [19] R. Cammi, J. Tomasi, J. Comp. Chem. 16 (1995) 1449. [20] J. Tomasi, M. Persico, Chem. Rev. 94 (1994) 2027. [21] R. Cammi, M. Cossi, J. Tomasi, J. Chem. Phys. 104 (1996) 4611. [22] R. Cammi, M. Cossi, B. Mennucci, J. Tomasi, J. Chem. Phys. 105 (1996) 10 556. [23] R. Cammi, B. Mennucci, J. Tomasi, J. Am. Chem. Soc. 34 (1998) 8834. [24] R. Cammi, B. Mennucci, J. Tomasi, J. Phys. Chem. A102 (1998) 870. [25] R. Cammi, J. Chem. Phys. 109 (1998) 3185.
[26] H.F. Hameka, Rev. Mod. Phys. 34 (1962) 87. [27] R. Ditchfield, Mol. Phys. 27 (1974) 789. [28] K. Wolinski, J.F. Hinton, P. Pulay, J. Am. Chem. Soc. 112 (1990) 8251. [29] E. Cances, B. Mennucci, J. Math. Chem. 23 (1998) 309. [30] E. Cances, B. Mennucci, J. Tomasi, J. Chem. Phys. 107 (1997) 3031. [31] B. Mennucci, E. Cances, J. Tomasi, J. Phys. Chem. B 101 (1997) 10 506. [32] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery, R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam, A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, C.S. Clifford, J. Ochterski, G.A. Petersson, P.Y. Ayala, Q. Cui, K. Morokuma, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J.V. Ortiz, A.G. Baboul, J. Cioslowski, B.B. Stefanov, G. Liu, C.A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, C.M.W. Wong, J.L. Andres, M. Head-Gordon, E.S. Replogle, J.A. Pople, Gaussian, Inc., Pittsburgh, PA, 1998; gaussian 99, Development Version (Revision A.10 ⫹ ). [33] V. Barone, M. Cossi, J. Tomasi, J. Chem. Phys. 107 (1997) 3210. [34] J.L. Pascual-Ahuir, E. Silla, J. Tomasi, R. Bonaccorsi, J. Comp. Chem. 8 (1987) 778. [35] J.M.L. Martin, C. Van Alsenoy, GAR2PED, University of Antwerp, 1995. [36] F.M. Floris, J. Tomasi, J. Comput. Chem. 10 (1989) 616. [37] F.M. Floris, J. Tomasi, J.L. Pascual-Ahuir, J. Comput. Chem. 12 (1991) 784. [38] P. Haberfield, J. Chem. Ed. 57 (1980) 347. [39] R.A.J. O’Hair, J.H. Bowie, S. Gronert, Int. J. Mass Spectrom. Ion Processes 117 (1992) 23. [40] G. Wada, E. Tamura, M. Okina, M. Makamura, Bull. Chem. Soc Jpn. 55 (1982) 3064. [41] M.J. Locke, R.T. McIver Jr., J. Am. Chem. Soc. 105 (1983) 4226. [42] P. Palla, C. Petrongolo, J. Tomasi, J. Phys. Chem. 84 (1980) 435. [43] V. Barone, F. Lelj, C. Adamo, J. Chem. Phys. 102 (1995) 364. [44] D.T. Nguyen, A.C. Scheiner, J.W. Andzelm, S. Sirois, D.R. Salahub, A.T. Hagler, J. Comput. Chem. 18 (1997) 1609. [45] S.G. Stepanian, I.D. Reva, E.D. Radchenko, M.T.S. Rosado, M.L.T.S. Duarte, R. Fausto, L. Adamowicz, J. Phys. Chem. A102 (1998) 1041. [46] S.G. Stepanian, I.D. Reva, E.D. Radchenko, M.T.S. Rosado, J. Phys. Chem. A102 (1998) 4623. [47] R. Bonaccorsi, P. Palla, J. Tomasi, J. Am. Chem. Soc. 106 (1984) 1945. [48] J.H. Jensen, M.S. Gordon, J. Am. Chem. Soc. 117 (1995) 8159. [49] J.R. Sambrano, A.R. de Sousa, J.J. Queralt, J. Andre´s, E. Longo, Chem. Phys. Lett. 294 (1998) 1.
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127 [50] Y. Ding, K. Krogh, Jespersen, Chem. Phys. Lett. 199 (1992) 261. [51] G.N. Lewis, M. Randall, Thermodynamics, 2nd ed., McGrawHill, New York, 1961 (revised by K.S. Pitzer, L. Brewer). [52] K.S. Pitzer, D.W. Gwinn, J. Chem. Phys. 10 (1942) 428. [53] O.N. Ventura, A. Lledos, R. Bonaccorsi, J. Bertran, J. Tomasi, Theor. Chim. Acta 72 (1987) 175. [54] P.Y. Ayala, H.B. Schlegel, J. Chem. Phys. 108 (1998) 2318. [55] J.S. Gaffney, R. Pierce, L. Friedman, J. Am. Chem. Soc. 99 (1977) 4293. [56] R. Bonaccorsi, E. Scrocco, J. Tomasi, J. Biosci. 8 (1985) 627 (Proceedings of the International Symposium on Biomolecular Structural Interactions, Suppl.) [57] J.S. Alper, H. Dothe, D.F. Coker, Chem. Phys. 153 (1991) 51. [58] J.S. Alper, H. Dothe, D.F. Coker, Chem. Phys. 161 (1992) 199. [59] M. Takeda, R.E.S. Iavazzo, D. Garfinkel, I.H. Scheinberg, J.T. Edsall, J. Am. Chem. Soc. 80 (1958) 3813.
127
[60] S.A.S. Ghazanfar, D.V. Myers, J.T. Edsall, J. Am. Chem. Soc. 86 (1964) 3439. [61] C. Cappelli, B. Mennucci, C.O. Da Silva, J. Tomasi, J. Chem. Phys. 12 (2000) (in press). [62] J.R. Cheeseman, M.J. Frisch, F.J. Devlin, P.J. Stephens, Chem. Phys. Lett. 252 (1996) 211. [63] F.J. Devlin, P.J. Stephens, J.R. Cheeseman, M.J. Frisch, J. Phys. Chem. A101 (1997) 6322. [64] F.J. Devlin, P.J. Stephens, J.R. Cheeseman, M.J. Frisch, J. Phys. Chem. A101 (1997) 9912. [65] J.R. Cheeseman, G.W. Trucks, T.A. Keith, M.J. Frisch, J. Chem. Phys. 104 (1996) 5497. [66] J.B. Stothers, Carbon-13 NMR Spectroscopy, Academic Press, London, 1986. [67] J.B. Stothers, Annual Reports in NMR Spectroscopy, vol. 18, Academic Press, London, 1986, p. 365.
Glycine and alanine: a theoretical study of solvent effects upon energetics and molecular response properties L. Gontrani, B. Mennucci, J. Tomasi* Dipartimento di Chimica e Chimica Industriale, Universita` di Pisa, Via del Risorgimento 35, 56126 Pisa, Italy
Abstract In the present work, we show that polarisable continuum model is able to reproduce the stability of zwitterionic forms of amino acids alanine and glycine in water solution at B3-LYP/6-31G (d) level of theory. The model is then extended to the calculation of vibrational frequencies, Vibrational circular dichroism spectra and nuclear magnetic resonance chemical shifts. The agreement with experimental data is good, except in the case of vibrations where specific hydrogen bond interactions are involved. In the latter case, a supermolecular approach may help in the predictions of some vibrational frequencies of the groups which form hydrogen bonds. 䉷 2000 Elsevier Science B.V. All rights reserved. Keywords: Polarisable continuum model; Solvation; Vibrational frequencies; Vibrational circular dichroism; Nuclear magnetic resonance chemical shifts
1. Introduction The polarisable continuum (PCM) solvation model [1] is widely used to accurately describe the solvent effect on the energy and on conformational barriers of various solutes, which do not undergo chemical transformation in solution. Recently, owing to the implementation of analytical first and second derivatives of the free energy [2–4], the calculation of various molecular properties has been made feasible, for instance vibrational/ Raman scattering frequencies and intensities (hyper) polarisabilities, nuclear shielding, ESR g-tensors and Vibrational circular dichroism (VCD) [5–7]. In this work, we present one of the first systematic applications of some of these new features to the compu* Corresponding author. Tel.: ⫹ 39-50-918244; fax: ⫹ 39-50502270. E-mail address: [email protected] (J. Tomasi).
tation of molecular properties in aqueous medium, selecting the two simplest molecules belonging to the class of amino acids, namely glycine and alanine, which show specific features in water solution, requiring a reappraisal of the solvation energy. The importance of amino acids in chemistry is well established [8]. As a matter of fact, they are among the most important biological molecules, since they constitute the “building blocks” of peptides and proteins. Most of their outstanding properties are due to their amphoteric nature, which derives from the presence, in their structure, of two functional groups of different polarity (the amino, NH2 and the carboxyl, COOH, groups). They are responsible, for example, of the existence of a tautomeric equilibrium between neutral and zwitterionic forms: in the gas phase (in case that the molecule vaporises without decomposing) the neutral form is energetically favoured, while in condensed phases, such as solutions in polar solvents and crystals, the ionic
0166-1280/00/$ - see front matter 䉷 2000 Elsevier Science B.V. All rights reserved. PII: S0166-128 0(00)00390-0
114
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
form is more stable [9,10]. These functional groups also play a fundamental role in the formation of hydrogen bonds, both with polar solvents and with other molecules in crystals; the presence of such interactions can be revealed through infrared (IR) or Raman spectroscopy and by nuclear magnetic resonance (NMR) experiments. Furthermore, all amino acids (except non-deuterated glycine) are chiral molecules and are naturally found in enantiomerically pure form (l). They exhibit, therefore, circular dichroism spectra. In particular, VCD has been extensively investigated recently, both experimentally and theoretically [11–16]. The choice of glycine and alanine, which, as stated before, are zwitterionic in water solutions, is rather challenging from the theoretical point of view because it is well known [17,18] that the definition of the cavity, which is at the basis of PCM model (see methodological section), turns out to be particularly crucial to get quantitative results on solvation energy in the case of anions and cations. Scope of this work is to verify and discuss how our model behaves in this particular case and to confront our findings with recent ab initio calculations that followed the more computationally expensive “supramolecular” approach [16].
2. Methodology and computational details The PCM [1,19] belongs to the category of solvation methods exploiting a continuum description of the solvent [20]. In brief, the physical system is represented by a charge distribution which describes the target molecule (or group of molecules) one has identified as solute. This charge occupies a given volume of space, called the molecular cavity, which is surrounded by an infinite continuum dielectric (the solvent) characterised by specific macroscopic properties (density, dielectric constant, refractive index, etc.). The interaction between the solute charge, or better the field it produces into the external volume occupied by the solvent, and the solvent itself, here described as a polarisable dielectric, is represented in terms of an apparent surface charge, s , spreading on the cavity surface. As a response, such charge produces a new field, called reaction field, which can modify the
previous state of the solute charge distribution (if the latter is not assumed to be composed of fixed point charges). The combined action of the solute field and the solvent reaction field can be analytically treated so as to obtain a final situation in which both solute and solvent distribution charges are mutually equilibrated. If this general picture is translated into a quantum mechanical language, the target solute is represented by its wavefunction and the interactions with the solvent are introduced in the Hamiltonian through specific perturbation, or reaction, operators. The Schro¨dinger equation to be solved becomes H 0 ⫹ V reac 兩C典 E兩C典
1
where H 0 is the Hamiltonian in the absence of the solvent. The solvent operator V reac acting on C is defined in terms of the surface apparent charge s , and thus, as said before, depends on the solute charge distribution (i.e. the solute wavefunction C ). This mutual interactions between C and V reac assures that the solution of Eq. (1) represents an equilibrated solute–solvent system. Here, we do not report any detail about the operators appearing in Eq. (1), and the procedures used to solve it, but we only recall that the solvent operator is defined as an integral over the cavity surface, to be solved analytically without any approximation. In the computational practice, the form of such integral, as well as of the whole apparent surface charge s , is simplified by discretising the cavity surface into small regions, called tesserae, of known area, ak. At each tessera k we can associate an apparent point charge, q k ak s , computed assuming that the surface charge s is constant in this small region. This discretisation allows to largely simplify all the operators involving the solvent. The main advantage of the PCM method is its almost immediate generalisability to different levels of the quantum mechanical description and to more complex physical systems involving, for example, external electromagnetic fields. In addition, the large freedom in the definition of the molecular cavity, that here has not to be limited to simple spherical or elliptical forms, but can be modelled according to the real shape of the solute geometry, allows a by far more realistic description of the solute–solvent system.
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
Thanks to all these peculiarities, in the last years the PCM model has been extended to the calculation of different molecular properties so as to make the solvated system a quantum mechanical problem, which can be studied with all the standard tools, originally defined for isolated molecules. This has first required the formulation of analytical gradients [2,3] and Hessians [4] taking into account the solvent effects, in order to allow the search, and the following characterisation, of all the interesting points on the potential energy surface of the system; geometry optimisations, definition of transition states, and vibrational analyses have then become available procedures exactly as in vacuo. Once this preliminary step is solved, the attention has been then shifted to more complex analyses, in particular, to the study of the combined action of the solvent and externally applied electric [21–24] and/or magnetic fields [25–28]. This research field is in fact characterised by numerous important applications, suffice it to cite the various spectroscopies, and thus it represents a big challenge for any theoretical solvation model. In the following analysis, besides to vibrational analysis, two specific external fields, and the related spectroscopies, will be considered; namely, an external magnetic field (to evaluate the nuclear shielding tensors determined in NMR measurements), and the combination of electric and magnetic fields (exploited in the VCD spectroscopy. To go a step further in the analysis of the molecular response properties of PCM solutes which are related to these processes, it becomes compulsory to define the free energy functional G. It is in fact this energetic quantity which has to be minimised with respect to the variational parameters defining the solute wavefunction so as to obtain the solution of Eq. (1); namely we can write G H0 ⫹
1 reac V 2
2
where the factor 1/2 takes into account the dependence of the solvent reaction operator on the solute wavefunction. For a PCM molecular solute G is the fundamental energetic quantity, as it univocally determines the behaviour of the molecular system in the presence of any internal and external perturbation. Thus, the
115
nuclear magnetic shielding which can be defined as the response of each nucleus X of the molecule to the contemporaneous effect of an external magnetic field and the nuclear magnetic moment, in the presence of the solvent becomes [24] X s ab
s 2G 2B a 2mXb
3
where Ba and mXb a; b x; y; z are the Cartesian components of the external magnetic field B, and of the nuclear magnetic moment of nucleus X. As in vacuo, also in condensed phases, the presence of a magnetic field introduces the problem of the definition of the origin of the related vector potential (historically called “gauge origin”). In ab initio calculations these problems can be dealt with introducing field dependent atomic orbitals (the so called GIAOs) [26–28]. The GIAO method has been recently applied to the PCM scheme [5], in order to obtain reliable nuclear shieldings for molecules in solutions. Still related to the GIAO approach, another recent extension of PCM is represented by the computation of VCD intensities for chiral molecules in solution [6]. In this case the external perturbation to consider is both magnetic and electric and the final effect is that the chiral sample responses to the left and right circularly polarised light in a different way, namely X Ri fi n ; n i 4 D1 n kn i
D1 (n ) is the differential absorption between the two differently polarised beams, n i and Ri are the excitation frequencies and the rotational strengths of the fundamental transition (0–1) in the ith normal mode, and the normalised band shape, fi, is assumed to be Lorentzian. Ab initio calculations of the frequencies and the rotational strengths are carried out within the harmonic approximation using the equation: Ri Imb具0兩mel 兩1典i ·具0兩mmag 兩1典c
5
where m el and m mag are the electric and the magnetic transition dipoles, respectively. These two quantities can be obtained by applying standard derivatives techniques (being the normal coordinates and the external magnetic field of the derivative parameters) that now, as we have shown
116
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
3. Results 3.1. Calculation of thermodynamic quantities for neutral and zwitterionic forms
Scheme 1.
above, are available also in the presence of a PCM solvent. In the present paper, we shall exploit a revised version of PCM called Integral Equation Formalism (IEF) [29–31] which represents a development of the original formulation [1] toward a more general and computationally efficient solvation model. In particular, in the following numerical study we shall report results obtained with the density functional extension of the PCM-IEF method, here implemented in a development version of Gaussian code [32]. All the features recalled above, geometry optimisation, vibrational analysis, nuclear shieldings and determination of VCD spectra, will be obtained in this framework by exploiting the Becke-3-LYP (B3LYP) hybrid functional with the 6-31G(d) basis set. As concerns the specificities of the solvation model, a cavity defined accordingly to the United Atom Topological Model (UATM) [33] will be exploited. In brief, this cavity is formed by interlocking spheres centred on selected nuclei (in our case C, N, and O) with radii defined according to the topological state of each nucleus (i.e. number and type of atoms bonded). The construction of the cavity and its tessellation are performed with the program gepol. [34] The solvent will be represented by a dielectric constant of 78.39 corresponding to the experimental value for liquid water at room temperature. To correctly assign normal modes of vibration, we performed a potential energy distribution (PED) calculation with the program gar2ped [35]. Table 1 Calculated free energy differences (DG) for selected amino acids processes (values in kcal/mol) Process
Glycine
Alanine
(1) NTaq ! ZWaq (2) NTg ! NTaq (3) NTaq ! ZWaq
⫺3.42 ⫺10.11 ⫺13.52
⫺4.48 ⫺7.42 ⫺11.90
Within the PCM framework, several thermochemical quantities can be calculated for amino acid systems. The processes we have selected may be summarised in the following Scheme 1: (the two symbols will identify the two tautomeric forms from now on). Core of any thermochemical calculation is the free energy of the system, since it is related (at constant pressure) to the equilibrium constant of a given chemical process and as it can be used to obtain other quantities through standard thermodynamics machinery. The free energy of a system in solution [20] may be decomposed into: electrostatic contribution Gel, cavitation term (related to the work necessary to create the cavity in the solvent) Gcav, dispersion free energy Gdis (due to dispersion interactions between solute and solvent), repulsion free energy Grep and a term due to molecular motions (thermal motions), GMm. The first term is provided by the PCM single point calculation (as explained in the methodological section); the Gcav, Gdis and Grep terms (often referred to as non-electrostatic contributions) are evaluated by exploiting semi-empirical relationships [36,37]; the last contribution (GMm) can be calculated applying the standard methods of statistical mechanics, once vibrational frequencies and moments of inertia for the system are known. The total free energy can therefore be expressed as follows: G Gel ⫹ Gcav ⫹ Gdis ⫹ Grep ⫺ RT ln qtr qrot qvib 6 In this formula, the subscripts identify the various components into which the molecular partition function can be factorised according to the Born–Oppenheimer approximation. Notice that the zero-point energy (ZPE) is included in the vibrational term. The computed results for all the processes indicated in Scheme 1 are collected in Table 1; the various contributions to G are listed in Table 2. The analysis of such data will be given in the following sections.
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127 Table 2 Calculated values of the various components of solvation free energy (kcal/mol)
Electrostatic (Gel)
Glycine NT
Glycine ZW
Alanine NT
Alanine ZW
⫺12.11
⫺14.15
⫺10.77
⫺13.94
(Values calculated with respect to neutral form optimised in the gas phase) Total non-electrostatic 1.41 ⫺0.14 2.13 0.63 Cavitation (Gcav) 11.83 10.45 14.31 12.99 Dispersion (Gdis) ⫺14.48 ⫺15.48 ⫺15.74 ⫺16.62 Repulsion (Grep) 4.06 4.89 3.56 4.25 Thermal motions at 32.02 32.19 48.38 48.58 298.15 K (GMm)
3.1.1. Calculation of D G for process 1 (NTaq ! ZWaq) As it was stated previously, amino acids are found to exist mainly in zwitterionic form in neutral water solutions. Unfortunately, the mole fraction of the neutral form in such media is so negligible that it cannot be revealed through experimental measurements. Therefore, the equilibrium constant of process 1 in solution is not known and thermodynamic quantities are not accessible. Nevertheless, they can be estimated. A possible approach was proposed by Haberfield [38], who suggested that, in most cases, thermodynamic quantities for a given compound (i.e. the free energy of formation, DfG) could be calculated from the values which the same property assumes in the “fragments” in which the molecule can be divided (“summation method”). Following this guideline, it is possible to devise a thermodynamic cycle which allows one to obtain DGZW-NT if gas phase acidities are known [39]. We report it briefly below; the reader is referred to the original article for further details. DGZW-NT, corresponding to step A, is equal to ⫺B ⫹ C ⫹ D ⫹ E ⫺ F (Scheme 2). Some of these quantities can be obtained by experimental measurements, others are estimated as
Scheme 2.
117
mentioned before. By following this procedure, Haberfield reported a value of ⫺7.67 kcal/mol for DGZW-NT of glycine. Another kind of estimate, which is performed through several equilibrium constant measurements, is based on the assumption that the esterification of the carboxylic group has limited effect on the basicity of the amino group [40]. In this case, DGZW-NT ⫺7:27 kcal=mol: Despite the results obtained by applying the “summation method” have been argued about [41] since it seems that the method induces significant errors (as large as 5–10 kcal/mol), we chose to refer our discussion to Haberfield’s findings because that procedure can be readily modified to obtain DGZW-NT values for every amino acid, provided that its gas phase acidity is known [39]. Moreover, in the case of glycine, the two approaches give very similar results and, as far as we know, the second method has not been applied to other amino acids yet. It should be clear, at this point, that estimation data in our possession are not very reliable as a test to check the validity of our computational results. We are aware that both “experimental” and theoretical results undergo large errors. Nonetheless, our aim is first of all to verify that zwitterions are stable molecules within the PCM computational scheme, so as to be able to use the geometries obtained theoretically for the evaluation of other properties, such as vibrational frequencies and chemical shifts. In the first step, the geometry of ZW and NT forms in water solution was optimised at the B3-LYP/631G(d) level of theory with the IEF version of PCM (see methodological section). The potential energy surface (PES) of neutral glycine and alanine has been studied in detail by several authors [42–46]. In spite of their apparent simplicity, these two molecules give rise in vacuo to quite complex potential surfaces which include several minima and transition states connecting them. In both cases, in the lowest energy structure, the two hydrogen atoms of the amino group lean towards a lone pair of carboxyl oxygen; the structure for glycine has Cs symmetry (see Fig. 1). The ZW case is rather more complex [47–49]. In fact, the zwitterionic form does not seem to exist in the gas phase as a stable structure; explicit water molecules [48], a crystal field or a continuum solvent distribution [36] is
118
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
Fig. 1. Optimised geometries of glycine and S-alanine B3-LYP/6-31G(d): left, neutral forms; right, Zwitterionic forms.
needed to stabilise it. Ding et al. [50] performed a detailed study of glycine ZW PES at the RHF/631G ⫹ (d) level of theory, in which they were able to characterise a local minimum with Cs symmetry, ˚ ) N–H···O hydrogen with a particularly short (1.55 A bond; this structure underwent proton transfer when the basis set was enlarged to Dunning’s DZP, thus moving to the lower energy neutral structure. Starting from that geometry, Jensen et al., in a more recent investigation, found that two water molecules are necessary to prevent the proton transfer at the MP2/ DZP⫹⫹//RHF/DZP level of theory [48]. As for alanine ZW, Sambrano et al. [49] recently dealt with intramolecular proton transfer in aqueous medium within the SCRF (self-consistent reaction field) approach by using an ellipsoidal cavity and multipolar expansion. The inclusion of both diffuse functions and electron correlation was deemed necessary to make the ZW a global minimum. Bearing in mind these apparently contradictory results, we chose in vacuo-type geometries (according
to Ding) as starting structures to perform our analysis, by applying the PCM method. For the neutral forms, we took the aforementioned lowest energy structures. The most significant optimised geometrical parameters for the two systems in water are reported in Table 3. Two features must be noted: the introduction of a methyl group in the side chain of the NT glycine induces a deviation from planarity of about 15⬚ in the O–C–C–N backbone; a similar distortion can be found passing to ZW too, owing to the formation of a ˚ long in glycine hydrogen bond. The latter is 1.834 A ˚ ZW and 1.800 A in alanine ZW. Such optimised geometries have been used for the vibrational analysis needed to calculate the thermal contributions to the free energies. As concerns this point some preliminary considerations have to be done. It is well known from statistical mechanics theory [51] that certain torsional low-frequency motions cannot be regarded as “vibrations” but should rather
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
119
Table 3 Relevant geometrical parameters for NT and ZW forms in water (see Fig. 1 for labels) a Parameter (NT) d N5–H6 d N5–H7 d C4–N5 d C1–C4 d C1–O2 d C1–O3 d O2–H8 a C1–C4–N5 d O3–C1–N4–N5 d H6–N5–C4–C1 d H7–N5–C4–C1 a
Glycine (Cs)
Alanine
1.019 1.019
1.020 1.021
1.452 1.525 1.355 1.211 0.976 112.62 0. 56.50 ⫺56.50
1.463 1.532 1.344 1.219 0.991 112.62 15.78 59.20 ⫺54.31
Parameter (ZW)
Glycine
Alanine
d N5–H6 d N5–H7 d N5–H8 d C4–N5 d C1–C4 d C1–O2 d C1–O3 d O3–H8 a C1–C4–N5 d O3–C1–N4–N5 d H6–N5–C4–C1 d H7–N5–C4–C1
1.030 1.029 1.053 1.503 1.553 1.247 1.265 1.834 102.59 8.53 ⫺12.67 103.31
1.030 1.030 1.054 1.515 1.561 1.247 1.266 1.800 105.60 ⫺9.13 12.52 ⫺103.51
˚ , angles in degrees. Bond lengths in A
be considered as hindered rotations. In polyatomic non-linear molecules, in fact, individual molecular groups may rotate relative to each other across the bond linking them. A barrier to internal rotation UIRMAX therefore exists. According to the magnitude of the barrier with respect to RT, the motion can be classified as free rotation UIR⫺MAX =RT 艑 1; hindered rotation UIR⫺MAX =RT 5–10 or ordinary normal mode of vibration UIR⫺MAX =RT ⬎ 10: The subject was studied by Pitzer and Gwinn [52], who compiled accurate tables containing the corrections to thermodynamic functions (the increase of the free energy from the free rotor value is reported) for the single hindered rotor. An extension of Pitzer procedure to multiple rotors, assuming negligible coupling between them, has recently been carried out by Ayala et al. [54] and it is now available in the latest version of the Gaussian package [32]. The procedure automatically identifies hindered rotation motions among the calculated normal modes and yields the necessary corrections to partition functions, by using simple expressions (of Lennard-Jones type) for the rotation potentials. Two hindered rotations were found for glycine, and three for alanine; in either case, they correspond to rotations around C–C and C–N bonds, whose harmonic frequencies occur between 70 and 320 cm ⫺1 (this result was predictable before running the calculation). Approximately the same result is obtained when the contribution for hindered rotation is evaluated through direct calculation of the various UIR-MAX by
performing a torsional energy scan at the same level of theory and then using Pitzer tables. The difference between the final computed DG is about 0.5 kcal/mol for alanine; we therefore decided to use the Ayala method. Still another contribution to free energy must be evaluated. In fact, even though the solvent is modelled as a continuum, it is made up of individual molecules, characterised by their internal and external motions. While moving, they can collide with the solute and exert a friction upon it; the solute–solvent interaction itself also plays a similar role. This kind of “medium hindrance” [53] can be described using Pitzer model as well. The increase in free energy was calculated to be about 2.3 kcal/mol by considering UIR⫺MAX =RT 5: By summing all the contributions to G, we found a value of ⫺3.42 kcal/mol for glycine DGZW-NT and ⫺4.48 kcal/mol for alanine (see Table 1). The DGZW-NT values reported in Table 1, although not in excellent agreement with “experimental” results (but see the former discussion upon the reliability of them), are able to correctly predict ZW stability without requiring explicit water molecules and to reproduce suitably the difference between alanine and glycine. In fact, the difference between the DGZW–NT of the molecules is equal to 0.68 kcal/mol from the estimates and 1.06 kcal/mol from our calculation. We intend to perform calculations on other amino acids to verify if the trend is respected and then calibrate our method against systematic errors that, according to us, may occur on both sides.
120
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
3.1.2. Calculation of D G for processes 2 (NTg ! ZWaq) and 3 (NTg ! NTaq) Gaffney et al. [55] performed an estimate of the enthalpy difference (DH) for the process (2) by calculating DHdiss ⫺ DHsub, the solution and sublimation heats of solid glycine, respectively, reporting a value of ⫺19.2 ^ 1 kcal/mol at 298 K (note that the uncertainty regards the sublimation heat only, while a relatively higher, but not quantifiable error, is to be associated with the solution heat). The calculation of this quantity within the PCM procedure was studied in the past by Bonaccorsi et al. [47,56], who appraised the purely electrostatic free energy (Gel) difference to be ⫺16.5 kcal/mol at the RHF/4-31G ⴱ level of theory; to obtain the required DH they evaluated the electrostatic entropic correction to G in solution with the relation 2 DG 8 DH DG ⫺ T 2T P Assuming that the temperature change of electrostatic DG is given only by modification of the dielectric constant 1 and of the cavity volume v, we obtain: 2 dv 2 d1 DG ⫺T DG DH DG ⫺ T 2v 21 P dT P dT 9 The values in brackets are calculated analytically, while the derivatives of v and 1 are taken from experimental data. They reported a value of 3.0 kcal/mol (see the original paper for details) for such correction, and ⫺19.5 kcal/mol for DH. The authors regarded the almost perfect numerical coincidence as due to error compensation, both from the theoretical point of view (neglect of thermal contributions to the free energy) and from the experimental one (the already cited large errors on solution enthalpy). We performed a free energy calculation for the neutral forms of glycine and alanine in gas phase, by optimising the geometry of the species and by subsequently calculating thermal contributions to G, as discussed previously for process (1). From the G values thus obtained for the gas phase species, we were able to appraise DG for processes (2) and (3). Similar calculations were also performed for alanine. A summary of all the quantities calculated is reported in Table 1.
A value of ⫺13.52 kcal/mol was calculated for DG of process 2 in the case of glycine. We also re-evaluated the electrostatic contribution (9), getting a value of ⫺4.07 kcal/mol; we computed, as well, the thermal motions entropy contribution from the calculated partition functions (⫺1.0 kcal/mol). Therefore, our estimate for DH of process (3) is ⫺18.57 kcal/mol. 3.2. Calculation of IR and VCD spectra The second subject we dealt with is the prediction of vibrational and VCD spectra of alanine and glycine in aqueous environment. As already cited in Section 1, the calculation of vibrational frequencies was performed at B3-LYP/6-31G (d) level on the geometries optimised with the same method. All the models were carried out following the harmonic approach and considering the solvent properties and the cavity unperturbed; moreover, we assume that the solvent is always equilibrated to the solute and that this equilibrium is preserved during vibrations. This hypothesis is reasonable for normal modes of intermediate frequency (bending region) while it may be unsatisfactory for high-frequency (stretching) and low-frequency (torsional) motions. In the former case, the solvent may not be able to attain equilibration during the solute motion; in the second case, for which large amplitudes are possible, there might be effects on the cavity size and the orientational motions of the solvent may perturb the solute displacements. In this study, we did not take into account any of these effects. Further details on the calculation of vibrational frequencies in solution can be found in a general review [20] and in a recent paper by our group [61]. 3.2.1. Vibrational spectra Among all the calculations we ran, we choose to report the aqueous ZW spectra only, because it is they that are worth discussing to assess the behaviour of the PCM method. The results are reported in Tables 4 and 5, together with the experimental data available and with a set of theoretical results obtained by Alper et al. [57,58] (glycine) and by Tajkhorshid et al. [16] for alanine. Experimental measurements go through several difficulties in water solution because of water absorptions interference (at least at the time the data we
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127 Table 4 Vibrational frequencies of Glycine ZW in water (cm ⫺1) PED HN–H2 ant. Str. HN–H2 sym. Str. C–H10 str. H2N–H str. ⫹ C–H9 str CH2 sym. str ⫹ H2N–H str. CO2⫺ ant. str. NH3 ant. bend. NH3 ant. bend. CH2 sciss. NH3 sym. bend. CH2 rock ⫹ NH3 sym. bend. CH2 wag. ⫹ NH7 bend. CH2 twist. ⫹ NH8 bend. CH2 wag. ⫹ CCN bend. HCC ⫹ HNC bend. C–N str. ⫹ CCN bend HN–H2 bend. ⫹ CHN bend. CH2 wag. ⫹ C–C str. ⫹ CO2 bend. CCN ⫹ CCH bend. Out of plane CO2 bend. ⫹ CCH bend. CCN ⫹ CCO bend. CCN bend. NH3 rotation NH3 and CH2 rotation a b c
Aa
Bb
Cc
3152 3058 2962 2968 1634 1615 1597 1533 1440 1410 1327 1315 1118 1119 1027 919 896
3397 3345 3126 3021 3011 1735 1659 1617 1454 1443 1388 1339 1289 1104 1086 983 906 868
3390 3366 3288 2924 2868 1831 1790 1700 1651 1547 1457 1437 1379 1257 1196 1032 1015 892
665 613
680 560
781 701
577 502 364 238
497 321 262 73
661 569 408 234
Takeda et al. [59], Ghazanfar et al. [60]. This work. Alper et al. [57,58]; scaled values.
possess [59,60] were collected); therefore, comparison with experiment is rather hampered. Whenever there are uncertainties in the datum, crystalline phase frequencies have been reported, as in other theoretical works [57,58]. The experimental spectrum of alanine was assigned by Diem et al. [11] by means of isotope substitution. To correctly assign the calculated frequencies, besides a graphical animation of vibrational modes with Gaussview package (included in Gaussian distribution [32]), a convenient way we encountered was to perform a PED calculation. We made it by employing the program gar2ped, a Gaussian output postprocessing utility written by Martin and Van Alsenoy [35] which uses a valence internal co-ordinate scheme proposed by Pulay and yields the fractional contribution of each defined co-ordinate within the calculated normal mode. Our assignment is reported in Tables 4
121
and 5 as well (only the major contributions to the mode are reported). As it can be seen, there is a fairly good agreement between calculated and experimental data (it is worth noticing that no scaling of the calculated values was applied). Yet, a couple of remarks are needed. We will limit our discussion to alanine spectrum, because only in this case an experimental assignment of normal modes is available (we consider it at all reliable, having been carried out by isotope substitution). First, the quality of the agreement is poorer for N– H stretching modes. In the experimental spectra, in fact, the three highest frequency modes (3080, 3060 and 3020 cm ⫺1) are attributed to the stretching vibrations of NH⫹ 3 group (two antisymmetrical and one symmetrical, respectively). The shift in frequency between these modes and the neutral gaseous form N– H stretching vibrations (occurring at about 3560 cm ⫺1[46]) is 500 cm ⫺1 on average. In our calculation, we obtained two N–H stretching frequencies in the upper part of the spectrum (3391 and 3338 cm ⫺1), which the PED calculation assigned to antisymmetrical and symmetrical stretching motions, respectively, of the two N–H bonds not involved in the intramolecular hydrogen bond (N5–H6 and N5–H7) and one N–H stretching frequency at 2992 cm ⫺1, which was assigned to a pure N5–H8 (that is, the H-bond acceptor fragment) stretching mode. On the contrary, the corresponding sixth normal mode (in decreasing order of frequency) of the experimental spectrum, which absorbs at 2962 cm ⫺1, was proved to be a methine stretching vibration. Therefore, it seems that the continuum alone is not able to account for such strong interactions with water (responsible of a shift of 500 wavenumbers, while only 150–200 are supplied by PCM). Moreover, in this case, the proximity of NH⫹ 3 and CO2 – groups is crucial. The resulting intramolecular H-bond (O3–H8) is much more effective than the interaction with the continuum solvent, even if the latter has a very high dielectric constant. This overestimate of a single N–H with respect to all the others is responsible for the permutation of the frequency order as well. We also performed a test calculation on a small cluster made up of an amino acid molecule and one water molecule in the vicinity of NH⫹ 3 and CO2 – groups; all the system was surrounded by a continuum
122
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
Table 5 A–C: Vibrational frequencies (cm ⫺1); D–E: VCD rotational strengths (esu 2 cm 2 × 10 44) for S-Alanine ZW in water PED a
Ab
Ba
Cc
Da
HN–H2 ant. str. HN–H2 sym. str. CH3 ant. str. CH3 ant. str. CH3 sym. str. H2N–H str. d Cm –H str. d CO2 ant. str. NH⫹ 3 bend. NH⫹ 3 bend. CH3 sym. bend. CH3 ant. bend. CH3 bend. ⫹ NH⫹ 3 bend. CH3 bend. ⫹ NH⫹ 3 bend. CH3 bend. ⫹ Cm –H bend. ⫹ CO2 s str. Cm –H bend. Cm –H bend. NH2 rock. ⫹ CH3 wag. C–N str. ⫹ CH3 bend. H2N–H bend. ⫹ Cm –H bend. NH⫹ 3 bend. ⫹ CH3 twist. HN–H2 bend. C–N str. ⫹ C–C str. CO2 bend. ⫹ OCC bend. Out of plane CO2 bend. C1–C–N bend. ⫹ C9–C–N bend. C–C–O bend. ⫹ C–C str. ⫹ C9–C–N bend. C–C–O bend. ⫹ C1–C–N bend. C–N tors. ⫹ C–CO2 tors. C–N tors. ⫹ C–CH3 tors. C–CH3 tors. CO2 –CH3 tors. ⫹ C–N tors.
3080 3060 3020 3003 2993 2962 2949 1645 1625 1607 1498 1459 1459 1410 1375
3391 3338 3149 3121 3051 2992 2952 1731 1664 1623 1520 1516 1445 1422 1388
3164 3142 3137 3135 3113 3093 3064 1759 1750 1695 1635 1536 1529 1441 1416
⫺5 19 0 0 ⫺3 34 ⫺20 ⫺36 42 ⫺100 ⫺11 ⫺2 ⫺19 5 104
⫺11 ⫺4 ⫺71 37 ⫺30 56 ⫺2 ⫺17 27 ⫺53 ⫺9 0 ⫺9 68 48
1351 1301 1220 1145 1110 1001 995 922 850 775 640
1364 1297 1213 1105 1084 1002 984 881 827 759 645
1390 1350 1290 1213 1133 1062 1028 918 840 775 635
⫺73 2 25 43 ⫺87 ⫺1 24 ⫺42 5 23 26
⫺171 38 50 1 ⫺11 3 ⫺2 ⫺22 20 22 ⫺8
527 477 399
520 381 342
625 528 432
⫺9 ⫺11 ⫺25
25 ⫺20 22
296 283 219 184
283 263 240 70
355 290 264 170
38 71 ⫺25 ⫺10
⫺33 50 ⫺30 ⫺9
a b c d
Ec
Present work. Diem et al. [11–13]. Tajkhorshid et al. [16]. See text.
distribution of water. The intramolecular · ··CO –H–bond was replaced by an interNH⫹ 2 3 molecular N–H· ··OH2 interaction, but the calculated frequency shift was only a little modified (2992– 2943 cm ⫺1). This result is compliant with recent findings by our group [61] about ROH···OR 0 interactions in water and methanol. In order to reproduce the observed frequency shifts, every H-bond site should
likely be allowed to interact directly with an explicit water molecule. In the rest of the spectrum, instead, the agreement of our frequencies is better than those calculated with the “super-molecular” approach [16]. Therefore, it seems that whenever the interaction is less specific, a continuum model is rather more efficient.
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
3.2.2. VCD spectra As it was already anticipated in the methodological section, the newly developed versions of PCM have recently made the calculation of several molecular response properties in solution feasible. Among them, VCD plays a fundamental role for the significant impact it has in the study of biomolecules. Amino acids, being the most easily available chiral compounds, have been the subject of extensive research for a couple of decades, both experimental and (more recently) theoretical [11–16]. Some amino acids VCD spectra show recurrent patterns (i.e. progression of positive and negative bands). One of the most important patterns can be found in the methine bending absorptions (between 1250 and 1450 cm ⫺1), which show a typical ( ⫹ , ⫺ , ⫹ ) pattern. It can be considered as a sort of “fingerprint” for the smallest amino acids (from deuterated glycine to proline) and for their polypeptides. Ab initio (or DFT) calculation of VCD spectra is a very delicate matter. In order to obtain good results, both for sign and for band intensities, optimisation and subsequent frequency calculation employing very large basis sets is recommended [62–64]. Such a design is hardly applicable in our case, since the stability of the zwitterionic form is largely dependent on the basis set and on the method used. Moreover, geometry convergence in these systems is already very difficult when a medium size basis is used, becoming almost prohibitive when the set is enlarged. Therefore, in the present work, we chose to calculate VCD spectra for l-alanine (corresponding to absolute S configuration) and S-Glycine-Cd1 at B3-LYP/6-31G (d) level of theory as well. Our findings are reported in Table 5 and Fig. 2a and b. In the former table, we confront our results with those of Tajkhorshid et al. [16], who employed a supermolecular model, while in Fig. 2a and b, the calculated VCD spectra of the two molecules in the methine bending region are drawn. As it can be seen (Table 5), the VCD spectra of lalanine yielded by the PCM calculation, is quite similar to that obtainable from a supermolecular model, with the exception of few vibrations (see also Section 3.2.1). Moreover, the spectra for both amino acids
123
(Fig. 2a and b) in the methine bending region reflects the experimental trend fairly well. 3.3. NMR chemical shift calculations Chemical shift calculations, in addition to being a theoretically interesting subject, may constitute nowadays a valid aid for complex NMR spectra interpretation. Recent works have shown that ab initio or density function calculations are able to give satisfactory results even for significantly large molecules, as taxol [65]. A great deal of experimental measurements is performed in solution; therefore, prevision of solvent effect is crucial. As it was mentioned in the methodological section, the newly implemented features of PCM allow one to calculate several properties, including nuclear shielding tensors. In this case, we decided to perform several calculations of shielding properties at the optimised B3-LYP/ 6-31G(d) geometry, employing different methods and basis sets. The reason for such a choice is that shielding properties are known to depend loosely on the geometry of the system and more tightly on the basis sets used; the quantum mechanical method applied also seems to play a role, but it has not been well established yet. We evaluated chemical shifts for B3-LYP/6-31G(d) optimised geometries by calculating the isotropic nuclear shielding of C and N atoms of glycine and alanine at the following levels of theory, namely HF/6-31G(d), HF/6-311⫹G(d,p), HF/6-311⫹G(2d,p), B3-LYP/6-31G(d), B3-LYP/6311⫹G(d,p), B3-LYP/6-311⫹G(2d,p) and then by subtracting these values from those calculated, at the same level of theory, for C atoms of neat tetramethylsilane and N atoms of neat nitromethane (data in our possession are referred to these standards). The values for neat liquids were obtained optimising the molecular geometry of TMS and CH3NO2 (at B3-LYP/6-31G(d) level), each embedded in the like liquid, described as a continuum according to the PCM scheme and then calculating the isotropic shielding with the required method. The calculated 13C and 15N chemical shift values for glycine and alanine in water solution are reported in Table 6, together with experimental values [66,67]. As it can be seen, in the case of 13C shifts, the best agreement was found using Hartree–Fock theory at the same basis set (smallest) used for geometry
124
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
Fig. 2. (a) Calculated VCD spectrum of S-Glycine-C-d1 in the methine bending region. (b) Calculated VCD spectrum of l-alanine in the methine bending region.
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
125
Table 6 13 C and 15N chemical shifts of glycine and alanine in water (in ppm with respect to TMS and nitromethane) Glycine
Expt. a
Calc.
Method
Alanine
Expt. a
Calc.
Method
C(yO)
173.2
171.1596 184.5958 184.1293 157.8906 178.1734 178.2840
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
C(yO)
176.5
173.6479 187.3207 186.7785 161.6340 182.1343 182.0875
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
C(–N)
42.2
42.8664 50.7837 51.4952 41.9151 47.1808 46.5477
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
C(–N)
51.3
50.4673 53.2342 52.5814 52.6163 59.0460 58.5760
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
350.2
468.6218 494.2359 488.4879 348.8527 380.1944 379.2182
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
Ca
17.0
16.9199 18.3043 17.7558 18.7769 20.6157 20.0401
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
458.2070 483.4378 477.7165 337.5204 367.7613 367.0281
HF/6-31G(d) HF/6-311⫹G(d,p) HF/6-311⫹G(2d,p) B3LYP/6-31G(d) B3LYP/6-311⫹G(d,p) B3LYP/6-311⫹G(2d,p)
N
N
a
338
Experimental data from [66,67].
optimisation at DFT level. The enlargement of the basis and the use of density functional did not lead to any improvement. When using DFT, the set augmentation seems to help only in the carbonyl case. The situation is reversed for 15N: Hartree–Fock methods yield poor values while B3-LYP behaves much better; even in this case, a larger basis set leads to worse results. The relative behaviour HF/ DFT with different nuclei has been already found to be unpredictable [65]. It is clear that our findings cannot easily be interpreted and we believe that the surprisingly good agreement of low-level calculations is due to a fortuitous cancellation of errors. Nevertheless, we can affirm that our calculated chemical shifts are, on average, in rather good agreement with the experiment, apart from 15N Hartree–Fock values.
4. Conclusions In this paper, we have described some of the capabilities of the PCM continuum solvation model, as well as some of its shortcomings. We have demonstrated that it is able to correctly reproduce solvent effects upon the energy of molecular systems, even in the “difficult” example of zwitterionic amino acids stability. A very good behaviour is found in the prediction of nuclear shielding properties as well, even though the results we obtained depend on the choice of the theoretical method and on the basis set. As regards vibrational properties and VCD, we have found that PCM yields a satisfactory description in all the cases in which the interactions between solute and solvent are not short-range and specific. In this case (i.e. the stretching modes of N–H fragments) a supermolecular model could give results in better agreement with experimental data.
126
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127
Comprehensively, we may state that the model is a convenient and versatile method to account for solvent effects. Acknowledgements Partial support of CNR is here acknowledged. References [1] S. Miertus, E. Scrocco, J. Tomasi, J. Chem. Phys. 55 (1981) 117. [2] E. Cances, B. Mennucci, J. Chem. Phys. 109 (1998) 249. [3] E. Cances, B. Mennucci, J. Tomasi, J. Chem. Phys. 109 (1998) 260. [4] B. Mennucci, R. Cammi, J. Tomasi, J. Chem. Phys. 110 (1999) 6858. [5] R. Cammi, B. Mennucci, J. Tomasi, J. Chem. Phys. 110 (1999) 9877. [6] B. Mennucci, R. Cammi, J. Tomasi, P. Stephens, J.R. Cheeseman, M.J. Frisch, in preparation. [7] N. Rega, M. Cossi, V. Barone, J. Am. Chem. Soc. 119 (1997) 12 962. [8] P.J. Greenstein, M. Winitz, Chemistry of the Amino Acids, Wiley, New York, 1961, pp. 435–522. [9] K. Ijima, K. Tanaka, S. Onuma, J. Mol. Struct. 246 (1991) 257. [10] J. Almlo¨f, A. Kuick, O.J. Thomas, J. Chem. Phys. 59 (1973) 3901. [11] M. Diem, P.L. Polavarapu, M. Oboodi, L.A. Nafie, J. Am. Chem. Soc. 104 (1982) 3329. [12] B.B. Lal, M. Diem, P.L. Polavarapu, M. Oboodi, T.B. Freedman, L.A. Nafie, J. Am. Chem. Soc. 104 (1982) 3337. [13] T.B. Freedman, M. Diem, P.L. Polavarapu, L.A. Nafie, J. Am. Chem. Soc. 104 (1982) 3343. [14] M. Diem, J. Am. Chem. Soc. 110 (1988) 6967. [15] W.M. Zuk, T.B. Freedman, L.A. Nafie, J. Phys. Chem. 93 (1989) 1771. [16] E. Tajkhorshid, K.J. Jalkanen, S. Suhai, J. Phys. Chem. B 102 (1998) 5899. [17] J. Mestres, M. Solar, M. Orozco, J. Phys. Chem. 100 (1996) 606. [18] F.J. Luque, M. Orozco, Chem. Phys. 182 (1994) 237. [19] R. Cammi, J. Tomasi, J. Comp. Chem. 16 (1995) 1449. [20] J. Tomasi, M. Persico, Chem. Rev. 94 (1994) 2027. [21] R. Cammi, M. Cossi, J. Tomasi, J. Chem. Phys. 104 (1996) 4611. [22] R. Cammi, M. Cossi, B. Mennucci, J. Tomasi, J. Chem. Phys. 105 (1996) 10 556. [23] R. Cammi, B. Mennucci, J. Tomasi, J. Am. Chem. Soc. 34 (1998) 8834. [24] R. Cammi, B. Mennucci, J. Tomasi, J. Phys. Chem. A102 (1998) 870. [25] R. Cammi, J. Chem. Phys. 109 (1998) 3185.
[26] H.F. Hameka, Rev. Mod. Phys. 34 (1962) 87. [27] R. Ditchfield, Mol. Phys. 27 (1974) 789. [28] K. Wolinski, J.F. Hinton, P. Pulay, J. Am. Chem. Soc. 112 (1990) 8251. [29] E. Cances, B. Mennucci, J. Math. Chem. 23 (1998) 309. [30] E. Cances, B. Mennucci, J. Tomasi, J. Chem. Phys. 107 (1997) 3031. [31] B. Mennucci, E. Cances, J. Tomasi, J. Phys. Chem. B 101 (1997) 10 506. [32] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery, R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam, A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, C.S. Clifford, J. Ochterski, G.A. Petersson, P.Y. Ayala, Q. Cui, K. Morokuma, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J.V. Ortiz, A.G. Baboul, J. Cioslowski, B.B. Stefanov, G. Liu, C.A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, C.M.W. Wong, J.L. Andres, M. Head-Gordon, E.S. Replogle, J.A. Pople, Gaussian, Inc., Pittsburgh, PA, 1998; gaussian 99, Development Version (Revision A.10 ⫹ ). [33] V. Barone, M. Cossi, J. Tomasi, J. Chem. Phys. 107 (1997) 3210. [34] J.L. Pascual-Ahuir, E. Silla, J. Tomasi, R. Bonaccorsi, J. Comp. Chem. 8 (1987) 778. [35] J.M.L. Martin, C. Van Alsenoy, GAR2PED, University of Antwerp, 1995. [36] F.M. Floris, J. Tomasi, J. Comput. Chem. 10 (1989) 616. [37] F.M. Floris, J. Tomasi, J.L. Pascual-Ahuir, J. Comput. Chem. 12 (1991) 784. [38] P. Haberfield, J. Chem. Ed. 57 (1980) 347. [39] R.A.J. O’Hair, J.H. Bowie, S. Gronert, Int. J. Mass Spectrom. Ion Processes 117 (1992) 23. [40] G. Wada, E. Tamura, M. Okina, M. Makamura, Bull. Chem. Soc Jpn. 55 (1982) 3064. [41] M.J. Locke, R.T. McIver Jr., J. Am. Chem. Soc. 105 (1983) 4226. [42] P. Palla, C. Petrongolo, J. Tomasi, J. Phys. Chem. 84 (1980) 435. [43] V. Barone, F. Lelj, C. Adamo, J. Chem. Phys. 102 (1995) 364. [44] D.T. Nguyen, A.C. Scheiner, J.W. Andzelm, S. Sirois, D.R. Salahub, A.T. Hagler, J. Comput. Chem. 18 (1997) 1609. [45] S.G. Stepanian, I.D. Reva, E.D. Radchenko, M.T.S. Rosado, M.L.T.S. Duarte, R. Fausto, L. Adamowicz, J. Phys. Chem. A102 (1998) 1041. [46] S.G. Stepanian, I.D. Reva, E.D. Radchenko, M.T.S. Rosado, J. Phys. Chem. A102 (1998) 4623. [47] R. Bonaccorsi, P. Palla, J. Tomasi, J. Am. Chem. Soc. 106 (1984) 1945. [48] J.H. Jensen, M.S. Gordon, J. Am. Chem. Soc. 117 (1995) 8159. [49] J.R. Sambrano, A.R. de Sousa, J.J. Queralt, J. Andre´s, E. Longo, Chem. Phys. Lett. 294 (1998) 1.
L. Gontrani et al. / Journal of Molecular Structure (Theochem) 500 (2000) 113–127 [50] Y. Ding, K. Krogh, Jespersen, Chem. Phys. Lett. 199 (1992) 261. [51] G.N. Lewis, M. Randall, Thermodynamics, 2nd ed., McGrawHill, New York, 1961 (revised by K.S. Pitzer, L. Brewer). [52] K.S. Pitzer, D.W. Gwinn, J. Chem. Phys. 10 (1942) 428. [53] O.N. Ventura, A. Lledos, R. Bonaccorsi, J. Bertran, J. Tomasi, Theor. Chim. Acta 72 (1987) 175. [54] P.Y. Ayala, H.B. Schlegel, J. Chem. Phys. 108 (1998) 2318. [55] J.S. Gaffney, R. Pierce, L. Friedman, J. Am. Chem. Soc. 99 (1977) 4293. [56] R. Bonaccorsi, E. Scrocco, J. Tomasi, J. Biosci. 8 (1985) 627 (Proceedings of the International Symposium on Biomolecular Structural Interactions, Suppl.) [57] J.S. Alper, H. Dothe, D.F. Coker, Chem. Phys. 153 (1991) 51. [58] J.S. Alper, H. Dothe, D.F. Coker, Chem. Phys. 161 (1992) 199. [59] M. Takeda, R.E.S. Iavazzo, D. Garfinkel, I.H. Scheinberg, J.T. Edsall, J. Am. Chem. Soc. 80 (1958) 3813.
127
[60] S.A.S. Ghazanfar, D.V. Myers, J.T. Edsall, J. Am. Chem. Soc. 86 (1964) 3439. [61] C. Cappelli, B. Mennucci, C.O. Da Silva, J. Tomasi, J. Chem. Phys. 12 (2000) (in press). [62] J.R. Cheeseman, M.J. Frisch, F.J. Devlin, P.J. Stephens, Chem. Phys. Lett. 252 (1996) 211. [63] F.J. Devlin, P.J. Stephens, J.R. Cheeseman, M.J. Frisch, J. Phys. Chem. A101 (1997) 6322. [64] F.J. Devlin, P.J. Stephens, J.R. Cheeseman, M.J. Frisch, J. Phys. Chem. A101 (1997) 9912. [65] J.R. Cheeseman, G.W. Trucks, T.A. Keith, M.J. Frisch, J. Chem. Phys. 104 (1996) 5497. [66] J.B. Stothers, Carbon-13 NMR Spectroscopy, Academic Press, London, 1986. [67] J.B. Stothers, Annual Reports in NMR Spectroscopy, vol. 18, Academic Press, London, 1986, p. 365.