- Email: [email protected]
Ab initio molecular orbital calculations including solvent effects by generalized Born formula. Conformation of zwitterionic forms of glycine, alanine and serine in water
79
Journal of Molecular Structure (Theochem), 305 (1994) 79-87 0166-1280/94/$07.00 0 1994 - Elsevier Science B.V. All rights reserved
Ab initio molecular orbital calculations including solvent effects by generalized Born formula. Conformation of zwitterionic forms of glycine, alanine and serine in water Osamu
Kikuchi*,
Department (Received
of 30
Takayuki
Matsuoka,
Hiroyuki
Sawahata,
Ohgi Takahashi
Chemistry, University of Tsukuba, Tsukuba 30.5, Japan
July 1993; accepted 13 August 1993)
Abstract The continuum model, in which the solvation energy is expressed by the generalized Born formula and the atomic fractional charges are evaluated by Liiwdin population analysis, has been incorporated into ab initio SCF calculations. The effective atomic radius was expressed as a function of the atomic charge and optimized in the variational procedure to obtain the Fock matrix elements. The method has been applied to conformational analysis of the zwitterionic forms of glycine, alanine and serine in water. The solvent effects on the molecular and electronic structures of these species were also examined.
Introduction Recently much attention has been paid to MO calculations including solvent effects. Among the many models which take into account solvent effects, the continuum mode1 using the generalized Born formula is a representative one [l-7]. In this mode1 the solvent is treated as a continuum with dielectric constant e, and the solute molecule is an assembly of spheres with fractional charges. The salvation energy is expressed by the generalized Born formula [S-lo],
Es“’= - ;
(1
-
;) x x A
QAQBrAB
B
(2) author.
SSDI 0166-1280(93)03460-O
1 rAB
(3)
=
Rig +
[crA
+
rB)/212
where RAB is the distance between atoms A and B. The energy of a molecule in solution, ES, is the sum of the energy of the molecule in its isolated state, E”, and the salvation energy, Es“‘. ES = E” + Es“’
where A and B are atoms in the molecule, and QA and Qn are fractional charges. FAA is a reciprocal of the atomic radius, rA.
* Corresponding
IAa is the two-center repulsion energy between unit charges on atoms A and B, and was evaluated by using an Ohno-Klopman type formula in our previous papers [3,4],
(4)
This mode1 includes only the electrostatic polarization factor and may be too simple to describe the complex solvent effects on molecular properties. This method, however, has the following advantages. (1) The mode1 is simple and can easily be incorporated into the SCF framework. The computation time required is very short. (2) The mode1 has no cavity in which the solute
80
0. Kikuchi et al./J. Mol. Struct. (Theochem)
molecule is embedded and can easily be applied to chemical reactions involving bond dissociation. (3) The model provides the wavefunction of the solute molecule in solution which can be used to evaluate physical properties of the molecule in solution. This continuum model has been incorporated with semiempirical MO methods and found to be very useful for discussing the molecular and electronic structures and energetics in solution [3-71. In this model, one-center solvation energy depends directly on the atomic radius, and the evaluation of atomic radii is important especially in the calculation of the chemical processes in which the atomic charges are expected to change significantly, as in the case of heterolysis in water. The hydration energies of the H+ and H- ions are significantly different from each other, and it is not appropriate to use a fixed radius for the hydrogen atom. In this paper, an ab initio SCF procedure is incorporated with the extended Born formula in which the atomic radius is expressed as a function of the atomic fractional charge. Although a number of quantum chemical calculations have been reported on the structure of amino acids in vacua [l 11, only a few have been reported on their preferred conformations in solution [IZ-141. The present ab initio method including solvent effects was applied to the conformational analysis of the zwitterionic forms of amino acids, glycine, t_-alanine, and L-serine. The solvent effect on the molecular and electronic structures of these species in water was examined.
energy depends directly on this definition. A most popular definition of electron population is the Mulliken gross population [ 151,
Theory A closed-shell molecule is considered for E” in Eq. (4). In order to evaluate the solvation energy, Eso’, we must define the electron population on atom, PA, and calculate the atomic charge, QA=~A--PA
(5)
where zA is the charge on the nucleus. The detinition of PA is very important since the solvation
pA = c cPs),, PEA
305 (1994) 79-87
(6)
where P is the bond-population matrix and S is the overlap-integral matrix. If the electron population is calculated by Eq. (6) and the variational theorem is applied to ES in Eq. (4), Fock matrix elements including solvent effect are obtained Fpy=cy+F$
(9) F$ is the Fock matrix element for an isolated molecule, and F$ describes the contribution of the solvent. Q, is the charge of the atom which the h basis set belongs to. The second and third terms in Eq. (9) describe the charge dependence of I’, If we use the fixed values for atomic radii, these terms vanish and the Fock matrix elements which correspond to those for the semiempirical MO method including the Born’s solvation energy [3,6] are obtained. One unfavorable point of the Mulliken population is that, when an electron in a molecular orbital is shared among atomic orbitals in the molecule, a negative amount of electrons may be assigned to a specific atomic orbital, since the overlap electron population, which can be negative, is apportioned equally between two atomic orbitals. This causes a serious problem in the present model where the atomic charge is involved in the energy functional as in Eq. (4). The Liiwdin population [16] is another definition of the electron population in an atom of a molecule
81
0. Kikuchi et al./J. Mol. Struct. (Theochem) 305 (1994) 79-87
and has been pointed out to be better than the Mulliken population [17]. It was also found in this study that if the solvent effect was estimated by the generalized Born formula, the Lowdin population gave more reasonable results than the Mulliken population. We then used the Lowdin population for PA in the present ab initio SCF calculation.
(10)
PA = ~(sl'*Psl'*),, PEA
GLYCINE
SERINE( 1)
In this case, the contribution of solvent to the Fock matrix elements becomes,
(s”2),x(s1’2)“x
(11)
The two-center interaction between atomic charges can be divided into different types of electrostatic interactions, the nucleus-nucleus, and electron-electron; these nucleus-electron, terms can be evaluated at the ab initio MO level of calculations. However, the solvation energy was
ALANINE
SERINE( 2)
Fig. 1. Definition of two rotational angles, q5and 0, for the zwitterionic forms of glycine, alanine and serine. 4 indicates the rotation of the NH, fragment, while 0 indicates the rotation of the CO2 fragment. The conformations shown are (q$@) = (0,O).
82
0. Kikuchi et al./J. Mol. Struct.
not appropriately evaluated by ab initio type calculations of these terms, since the one-center term is evaluated empirically by the Born formula, Eq. (2). In the present study, the solvation model was incorporated into ab initio calculations, where the solvation terms are evaluated semiempirically
charge dependence of the atomic radius.
by Eq. (3). The second and third terms in Eq. (9) vanish if the atomic radius is fixed in the calculations. As was pointed out, however, the charge dependence of the atomic radius is very important in calculations of ionic species and the atomic radius should be determined as a function of atomic charge in each SCF iteration. There have been some reports [3,6] in which the atomic radii in the Born formula are varied with the atomic charge. However, in these studies, Fock matrix elements have not been derived by correct variational procedures; the present study gives the first correct variational results. We adopted the following formula with two parameters, o and ,D, for the
(Theochem)
305 (1994) 79-87
rA = ae-4Q~
(12)
In order to examine how the atomic radius depends on the atomic charge, the radius of the atomic sphere in which 95% of electrons are included was calculated for a neutral atom, a mono-cation and a mono-anion for the second-row and third-row atoms by using Slater atomic orbitals. Equation (12) was found to represent appropriately the charge dependence of such theoretical radii; p = 0.074 and 0.075 were obtained for F and Cl atoms, respectively. The o parameters for F and Cl were determined empirically from the hydration energies of the F- and Cl- ions. The a and ,D parameters for other atoms except hydrogen were determined from the theoretical atomic radii and those of F and Cl atoms. The parameters for hydrogen atoms were determined from the experimental hydration energies of H+ and H-. These are given elsewhere [ 181.
lS0-
0
30
60
90
120
0
30
60
90
120
4
Fig. 2. STO-3G potential energy map of glycine zwitterion in water (E = 79). The values cited are in kcal mol-’ which are relative to the most stable conformation.
Fig. 3. 6-31G potential energy map of glycine zwitterion in water (C = 79). The values cited are in kcal mol-’ which are relative to the most stable conformation.
0. Kikuchi et al./J. Mol. Struct.
(Theochem)
305 (1994)
83
79-87
Calculations
The amino acids considered in this study are shown in Fig. 1. Experimental [ 19,201and theoretical [l l-141 studies have indicated that amino acids exist in zwitterionic form in water and in the neutral form in the gas phase. In the present study, the zwitterionic forms of glycine, L-alanine, and L-serine were calculated in water (E = 79). The conformation is described by two rotational angles, 4 and B, shown in Fig. 1. The (q$@)potential energy maps were created with the STO-3G and 6-31G basis sets by changing the q!~and 8 angles in steps of lo”, the other molecular parameters being fixed at the standard values reported [ 12,14,21]: bond lengths C-C = 1.53 A, C-O = 1.25& C-N = 1.48& C-H = 1.09& N-H = 1.03& bond angles OCO = 125”, CCN = 109.47”, CCH = 109.47”, CNH = 109.47”, and the standard parameters were used for the CH3 and CH*OH fragments in alanine and serine. Two
conformations were considered for the CH20H group in serine. In serine(l), CC”C0 = 180” and C”COH = 180”, while in serine(2), CCC0 = 60” and C*COH = 0”; these were considered as the of the serine conformations representative zwitterion [21]. Thus the NH3 and CO2 groups have the C, and C, local symmetry axes, respectively, and the (4, 0) contour map shows 120” and 180” periodicity for 4 and 8, respectively. The molecular structures of the (O,O), (60,0), (60,90), and (0,90) conformations of glycine in water were optimized with the 6-31G basis set under the C, symmetry restriction.
Results and discussion
The STO-3G (q!@) potential energy map of the glycine zwitterion (GZW) in water (6 = 79) is shown in Fig. 2. The lowest-energy conformation is the (0,O) conformation. The (60,O) conformation
1802
0
30
60
90
120
4
Fig. 4. 6-31G potential energy map of alanine zwitterion in water (e = 79). The values cited are in kcal mol-’ which are relative to the most stable conformation.
Fig. 5.6-3 IG potential energy map of serine(1) zwitterion in water (e = 79). The values cited are in kcal mol-’ which are relative to the most stable conformation.
84
0. Kikuchi et al./J. Mol. Struct. (Theochem)
Table 1 6-31G optimized structures of four conformations Geometrical parametersa
(090)
305 (1994) 79-87
of the glycine zwitterion in water
(60,O)
(60,90)
(0,90)
Experimentalb
Bond length (A) P--C
c-oc c-o’ c-o* CO-N C-H N-HC N-H’
1.515 1.269 1.270 1.473 1.082 1.024 -
1.513 1.269 1.271 _ 1.469 1.083 _ 1.031
1.520 1.271 1.481 1.079 _ 1.036
1.520 _ 1.270 1.488 1.078 1.020 _
N-Hg
1.028
1.022
1.023
1.031
1.526 1.250 1.251 1.476
Bond angle (deg) NC”C
CQCOC C”C0’
112.9 117.8 115.9
113.0 117.6 116.1
109.9 _ _
110.7 _ _
C”COf CC”H CQNHC C’NH’
108.4 110.4 -
108.8 _ 111.1
116.9 109.8 _ 110.3
117.0 109.3 113.4 _
C”NHK
112.2
111.8
111.8
111.3
122.1 _
121.7
_
_
_ 58.9
29.1 60.7
29.7 120.2
111.9 117.5 117.1
Dihedral angle (deg)
OCCC”H O*CC”H CC*NHg
119.6
a 0’ and 0’ are the cis and trans 0 atoms with respect to the NH3 group, respectively. 0* is the 0 atom in the perpendicular (0 = 90”) conformation. Superscripts c, t, and g of the H atom in the NH3 group indicate cis, trans, and gauche with respect to the C”-C bond. b The NCCOz skeleton in the crystalline state [26]; glycine zwitterion has the conformation which is very close to the (60,O) conformation.
has the energy 3.5 kcalmol-’ relative to the global minimum, and the (60,90) conformation has the energy 7.5 kcalmol-*. The (0,90) conformation corresponds to the energy maximum point in the (4, 0) map. The 6-31G potential energy map of GZW in water is shown in Fig. 3. The global minimum is found at (60,O) and the (0,O) conformation has an energy slightly higher than the (60,O) conformation. The shapes of the potential energy surface near (60,90) and (0,90) are similar to those in the STO-3G map. The 6-3 1G (q+,0) maps for L-alanine and L-serine are shown in Figs. 4 and 5, respectively. The overall
shapes of the maps resemble the 6-31G map of GZW. The global minima in these maps, (58, -7) for alanine and (52, 5) for serine, are located near the (60,O) conformation, and the energies for the (0,O) conformation are close to the global minima. These maps indicate that 8 = 0 is the most stable conformation for all 4 values. There are only a few theoretical studies on the conformational stability of GZW in water [ 12- 141. Bonaccorsi et al. [12] calculated the (c#J,I~) map of GZW in water by the continuum model with a cavity on the solute molecule that was developed for the ab initio calculations of physicochemical properties of molecules in solution. Their
0. Kikuchi et al/J.
Mol. Struct. (Theochem)
305 (1994)
85
79-87
Table 2 Calculated relative energies (keal mol-‘) of representative
confo~ations
of the zwitterions
in water
Zwitterion
Method
(W)
(60,O)
W&90)
(0,90)
Glycine
6-31G 6-31Ga STO-3G 4-31G STO-3G b _c 6-31G STO-3G 6-31G STO-3G 6-31G
0.7 0.7 0.0 0.7 0.0 6.6 0.0 0.7 0.0 1.3 0.0 0.6
0.0 0.0 3.2 0.0 2.0 0.0 2.4 0.0 3.0 0.0 2.7 0.0
5.7 6.2 7.2 5.4 4.4
9.0 9.2 10.6
Alanine Serine( 1) Serine(2)
0.9 5.9 1.2 4.8 6.6 7.1
12 12 14 24, 25
9.1 10.6 8.6 IO.4 10.3
a Structure was optimized for each conformation. b Estimated from minimal basis set calculation for the isolated state and the solvation energy estimated separately. ’ Relative energies for two conformations in the gas phase were calculated with the 6-31G* basis set by about 5 kcalmol-‘, the solvation energieswere calculated by Monte Carlo simulations.
STO-3G calculation predicted the (0,O) conformation as the most stable conformation in water, while the 4-31G calculation predicted the (60,O) The 4-31G calculations also conformation. suggest that the energy difference between the two conformations, ~O,O~and (60,0), is very small. These results agree well with the present (&?) maps in Figs. 2 and 3. Our previous MNDO continuum model calculations predicted the (60,O) conformation as the most stable conformation of GZW in water [13]. Deprick-Cote et al. [14] estimated the relative energies of three conformations of GZW in water from the energy in the gas phase calculated by a minimal basis set and the solvation energy calculated separately by a continuum model. The (60,O) conformation is predicted as the most stable conformation in water. Their results also predicted that the (0,90) conformation is very stable and higher than the (60,O) global minimum by only 0.9 kcal mol-‘. This is against the present results and those of Bonaccorsi et al. [12] in which the (0,90) conformation corresponds to the maximum point in the (&9) potential energy surface. The (0,90) conformation has been considered in the ab initio calculation of GZW-6HzO and GZW12H20 complexes [22]. However, the conforma-
Ref.
and
tional dependence of GZW stability in water was not studied there. The first MO study of a GZWHz0 system was performed by Imamura et al. [23] in which only the (60,O) conformation was ~l~ulated. Monte Carlo (MC) simulation has been apphed to GZW in water by Carozzo et al. and the solvation structure was discussed [21]. Alagona and co-workers applied their MC simulation to GZW in water and estimated the solvation energies for (0,O) and (60,O) conformations [24,25]. The solvation energy of the (60,O) conformation was calculated to be larger than that for the (0,O) conformation by 2.6 kcalmol-‘. Since this value is small, they concluded that the lowest-energy confo~ation of GZW in water is the (0,O) confo~ation, Although more theoretical efforts may be required to solve the energetic problems of GZW conformers in water, the present results suggest that the (60,O) and (0,O) conformations have almost the same stability in water. The 6, rotation requires a small activaiton energy, while the 8 rotation has an appreciable activation energy; the 0 = 0 conformation is stable. In the calculations of amino acids in water, the standard geometrical parameters have been
86
0. Kikuchi et al./J. Mol. Struct. (Theochem)
Table 3 Solvation energies (kcalmol-‘) for four conformations glycine, alanine and serine zwitterionsa
of
Zwitterion
(070)
(60,O)
(6090)
(0,90)
Glycine Alanine Serine( 1) Serine(2)
-86.4 -85.1 -88.0 -93.8
-89.3 -87.7 -90.5 -96.3
-93.1 -91.0 -93.1 -102.3
-94.4 -92.5 -95.0 -103.8
a The difference between energies in vacua (E = 1) and in solution (E = 79) which were calculated with 6-31G basis set using the standard geometrical parameters.
employed and geometry optimization by quantum chemical calculations including solvent effects have not been reported so far. It is thus interesting to present the calculated molecular structures in water. The four representative conformations of GZW in water (e = 79) were optimized with the 6-3 1G basis set; the results are listed in Table 1. The optimized molecular parameters appear to be reasonable in comparison with the experimental values for GZW in the crystalline state [26]. The relative energies for these optimized structures are listed in Table 2, and have a trend which is the same as that observed for the energies estimated using standard geometry. The solvation energy is defined here as the difference between energies in the gas phase and solution. The solvation energy depends thus on the structure adopted for the calculations in the gas phase and in solution. The molecular structure of GZW in vacua is difficult to determine by optimization, since the neutral form is more stable in the gas phase and geometry optimization of GZW gives the neutral form for some conformations of GZW [27]. In this study, the standard geometry was used and the energies for (O,O), (60,0), (60,90), and (0,90) conformations were calculated for e = 1 and E = 79, and the solvation energy was estimated as the energy difference between them; they are listed in Table 3. Two characteristics are observed in the solvation energies for different conformations: (1) the twisted forms (4 = 90) are stabilized by the solvent more than the planar (0 = 0) forms; and (2) the
305 (1994) 79-87
solvation energy of the (60,O) conformation is larger than the (0,O) conformation by 2.5 2.9 kcal mol-' . These agree well with the solvation energies for GZW estimated by Bonaccorsi et al. [12] and those obtained by MC simulations [24,25]. The hydration energy for GZW has been estimated to be -5O--99 kcalmol-’ [12,24,27,28]. Absolute values of the solvation energies calculated by the present model (Table 3) are in this range and the relative values are expected to be reasonable. References 1 0. Tapia, in R. Daudel, A. Pullman, L. Salem and A. Veillard (Eds.), Quantum Theory of Chemical Reactions, Vol. 2, Reidel, Dordrecht, 1981, p. 25. 2 R. Constanciel and R. Contreras, Theor. Chim. Acta, 65 (1984) 1. 3 T. Kozaki, K. Morihashi and 0. Kikuchi, J. Mol. Struct. (Theochem), 168 (1988) 265. 4 T. Kozaki, K. Morihashi and 0. Kikuchi, J. Am. Chem. Sot., 111 (1989) 1547. 5 S.C. Tucker and D.G. Truhlar, Chem. Phys. Lett., 157 (1989) 164. 6 C.J. Cramer and D.G. Truhlar, J. Am. Chem. Sot., 113 (1991) 8305, 8552. I C.J. Cramer and D.G. Truhlar, J. Comput. Chem., 13 (1992) 1089. 8 M. Born, Z. Phys., 45 (1920) 1. 9 G.J. Hoijtink, E. de Boer, P.H. van der Meij and W.P. Weijland, Reel. Trav. Chim. Pays-Bas, 75 (1956) 487. 10 I. Jano, C.R. Acad. Sci., 261 (1965) 103. 11 C.-H. Hu, M. Shen and H.F. Schaefer III, J. Am. Chem. Sot., 115 (1993) 2923, and references cited therein. A.G. Csaszar, J. Am. Chem. Sot., 114 (1992) 9568. 12 R. Bonaccorsi, P. Palla and J. Tomasi, J. Am. Chem. Sot., 106 (1984) 1945. 13 0. Kikuchi, T. Natsui and T. Kozaki, J. Mol. Struct. (Theochem), 207 (1990) 103. J. Langlet, J. Caillet, J. Berges, 14 B. Deprick-Cote, E. Kassab and R. Constanciel, Theor. Chim. Acta, 82 (1992) 435. 15 R.S. Mulliken, J. Chem. Phys., 23 (1955) 1833. 16 P.-O. Lowdin, J. Chem. Phys., 18 (1950) 365. 17 E.W. Stout, Jr., and P. Politzer, Theor. Chim. Acta, 12 (1968) 379. 18 0. Kikuchi and 0. Takahashi, JCPE Newslett., 5 (1993) 24. 19 R.D. Brown, P.D. Godfrey, J.W.V. Storey and M.-P. Bassez, J. Chem. Sot., Chem. Commun., (1978) 541.
79-87
87
J.S. Gaffney, R.C. Pierce and L. Friedman, J. Am. Chem. Sot., 99 (1977) 4293. L. Carozzo, G. Corongiu, C. Petrongolo and E. Clementi, J. Chem. Phys., 68 (1978) 787. W. Fiirner, P. Otto, J. Bernhardt and J.J. Ladik, Theor. Chim. Acta, 60 (1981) 269. A. Imamura, H. Fujita and C. Nagata, Bull. Chem. Sot. Jpn., 42 (1969) 3118. G. Alagona, C. Ghio and P.A. Kollman, J. Mol. Struct. (Theochem), 166 (1988) 385.
25 G. Alagona and C. Ghio, J. Mol. Liq., 47 (1990) 139. 26 P.-G. JGnsson and A. Kvick, Acta Crystallogr., Sect. B, 28 (1972) 1822. 27 Y.-C. Tse, M.D. Newton, S. Vishveshwara and J.A. Pople, J. Am. Chem. Sot., 100 (1978) 4329. 28 M. Mezei, P.K. Mehrotra and D.L. Beveridge, J. Biomol. Struct. Dyn., 2 (1984) 1.
0. Kikuchi et al./J. Mol. Struct.
20
21 22 23 24
(Theochem)
305 (1994)
Journal of Molecular Structure (Theochem), 305 (1994) 79-87 0166-1280/94/$07.00 0 1994 - Elsevier Science B.V. All rights reserved
Ab initio molecular orbital calculations including solvent effects by generalized Born formula. Conformation of zwitterionic forms of glycine, alanine and serine in water Osamu
Kikuchi*,
Department (Received
of 30
Takayuki
Matsuoka,
Hiroyuki
Sawahata,
Ohgi Takahashi
Chemistry, University of Tsukuba, Tsukuba 30.5, Japan
July 1993; accepted 13 August 1993)
Abstract The continuum model, in which the solvation energy is expressed by the generalized Born formula and the atomic fractional charges are evaluated by Liiwdin population analysis, has been incorporated into ab initio SCF calculations. The effective atomic radius was expressed as a function of the atomic charge and optimized in the variational procedure to obtain the Fock matrix elements. The method has been applied to conformational analysis of the zwitterionic forms of glycine, alanine and serine in water. The solvent effects on the molecular and electronic structures of these species were also examined.
Introduction Recently much attention has been paid to MO calculations including solvent effects. Among the many models which take into account solvent effects, the continuum mode1 using the generalized Born formula is a representative one [l-7]. In this mode1 the solvent is treated as a continuum with dielectric constant e, and the solute molecule is an assembly of spheres with fractional charges. The salvation energy is expressed by the generalized Born formula [S-lo],
Es“’= - ;
(1
-
;) x x A
QAQBrAB
B
(2) author.
SSDI 0166-1280(93)03460-O
1 rAB
(3)
=
Rig +
[crA
+
rB)/212
where RAB is the distance between atoms A and B. The energy of a molecule in solution, ES, is the sum of the energy of the molecule in its isolated state, E”, and the salvation energy, Es“‘. ES = E” + Es“’
where A and B are atoms in the molecule, and QA and Qn are fractional charges. FAA is a reciprocal of the atomic radius, rA.
* Corresponding
IAa is the two-center repulsion energy between unit charges on atoms A and B, and was evaluated by using an Ohno-Klopman type formula in our previous papers [3,4],
(4)
This mode1 includes only the electrostatic polarization factor and may be too simple to describe the complex solvent effects on molecular properties. This method, however, has the following advantages. (1) The mode1 is simple and can easily be incorporated into the SCF framework. The computation time required is very short. (2) The mode1 has no cavity in which the solute
80
0. Kikuchi et al./J. Mol. Struct. (Theochem)
molecule is embedded and can easily be applied to chemical reactions involving bond dissociation. (3) The model provides the wavefunction of the solute molecule in solution which can be used to evaluate physical properties of the molecule in solution. This continuum model has been incorporated with semiempirical MO methods and found to be very useful for discussing the molecular and electronic structures and energetics in solution [3-71. In this model, one-center solvation energy depends directly on the atomic radius, and the evaluation of atomic radii is important especially in the calculation of the chemical processes in which the atomic charges are expected to change significantly, as in the case of heterolysis in water. The hydration energies of the H+ and H- ions are significantly different from each other, and it is not appropriate to use a fixed radius for the hydrogen atom. In this paper, an ab initio SCF procedure is incorporated with the extended Born formula in which the atomic radius is expressed as a function of the atomic fractional charge. Although a number of quantum chemical calculations have been reported on the structure of amino acids in vacua [l 11, only a few have been reported on their preferred conformations in solution [IZ-141. The present ab initio method including solvent effects was applied to the conformational analysis of the zwitterionic forms of amino acids, glycine, t_-alanine, and L-serine. The solvent effect on the molecular and electronic structures of these species in water was examined.
energy depends directly on this definition. A most popular definition of electron population is the Mulliken gross population [ 151,
Theory A closed-shell molecule is considered for E” in Eq. (4). In order to evaluate the solvation energy, Eso’, we must define the electron population on atom, PA, and calculate the atomic charge, QA=~A--PA
(5)
where zA is the charge on the nucleus. The detinition of PA is very important since the solvation
pA = c cPs),, PEA
305 (1994) 79-87
(6)
where P is the bond-population matrix and S is the overlap-integral matrix. If the electron population is calculated by Eq. (6) and the variational theorem is applied to ES in Eq. (4), Fock matrix elements including solvent effect are obtained Fpy=cy+F$
(9) F$ is the Fock matrix element for an isolated molecule, and F$ describes the contribution of the solvent. Q, is the charge of the atom which the h basis set belongs to. The second and third terms in Eq. (9) describe the charge dependence of I’, If we use the fixed values for atomic radii, these terms vanish and the Fock matrix elements which correspond to those for the semiempirical MO method including the Born’s solvation energy [3,6] are obtained. One unfavorable point of the Mulliken population is that, when an electron in a molecular orbital is shared among atomic orbitals in the molecule, a negative amount of electrons may be assigned to a specific atomic orbital, since the overlap electron population, which can be negative, is apportioned equally between two atomic orbitals. This causes a serious problem in the present model where the atomic charge is involved in the energy functional as in Eq. (4). The Liiwdin population [16] is another definition of the electron population in an atom of a molecule
81
0. Kikuchi et al./J. Mol. Struct. (Theochem) 305 (1994) 79-87
and has been pointed out to be better than the Mulliken population [17]. It was also found in this study that if the solvent effect was estimated by the generalized Born formula, the Lowdin population gave more reasonable results than the Mulliken population. We then used the Lowdin population for PA in the present ab initio SCF calculation.
(10)
PA = ~(sl'*Psl'*),, PEA
GLYCINE
SERINE( 1)
In this case, the contribution of solvent to the Fock matrix elements becomes,
(s”2),x(s1’2)“x
(11)
The two-center interaction between atomic charges can be divided into different types of electrostatic interactions, the nucleus-nucleus, and electron-electron; these nucleus-electron, terms can be evaluated at the ab initio MO level of calculations. However, the solvation energy was
ALANINE
SERINE( 2)
Fig. 1. Definition of two rotational angles, q5and 0, for the zwitterionic forms of glycine, alanine and serine. 4 indicates the rotation of the NH, fragment, while 0 indicates the rotation of the CO2 fragment. The conformations shown are (q$@) = (0,O).
82
0. Kikuchi et al./J. Mol. Struct.
not appropriately evaluated by ab initio type calculations of these terms, since the one-center term is evaluated empirically by the Born formula, Eq. (2). In the present study, the solvation model was incorporated into ab initio calculations, where the solvation terms are evaluated semiempirically
charge dependence of the atomic radius.
by Eq. (3). The second and third terms in Eq. (9) vanish if the atomic radius is fixed in the calculations. As was pointed out, however, the charge dependence of the atomic radius is very important in calculations of ionic species and the atomic radius should be determined as a function of atomic charge in each SCF iteration. There have been some reports [3,6] in which the atomic radii in the Born formula are varied with the atomic charge. However, in these studies, Fock matrix elements have not been derived by correct variational procedures; the present study gives the first correct variational results. We adopted the following formula with two parameters, o and ,D, for the
(Theochem)
305 (1994) 79-87
rA = ae-4Q~
(12)
In order to examine how the atomic radius depends on the atomic charge, the radius of the atomic sphere in which 95% of electrons are included was calculated for a neutral atom, a mono-cation and a mono-anion for the second-row and third-row atoms by using Slater atomic orbitals. Equation (12) was found to represent appropriately the charge dependence of such theoretical radii; p = 0.074 and 0.075 were obtained for F and Cl atoms, respectively. The o parameters for F and Cl were determined empirically from the hydration energies of the F- and Cl- ions. The a and ,D parameters for other atoms except hydrogen were determined from the theoretical atomic radii and those of F and Cl atoms. The parameters for hydrogen atoms were determined from the experimental hydration energies of H+ and H-. These are given elsewhere [ 181.
lS0-
0
30
60
90
120
0
30
60
90
120
4
Fig. 2. STO-3G potential energy map of glycine zwitterion in water (E = 79). The values cited are in kcal mol-’ which are relative to the most stable conformation.
Fig. 3. 6-31G potential energy map of glycine zwitterion in water (C = 79). The values cited are in kcal mol-’ which are relative to the most stable conformation.
0. Kikuchi et al./J. Mol. Struct.
(Theochem)
305 (1994)
83
79-87
Calculations
The amino acids considered in this study are shown in Fig. 1. Experimental [ 19,201and theoretical [l l-141 studies have indicated that amino acids exist in zwitterionic form in water and in the neutral form in the gas phase. In the present study, the zwitterionic forms of glycine, L-alanine, and L-serine were calculated in water (E = 79). The conformation is described by two rotational angles, 4 and B, shown in Fig. 1. The (q$@)potential energy maps were created with the STO-3G and 6-31G basis sets by changing the q!~and 8 angles in steps of lo”, the other molecular parameters being fixed at the standard values reported [ 12,14,21]: bond lengths C-C = 1.53 A, C-O = 1.25& C-N = 1.48& C-H = 1.09& N-H = 1.03& bond angles OCO = 125”, CCN = 109.47”, CCH = 109.47”, CNH = 109.47”, and the standard parameters were used for the CH3 and CH*OH fragments in alanine and serine. Two
conformations were considered for the CH20H group in serine. In serine(l), CC”C0 = 180” and C”COH = 180”, while in serine(2), CCC0 = 60” and C*COH = 0”; these were considered as the of the serine conformations representative zwitterion [21]. Thus the NH3 and CO2 groups have the C, and C, local symmetry axes, respectively, and the (4, 0) contour map shows 120” and 180” periodicity for 4 and 8, respectively. The molecular structures of the (O,O), (60,0), (60,90), and (0,90) conformations of glycine in water were optimized with the 6-31G basis set under the C, symmetry restriction.
Results and discussion
The STO-3G (q!@) potential energy map of the glycine zwitterion (GZW) in water (6 = 79) is shown in Fig. 2. The lowest-energy conformation is the (0,O) conformation. The (60,O) conformation
1802
0
30
60
90
120
4
Fig. 4. 6-31G potential energy map of alanine zwitterion in water (e = 79). The values cited are in kcal mol-’ which are relative to the most stable conformation.
Fig. 5.6-3 IG potential energy map of serine(1) zwitterion in water (e = 79). The values cited are in kcal mol-’ which are relative to the most stable conformation.
84
0. Kikuchi et al./J. Mol. Struct. (Theochem)
Table 1 6-31G optimized structures of four conformations Geometrical parametersa
(090)
305 (1994) 79-87
of the glycine zwitterion in water
(60,O)
(60,90)
(0,90)
Experimentalb
Bond length (A) P--C
c-oc c-o’ c-o* CO-N C-H N-HC N-H’
1.515 1.269 1.270 1.473 1.082 1.024 -
1.513 1.269 1.271 _ 1.469 1.083 _ 1.031
1.520 1.271 1.481 1.079 _ 1.036
1.520 _ 1.270 1.488 1.078 1.020 _
N-Hg
1.028
1.022
1.023
1.031
1.526 1.250 1.251 1.476
Bond angle (deg) NC”C
CQCOC C”C0’
112.9 117.8 115.9
113.0 117.6 116.1
109.9 _ _
110.7 _ _
C”COf CC”H CQNHC C’NH’
108.4 110.4 -
108.8 _ 111.1
116.9 109.8 _ 110.3
117.0 109.3 113.4 _
C”NHK
112.2
111.8
111.8
111.3
122.1 _
121.7
_
_
_ 58.9
29.1 60.7
29.7 120.2
111.9 117.5 117.1
Dihedral angle (deg)
OCCC”H O*CC”H CC*NHg
119.6
a 0’ and 0’ are the cis and trans 0 atoms with respect to the NH3 group, respectively. 0* is the 0 atom in the perpendicular (0 = 90”) conformation. Superscripts c, t, and g of the H atom in the NH3 group indicate cis, trans, and gauche with respect to the C”-C bond. b The NCCOz skeleton in the crystalline state [26]; glycine zwitterion has the conformation which is very close to the (60,O) conformation.
has the energy 3.5 kcalmol-’ relative to the global minimum, and the (60,90) conformation has the energy 7.5 kcalmol-*. The (0,90) conformation corresponds to the energy maximum point in the (4, 0) map. The 6-31G potential energy map of GZW in water is shown in Fig. 3. The global minimum is found at (60,O) and the (0,O) conformation has an energy slightly higher than the (60,O) conformation. The shapes of the potential energy surface near (60,90) and (0,90) are similar to those in the STO-3G map. The 6-3 1G (q+,0) maps for L-alanine and L-serine are shown in Figs. 4 and 5, respectively. The overall
shapes of the maps resemble the 6-31G map of GZW. The global minima in these maps, (58, -7) for alanine and (52, 5) for serine, are located near the (60,O) conformation, and the energies for the (0,O) conformation are close to the global minima. These maps indicate that 8 = 0 is the most stable conformation for all 4 values. There are only a few theoretical studies on the conformational stability of GZW in water [ 12- 141. Bonaccorsi et al. [12] calculated the (c#J,I~) map of GZW in water by the continuum model with a cavity on the solute molecule that was developed for the ab initio calculations of physicochemical properties of molecules in solution. Their
0. Kikuchi et al/J.
Mol. Struct. (Theochem)
305 (1994)
85
79-87
Table 2 Calculated relative energies (keal mol-‘) of representative
confo~ations
of the zwitterions
in water
Zwitterion
Method
(W)
(60,O)
W&90)
(0,90)
Glycine
6-31G 6-31Ga STO-3G 4-31G STO-3G b _c 6-31G STO-3G 6-31G STO-3G 6-31G
0.7 0.7 0.0 0.7 0.0 6.6 0.0 0.7 0.0 1.3 0.0 0.6
0.0 0.0 3.2 0.0 2.0 0.0 2.4 0.0 3.0 0.0 2.7 0.0
5.7 6.2 7.2 5.4 4.4
9.0 9.2 10.6
Alanine Serine( 1) Serine(2)
0.9 5.9 1.2 4.8 6.6 7.1
12 12 14 24, 25
9.1 10.6 8.6 IO.4 10.3
a Structure was optimized for each conformation. b Estimated from minimal basis set calculation for the isolated state and the solvation energy estimated separately. ’ Relative energies for two conformations in the gas phase were calculated with the 6-31G* basis set by about 5 kcalmol-‘, the solvation energieswere calculated by Monte Carlo simulations.
STO-3G calculation predicted the (0,O) conformation as the most stable conformation in water, while the 4-31G calculation predicted the (60,O) The 4-31G calculations also conformation. suggest that the energy difference between the two conformations, ~O,O~and (60,0), is very small. These results agree well with the present (&?) maps in Figs. 2 and 3. Our previous MNDO continuum model calculations predicted the (60,O) conformation as the most stable conformation of GZW in water [13]. Deprick-Cote et al. [14] estimated the relative energies of three conformations of GZW in water from the energy in the gas phase calculated by a minimal basis set and the solvation energy calculated separately by a continuum model. The (60,O) conformation is predicted as the most stable conformation in water. Their results also predicted that the (0,90) conformation is very stable and higher than the (60,O) global minimum by only 0.9 kcal mol-‘. This is against the present results and those of Bonaccorsi et al. [12] in which the (0,90) conformation corresponds to the maximum point in the (&9) potential energy surface. The (0,90) conformation has been considered in the ab initio calculation of GZW-6HzO and GZW12H20 complexes [22]. However, the conforma-
Ref.
and
tional dependence of GZW stability in water was not studied there. The first MO study of a GZWHz0 system was performed by Imamura et al. [23] in which only the (60,O) conformation was ~l~ulated. Monte Carlo (MC) simulation has been apphed to GZW in water by Carozzo et al. and the solvation structure was discussed [21]. Alagona and co-workers applied their MC simulation to GZW in water and estimated the solvation energies for (0,O) and (60,O) conformations [24,25]. The solvation energy of the (60,O) conformation was calculated to be larger than that for the (0,O) conformation by 2.6 kcalmol-‘. Since this value is small, they concluded that the lowest-energy confo~ation of GZW in water is the (0,O) confo~ation, Although more theoretical efforts may be required to solve the energetic problems of GZW conformers in water, the present results suggest that the (60,O) and (0,O) conformations have almost the same stability in water. The 6, rotation requires a small activaiton energy, while the 8 rotation has an appreciable activation energy; the 0 = 0 conformation is stable. In the calculations of amino acids in water, the standard geometrical parameters have been
86
0. Kikuchi et al./J. Mol. Struct. (Theochem)
Table 3 Solvation energies (kcalmol-‘) for four conformations glycine, alanine and serine zwitterionsa
of
Zwitterion
(070)
(60,O)
(6090)
(0,90)
Glycine Alanine Serine( 1) Serine(2)
-86.4 -85.1 -88.0 -93.8
-89.3 -87.7 -90.5 -96.3
-93.1 -91.0 -93.1 -102.3
-94.4 -92.5 -95.0 -103.8
a The difference between energies in vacua (E = 1) and in solution (E = 79) which were calculated with 6-31G basis set using the standard geometrical parameters.
employed and geometry optimization by quantum chemical calculations including solvent effects have not been reported so far. It is thus interesting to present the calculated molecular structures in water. The four representative conformations of GZW in water (e = 79) were optimized with the 6-3 1G basis set; the results are listed in Table 1. The optimized molecular parameters appear to be reasonable in comparison with the experimental values for GZW in the crystalline state [26]. The relative energies for these optimized structures are listed in Table 2, and have a trend which is the same as that observed for the energies estimated using standard geometry. The solvation energy is defined here as the difference between energies in the gas phase and solution. The solvation energy depends thus on the structure adopted for the calculations in the gas phase and in solution. The molecular structure of GZW in vacua is difficult to determine by optimization, since the neutral form is more stable in the gas phase and geometry optimization of GZW gives the neutral form for some conformations of GZW [27]. In this study, the standard geometry was used and the energies for (O,O), (60,0), (60,90), and (0,90) conformations were calculated for e = 1 and E = 79, and the solvation energy was estimated as the energy difference between them; they are listed in Table 3. Two characteristics are observed in the solvation energies for different conformations: (1) the twisted forms (4 = 90) are stabilized by the solvent more than the planar (0 = 0) forms; and (2) the
305 (1994) 79-87
solvation energy of the (60,O) conformation is larger than the (0,O) conformation by 2.5 2.9 kcal mol-' . These agree well with the solvation energies for GZW estimated by Bonaccorsi et al. [12] and those obtained by MC simulations [24,25]. The hydration energy for GZW has been estimated to be -5O--99 kcalmol-’ [12,24,27,28]. Absolute values of the solvation energies calculated by the present model (Table 3) are in this range and the relative values are expected to be reasonable. References 1 0. Tapia, in R. Daudel, A. Pullman, L. Salem and A. Veillard (Eds.), Quantum Theory of Chemical Reactions, Vol. 2, Reidel, Dordrecht, 1981, p. 25. 2 R. Constanciel and R. Contreras, Theor. Chim. Acta, 65 (1984) 1. 3 T. Kozaki, K. Morihashi and 0. Kikuchi, J. Mol. Struct. (Theochem), 168 (1988) 265. 4 T. Kozaki, K. Morihashi and 0. Kikuchi, J. Am. Chem. Sot., 111 (1989) 1547. 5 S.C. Tucker and D.G. Truhlar, Chem. Phys. Lett., 157 (1989) 164. 6 C.J. Cramer and D.G. Truhlar, J. Am. Chem. Sot., 113 (1991) 8305, 8552. I C.J. Cramer and D.G. Truhlar, J. Comput. Chem., 13 (1992) 1089. 8 M. Born, Z. Phys., 45 (1920) 1. 9 G.J. Hoijtink, E. de Boer, P.H. van der Meij and W.P. Weijland, Reel. Trav. Chim. Pays-Bas, 75 (1956) 487. 10 I. Jano, C.R. Acad. Sci., 261 (1965) 103. 11 C.-H. Hu, M. Shen and H.F. Schaefer III, J. Am. Chem. Sot., 115 (1993) 2923, and references cited therein. A.G. Csaszar, J. Am. Chem. Sot., 114 (1992) 9568. 12 R. Bonaccorsi, P. Palla and J. Tomasi, J. Am. Chem. Sot., 106 (1984) 1945. 13 0. Kikuchi, T. Natsui and T. Kozaki, J. Mol. Struct. (Theochem), 207 (1990) 103. J. Langlet, J. Caillet, J. Berges, 14 B. Deprick-Cote, E. Kassab and R. Constanciel, Theor. Chim. Acta, 82 (1992) 435. 15 R.S. Mulliken, J. Chem. Phys., 23 (1955) 1833. 16 P.-O. Lowdin, J. Chem. Phys., 18 (1950) 365. 17 E.W. Stout, Jr., and P. Politzer, Theor. Chim. Acta, 12 (1968) 379. 18 0. Kikuchi and 0. Takahashi, JCPE Newslett., 5 (1993) 24. 19 R.D. Brown, P.D. Godfrey, J.W.V. Storey and M.-P. Bassez, J. Chem. Sot., Chem. Commun., (1978) 541.
79-87
87
J.S. Gaffney, R.C. Pierce and L. Friedman, J. Am. Chem. Sot., 99 (1977) 4293. L. Carozzo, G. Corongiu, C. Petrongolo and E. Clementi, J. Chem. Phys., 68 (1978) 787. W. Fiirner, P. Otto, J. Bernhardt and J.J. Ladik, Theor. Chim. Acta, 60 (1981) 269. A. Imamura, H. Fujita and C. Nagata, Bull. Chem. Sot. Jpn., 42 (1969) 3118. G. Alagona, C. Ghio and P.A. Kollman, J. Mol. Struct. (Theochem), 166 (1988) 385.
25 G. Alagona and C. Ghio, J. Mol. Liq., 47 (1990) 139. 26 P.-G. JGnsson and A. Kvick, Acta Crystallogr., Sect. B, 28 (1972) 1822. 27 Y.-C. Tse, M.D. Newton, S. Vishveshwara and J.A. Pople, J. Am. Chem. Sot., 100 (1978) 4329. 28 M. Mezei, P.K. Mehrotra and D.L. Beveridge, J. Biomol. Struct. Dyn., 2 (1984) 1.
0. Kikuchi et al./J. Mol. Struct.
20
21 22 23 24
(Theochem)
305 (1994)