Vibrational structure of the solvated glycine zwitterion

Chemical Physics 153 ( 199 1) 5 l-62 North-Holland

Vibrational structure of the solvated glycine zwitterion Joseph S. Alper Department of Chemistry University of Massachusetts-Boston, Boston, MA 02125, USA

Hoang Dothe and David F. Coker Department of Chemistry, Boston University, Boston, MA 02215, USA Received 20 April 1990; in final form 12 November

1990

The vibrational frequencies of glycine in aqueous solution are calculated. The effects of the water-glycine interactions on the forces and force constants of glycine, calculated using an ab initio SCF method, are treated by means of a molecular dynamics simulation of a system consisting of a glycine molecule surrounded by 253 water molecules. These interactions cannot be ignored because the isolated glycine zwitterion is unstable and is stabilized by the solvation. The fundamental frequencies of the glycine zwitterion and three of its deuterated isotopomers were calculated from the molecular dynamics forces and force constants. The resulting frequencies, after scaling with a single factor for the four isotopomers, are in good agreement with the experimental

1. Introduction

Ab initio calculations of vibrational structure and spectra generally assume that the molecule of interest is in the gas phase. The sample used for the experimental spectra, on the other hand, is often a solution of the molecule in some solvent. In order to relate the calculated to the experimental spectra, it is usually assumed that the effect of the solvent on the molecule is small. For many systems this approximation is a good one, especially in view of the inaccuracies in the ab initio calculations themselves. In fact, in a typical ab initio study, the amount of discussion of the discrepancies between theory and experiment caused by solvent effects is very small compared with the discussion of the errors caused by, for example, the limited size of the basis set used in the calculation of the molecular orbitals or the failure of the harmonic approximation. For molecules such as amino acids and glycine in particular, the effect of the solvent cannot be ignored; the structure of the molecule is completely different in the gas and liquid phases. In the gas phase glycine is a nonionic molecule, NHICHICOOH, while in water solution it is a zwitterion, NH: CH2COO-. If 0301-0104/91/S

an ab initio calculation of the vibrational structure of glycine is performed assuming that the gas phase form of the molecule is the zwitterion, several difficulties arise. First, the force on each of the glycine atoms is not zero because the geometry of the zwitterion is not the geometry at which the theoretical calculation gives the minimum energy. These nonzero forces must be included in the expression for the transformation of the force constant matrix expressed in Cartesian coordinates into the corresponding matrix expressed in internal coordinates. If the forces were zero, the transformation would depend on the B matrix alone which transforms the Cartesian coordinates into a set of internals [ 11. Pulay [ 2 ] has shown that the appropriate transformation for the case of nonzero forces is given by P=B”m-‘-

C @IB*ic’B-‘,

(1)

where, in the notation of Pulay, K and F are the harmonic force constant matrices in Cartesian and internal coordinates, respectively, @is the column vector of the forces expressed in internal coordinates, C’ is the second-order transformation matrix relating the

03.50 0 1991 - Elsevier Science Publishers B.V. (North-Holland)

52

J.S. Alper et al. / Vlbrattonal structure of the solvated glycrne zwtterlon

Cartesian and internal coordinates, poseofB”,andB”isgivenby B”=(BmB+)-‘Bm,

B-’ is the trans-

(2)

where m is any matrix for which (BmB+ ) is not singular. Second, and of crucial importance, is the difficulty that because of the existence of repulsive forces on several of the atoms of the isolated glycine zwitterion, one of the fundamental frequencies is imaginary. This imaginary frequency corresponds to the torsion normal mode involving the COO- group. It is as if this mode is trying to twist the molecule so that it will be able to form hydrogen bonds between the oxygen atoms and the hydrogens on the nitrogen atom. Clearly, in aqueous solution, the water molecules surrounding the glycine stabilize the zwitterion by reducing the magnitude of the average force on each of the glycine atoms. In this paper we study the vibrational spectrum of the glycine zwitterion in aqueous solution. We simulate the effect of the solvent by means of a molecular dynamics (MI?) calculation of a glycine molecule surrounded by 253 water molecules. The forces and force constants of an isolated glycine calculated by means of an ab initio calculation of the “gas phase” zwitterion are modified by the effects of water-glytine interactions. Once the effective forces and force constants for glycine have been obtained, the vibrational analysis proceeds in exactly the same manner as it would in a gas phase calculation. The paper is organized as follows: First we discuss the various potential energy functions used in the molecular dynamics calculation. Then we report those aspects of the molecular dynamics procedure specific to this study. After describing the determination of the effective forces and force constants, we present the vibrational analysis, discuss the results, and compare the results with various experimental spectra.

2. Potential energy functions

The total potential energy consists of the sum of four terms: the intramolecular glycine potential energy, the water-glycine potential, the intramolecular

water potential, and the water-water potential. We discuss each in turn.

intermolecular

2. I. Glycine intramolecular potential energy Our calculation of the glycine intramolecular potential surface for use in the condensed phase calculations is somewhat nonstandard. Usually such a surface is obtained by constructing an empirical force field from the experimental gas phase vibrational frequencies. This surface is then combined with an intermolecular potential function to predict solution phase properties. This procedure is not applicable here since, as mentioned above, the molecule vibrates about a stable nonionic geometry in the gas phase and about a zwitterionic form in the liquid phase. The experimental vibrational frequencies of the solution phase zwitterionic form cannot be used to determine the purely intramolecular surface because the experimental frequencies already include the influence of the intermolecular forces on the vibrational frequencies. If these experimental frequencies were used, the influence of intermolecular forces would be counted twice. The procedure we follow here involves the assumption that the solvated molecule vibrates about some zwitterionic geometry. This initially unknown geometry is approximated by an appropriately chosen zwitterionic geometry. We refer to the former geometry as the reference geometry and to the latter as the original geometry. Ab initio electronic structure calculations of the potential energy and its first and second derivatives with respect to Cartesian displacements of each of the atoms are then performed on this original geometry. This zwitterionic geometry is not at a potential minimum in the gas phase; the first derivatives are nonzero. In fact the molecule is unstable with at least one of the vibrational frequencies being imaginary. If the glycine molecule in solution does not make substantial excursions away from our original geometry, then we can use the ab initio second-order Taylor series expansion for the intramolecular potential and combine it with an intermolecular potential. This composite potential will result in a stable solution phase system, i.e., the intermolecular solvent forces will balance the ab initio intramolecular forces, and the average force on each glycine atom will be zero.

J.S. Alper et al. / Vibrational structure ofthe solvated glycine zwitterion

Our molecular dynamics calculations probe the process by which the solvent forces stabilize the original zwitterionic geometry producing the reference geometry about which the solvated glycine vibrates. The first step in obtaining the intramolecular glytine potential is the determination of the appropriate original geometry. As mentioned earlier, the optimized theoretical geometry describes the structure of the nonionic glycine so it cannot be used. The geometries of the glycine zwitterion found in the literature are based on the crystalline form [ 3 1. There is no reason to believe that the geometry in solution is the same as that in the crystal. However, it is probably not necessary to determine the original geometry of the isolated molecule with absolute precision since the atoms of the glycine molecule are free to move in the course of the molecular dynamics calculation. In view of these considerations, we chose to symmetrize the glycine crystalline geometry by using equal bond lengths for bonds of the same type, tetrahedral bond angles, and by setting the dihedral angle between the COO and NCC plane to be zero. Two possible original zwitterionic geometries must be considered: the CHz can be either staggered or eclipsed with respect to the NH3 group. Although the eclipsed form is the lower energy one in the gas phase, we found that the staggered form is more stable in aqueous solution. Fig. 1 shows the geometry of the staggered form of zwitterionic glycine obtained from the molecular dynamics calculation reported in this paper. Table 1 lists the various bond lengths, bond angles, and torsions for this geometry, for the geometry of the crystal [ 3 1, and for the original symme-

H6

01

Fig. 1. Molecular zwitterion.

H3

dynamics

geometry

of the solvated

glycine

53

trized geometry. Using this original geometry, we calculated the forces, h, and the force constants, f;i, for glycine using GAUSSIAN82 with a 6-31G basis set [ 41. The subscripts i and j label the Cartesian coordinates of the glycine atoms. Since -A and f;l are the first and second derivatives, respectively, of the intramolecular potential energy of glycine, we can use the forces and force constants to approximate the potential by means of the first three terms of a Taylor series based at the original geometry: I

(3) where the s, are the 30 Cartesian coordinates of glytine, the sp are the coordinates of glycine at the original geometry, and V. is the potential energy at the original geometry. Since V. is a constant and merely defines a zero of energy, it can be set equal to zero with no loss of generality. In using this ab initio intramolecular glycine potential for calculating the vibrational frequencies of the solvated glycine, we are assuming that the reference geometry of the solvated glycine is close to that of the staggered symmetric original geometry used to calculate the intramolecular glycine potential. Our MD calculations show that this assumption is justitied. We must also assume that the vibrational motion of the solvated molecule can be described by small harmonic oscillations about the reference glytine geometry. This assumption is reasonable since the strong hydrogen bonding network of the glycine and the surrounding solvent molecules restricts the motion about the reference geometry to small oscillations. This work also assumes that the staggered conformer of glycine is sufficiently more stable than the eclipsed so that the existence of the eclipsed conformer can be ignored, i.e., exp [ - (E, -Es) /kT] CK 1, where E, and Es are the energies of the eclipsed and staggered systems, respectively. If the eclipsed form of the zwitterion were accessible by thermal fluctuations, we would need to consider harmonic motions over a potential surface based on this conformer as well as motions over the surface derived for the staggered conformer. It would be of interest

54

J.S. Alper et al / Vibrational structure of the solvated glycme zwltterlon

Table 1 Geometry

of the glycine zwitterion Internal coordinate

Bond length (A), bond angle, or torsion symmetrized

Nl-C2 Nl-H3 NI-H4 Nl-HS C2-H6 C2-H7 C2-C8 CS-09 C8-010 C2-Nl-H3 C2-Nl-H4 C2-N I -H5 Nl-C2-H6 Nl-C2-H7 N 1-C2-C8 C2-C8-09 C2-C8-010 H3-Nl-C2-H6 H4-Nl-C2-H6 H5-N l-C2-H6 H3-N 1-C2-H7 H3-Nl-C2-C8 N 1-C2-C8-09 N 1-C2-C8-0 10 ‘) Table 7 ofref.

1.48 1.03 1.03 1.03 1.10 1.10 1.53 1.25 1.25 109.5 109.5 109.5 109.5 109.5 109.5 117.5 117.5 60.0 180.0 -60.0 -60.0 180.0 0.0 - 180.0

molecular dynamics 1.510

1.064 1.051 1.033 1.081 1.080 1.536 1.245 1.260 115.82 117.01 113.76 105.28 105.42 112.16 114.73 113.84 48.50 169.38 - 72.07 - 65.08 170.30 21.95 - 158.03

crystal a)

1.476 1.054 1.037 1.025 1.089 1.090 1.526 1.250 1.250 112.09 111.73 110.37 109.05 108.51 111.85 117.46 117.09 55.6 177.9 -63.7 -62.8 182.7 18.9 -161.8

[ 181.

to investigate a model tween conformers.

for the interconversion

2.2. Water-glycinepotential

be-

energy

The expression for the water-glycine intermolecular potential function used here was derived by Carozzo et al. [ 5 1. This potential was obtained by fitting the results of an SCF-MO calculation to the following form:

+C,qtq,lr, +D~jqiq,lr~) )

charges and potential parameters are tabulated in [ 5 1. We might note that the numerical values of terms involving C and D given in table 5 of ref. must be divided by 47~~ in accordance with the mensions of C and D specified in that table.

ref. the [ 5] di-

2.3. Intramolecular water potential The intramolecular water potential used here was derived by Reimers and Watts [ 6 ] as modified by Coker et al. [ 71 and described in Coker and Watts [ 8 1. This potential takes the form

(4)

where r,, is the distance between the ith atom of glytine and thejth atom of water, qi and q, are effective charges on the atoms of glycine and water, and A,,, B,,, C,, and D, are parameters. The sum is over the 10 glycink atoms and 3 water atoms. Values of the

where the sum is over the three local modes sj which are functions of the two O-H bond lengths and the included angle. The functions M( s ) are Morse oscillator potentials and fiZ is a coupling parameter. A

J.S. Alper et al. / Vibrational structure of the solvatedglycine zwitterion

complete discussion of this potential, with expressions for the local modes and Morse functions together with the values of the parameters characterizing the potential are given by Coker and Watts [ 8 1. 2.4. Water-water intermolecular potential For the water-water intermolecular potential we use the semiempirical RWK2 model of Reimers et al. [ 91. This potential contains terms involving both atom-atom interactions and Coulomb terms involving a negative charge located on the bisector of the HOH angle of each water molecule. These Coulomb terms involve the attraction of the point charge on one molecule with the hydrogen atoms on another and the repulsive interaction between the negative charges on two different molecules. The negative charges are included in order to reproduce the experimental dipole and quadrupole moments and to obtain good values for the dispersion coefficients. Because the potential is rather complicated we do not reproduce it here. Full details are given by Coker and Watts [ 8 1.

3. Molecular dynamics calculations 3.1. Implementation of the molecular dynamics procedure

We first performed a molecular dynamics calculation on water alone. This was done for three reasons. First, we wished to check results obtained using our MD procedure with the results obtained in the calculations of water performed by Reimers and coworkers [ 9, lo] who used potential functions similar to those used here; second, to obtain a library of configurations of pure water molecules at equilibrium; and third to provide a reasonable starting point for the calculation including the glycine molecule. The pure water molecular dynamics calculations were carried out using 256 water molecules at constant energy with an average temperature of approximately 300 K. The initial configuration was that of a face centered cubic crystal and the density was equal to the experimentally determined density of liquid water. Periodic boundary conditions with the minimum imaging convention were employed [ 111. The water-water potential was truncated at a spherical

55

cutoff of 7.3 A. The average total internal energy per water molecule was - 6.99 kcal mol- I compared with an experimental energy of - 8.1 kcal mol - ’ [ 9 1. The average internal energy reported by Reimers and Watts [ lo] for a simulation using 64 water molecules was -7.4 kcal mol-I. This discrepancy of 6% may be due to slight differences in the intramolecular potential or, more probably, to the different number of water molecules used in the two simulations. In an earlier study using a somewhat different potential, Watts [ 121 found that increasing the number of molecules from 64 to 500 changed the internal energy from - 7.8 to - 6.6 kcal mol- ‘, respectively. The radial distribution functions for the hydrogen-hydrogen, hydrogen-oxygen, and oxygen-oxygen separations were indistinguishable from those obtained in ref. [9]. Unlike the calculation for pure water which used periodic boundary conditions with the minimum image convention, the water-glycine calculation was performed by treating the system as a cluster. The reason for this choice is as follows: our focus is on the vibrational structure of the glycine which depends primarily on the intramolecular glycine potential and secondarily on the water-glycine potential. The magnitude of the intramolecular glycine potential is of the order of (3N- 6 )RT/2 where N= 10, the number of glycine atoms and is equal to approximately 7 kcal mol- ‘. This value is less than 1% of the total waterwater intermolecular potential energy. In order to maintain conservation of energy to within a few percent of the glycine intramolecular potential, the water-water potential and the water-glycine must be calculated extremely accurately. The high degree of accuracy is necessary because forces and force constants of glycine depend on the first and second derivatives of the glycine intramolecular and waterglycine intermolecular potentials. The failure of energy conservation is a result of the presence of long range terms in the intermolecular potential functions. These long range terms are nonnegligible even beyond a distance of one-half the cell length. Since the maximum cutoff using periodic boundary conditions with the minimum image convention is one-half the cell length, we felt obliged to treat the system as a cluster. For the purposes of studying the vibrational structure of glycine, we believe that it makes little differ-

56

J.S. Alper et al. / Vibratronal structure of the solvated giyme zwlttenon

ence whether a cluster or periodic boundary conditions are used. The vibrational calculation involves the glycine intra- and intermolecular potentials only. In a Monte Carlo calculation of glycine in 200 water molecules by Roman0 and Clementi [ 13 ] in which the system was treated as a cluster, the water-glycine potential energy was found to be - 99.19 kcal mol- ‘. In a more recent Monte Carlo calculation involving 2 15 water molecules using periodic boundary conditions, Mezei et al. [ 141 found the water-glycine potential energy to be -99.28 kcal mol-‘. In addition, as will be discussed below, the water-glycine potential energy remains essentially constant when the number of waters is reduced from 253 to 200. We might also note that while structural changes in the cluster might be induced by the free boundary, these effects are probably small compared with the reorganization around the solute molecule. We removed three water molecules from the center of the cluster of 256 water molecules obtained from the pure water simulation in order to form a cavity in which to insert a glycine molecule. The choice of orientation of the glycine molecule in the cavity was arbitrary. The incorporation of the water-glycine interaction into the molecular dynamics procedure for pure water is completely straightforward; however, the inclusion of the intramolecular glycine potential presents a few technical difficulties. The intramolecular glycine potential is expressed in terms of the positions of the glycine atoms, the forces and the force constants obtained from GAUSSIANSZ. All of these quantities are defined with respect to the original glycine coordinate system. In the course of the simulation of the dynamics of the system, the glycine molecule suffers translations and rotations in addition to vibrations. The coordinates of the glycine atoms at each step (the molecular dynamics coordinates) must be transformed back into the frame (the original frame) with respect to which the original coordinates, forces, and force constants are defined. The transformation for the translational motion is trivial; the center of mass of the glycine molecule calculated at each step is translated until it coincides with the center of mass of the glycine in the original frame. The rotational motion requires more effort. At each step we determine the rotation matrix needed to transform the molecular dynamics coordinates into

the original frame. In the original frame, the Eckart conditions [ 1 ] must be satisfied. If we let R be the rotation matrix, tf the vector of the equilibrium coordinate of atom i in the original frame, r: the vector of the molecular dynamics coordinates, and m, the mass of the ith atom, then the rotation is determined by the vector equation C m,rpx(Rr:-rp)=O.

(6)

In essence this rotation transforms the glycine into the original frame because this is the frame in which the angular momentum is zero, i.e. the frame in which there is no rotation. Finite rotations involve complicated functions of the Euler angles. Infinitesimal rotations, on the other hand, take the simple form [ 15 ]

R mf

=

1 -d&

W 1

- d& de,

( de,

-de,

1

,

(7)

)

which is linear in the infinitesimal rotation angles d&, de,, and d& At the nth step of the MD calculation, we determine Rlnf by means of the equation R, =

R,,rR,- , 2

(8)

where R, is the total rotation required at the nth step. This approximation for R, is an excellent one because the d& are on the order of 1Oe3 rad. The recursive procedure defined by eq. (8 ) is initiated by noting that at the beginning of the calculation, the molecular dynamics and original frames coincide and thus R. is the identity matrix. The rotation matrix R, is then used to rotate the molecular dynamics glycine coordinates into the original frame enabling us to calculate the intramolecular potential energy. Next, the forces calculated in the original frame must be rotated back to the molecular dynamics frame in which the rest of the calculation is performed. However, we cannot simply use R - 1 to effect this transformation. In the original frame the Cartesian components of the force, F,, are given by J’,=J;-

C_&(s,-$‘), J

where i andj label the 30 Cartesian components the glycine in the original frame. The potential ergy is given by

(9) s of en-

J.S. Alper et al. / Vibrational structure of the solvatedglycine zwitterion

(10) where S’ is the Cartesian coordinate in the molecular dynamics frame. The components of the force in the molecular dynamics frame, F: , are derivatives of I’, with respect to the coordinates in the molecular dynamics frame. We have

(11) The first factor on the right hand side is F,. The second factor consists of the sum of two terms. The first is simply R. The second involves the derivatives of R with respect to each Cartesian component. The matrix R is the same for each of the glycine atoms, consequently it is convenient to change labels and consider the (Y,/I and y components of the various atoms labelled i and j in the MD frame, e.g., ra( i). Denoting the forces in the original and MD frames by F,(i) and F i(j) respectively, we find after some algebra: F;(i) = 1 F,(i)& OI (12)

We now have expressions for the contribution to the forces and potential energy resulting from the intramolecular glycine interactions. 3.2. Thermodynamic and structural results All our studies of the water-glycine cluster used a time step-size of 2.0~ 10m4 ps. This step size is roughly 20 times smaller than the period of the highest vibrational frequency of glycine ( = 3000 cm- ’ ) . Approximately 15000 steps were used to equilibrate the water-glycine system at an average temperature of 300 K, starting from the configurations of the water molecules obtained from the simulation of pure water, and then another 10000 steps at constant energy were used to obtain the equilibrium results used in this work. An additional 5000 steps were subsequently run for further verification that equilibrium had been reached. Over the course of the calculation, the total energy was conserved to a part in 10000. After the equilibrating 15000 steps, no monotonic trends were ob-

57

served in any of the components of the total energy, e.g. water-glycine potential energy. The fluctuations in the total potential energy were less than 1%, while the fluctuations in the water-glycine potential energy were of the order of 3.5%. The fluctuations in the total (kinetic plus potential) energy amounted to only a few percent of the intramolecular glycine potential energy. Average values and their standard deviations for the various components of the energy are given in table 2. We note here that our value for the average waterglycine potential, - 172 kcal mol-’ is considerably lower than those obtained by Roman0 and Clementi [ 13 ] and by Mezei et al. [ 14 1, approximately - 99.2 kcal mol- ‘. (For the purposes of this discussion, we will refer to these latter calculations as RCM. ) We can attribute this difference in potential energy to one or more of the differences between our calculations and theirs. First, the form of the water-glycine potential differs somewhat from that used by RCM. The potential used by RCM was based on a three parameter fit using only the first three terms of eq. (4). Second, the intermolecular water potential used in this calculation is also different. As mentioned above, we used the RWK2 potential of Reimers et al. [ 91, while RCM used one based on a SCF with CI calculation [ 131. Third, we are including the intramolecular motions of the glycine and the waters. These motions affect the positions of the glycine and water atoms and so change the value of the potential energy. Fourth, our system contains a somewhat larger number of water molecules, we used 253, Roman0 and CleTable 2 Glycine-water temperature (K) and energies culated using molecular dynamics Average temperature intramolecular water potential energy water-water potential energy intramolecular glycine potential energy water-glycine potential energy total potential energy total energy

(kcal mol-‘)

Standard

309.0 174.0

12.0 19.0

- 1569.0

44.0

6.2

3.7

- 172.0

8.0

- 1561.0 -853.7

28.0 0.2

cal-

deviation

58

J.S. Alper et al. 1 Vlbratlonalstructure of the solvated glycme zwrtterron

menti used 200, while Mezei et al. used 2 15. Finally, we used the conformation of glycine in which the CH2 group is staggered with respect to the NH3 group; RCM used the eclipsed form [ 13 1. The expe~mental crystal structure of the glycine zwitterion has the staggered geometry. Furthermore, our molecular dynamics calculations on each of the conformers indicated that the staggered geometry is more stable in solution. In order to compare our results with the previous calculations and to assess the contributions that each of the factors listed in the previous paragraph made to the difference in water-glycine potential energy, we performed various Monte Carlo calculations using rigid molecules. Our first calculation treated the rigid eclipsed conformer with 200 water molecules. The potentials were the same as those used by RCM. We performed our calculation starting from a configuration of 200 waters chosen from the con~guration of 253 waters obtained from our molecular dynamics calculation. The bond lengths and angles of the water molecules were reset to the equilibrium values. We obtained a water-glycine potential energy of - 105 kcal mol- ‘. A second calculation using all 253 waters resulted in a water-glycine potential energy of - 109 kcal mol- ‘. The next Monte Carlo calculation, again using 253 water molecules, replaced the RCM water intermolecular potential with the RWK2 potential. The water-glycine potential decreased to - 147 kcal mol-‘. Finally, the eclipsed configuration was replaced by the staggered, again using the RWK2 potential. The water-glycine potential decreased to - 157 kcal mol- ‘. Since the water-glycine potential calculated by MD is - 172 kcal mol- ‘, the difference of 15 kcal mol- ’ can be accounted for by the inclusion of the intramolecular potentials and the use of the 4-term instead of the 3-term water-glycine potential. We conclude that the choice of intermolecular water-water potential is the single most important factor in determining the water-glycine potential energy, accounting for more than one-half the difference between the RCM and our values of the waterglycine potential energy. No single one of the other factors is responsible for more than approximately 15% of the difference. Since the RWK2 water-water intermolecular potential was constructed to be an

improvement over earlier potentials, including the one used by RCM, we believe that our value for the water-glycine potential energy is more reliable than the earlier values.

4. Vibrational calculations and results The basic premise of this paper is that the effect of the solvent on the forces and force constants of the solute molecule can be estimated by means of corrections determined by a molecular dynamics calculation of the entire system. These corrections can then be added to the ab initio quantum mechanical values of the forces and force constants of the solute molecule alone. Once equilibrium in! the molecular dynamics procedure has been reaqhed, the total force and total force constant for each Cartesian component of a solute atom position are given by autot) as I

F,(total)=-

=F _ avlv, , -g-, I

(13)

and FJtot)

= g

(14) 1

J

where F, was defined in eq. (9 ). Note that the derivatives are evaluated in the glycine original frame. Thus the second term in each of eqs. ( 13) and ( 14) must be transformed from the molecular dynamics frame using the rotation matrix described above. As mentioned above, we expect that the solvent molecules stabilize the glycine zwitterion by reducing the average forces on each glycine atom, where the average Cartesian component of the force on each atom of glycine is the mean of the values calculated at each step of the molecular dynamics calculation. At equilib~um, the average force should be reduced to zero. A comparison of the ab initio forces calculated for the original isolated glycine zwitterion, which would be zero if the zwitterion were in an equilibrium geometry, with the final average forces determined by the molecular dynamics calculation shows that this expectation is largely realized. The absolute magnitude of the nonzero forces on the atoms of the original isolated glycine zwitterion are of the order of 0.03-0.3 mdyn. The final average forces, reflecting the contribution from the solvent, are on the order of

J.S. Alper et al. / Vibratzonal structure of the solvatedglycine zwitterion

100 times smaller, with no component being larger than 4.3 x 10m3 mdyn. We have also calculated the forces on the atoms of an isolated glycine zwitterion whose geometry is that of the average geometry obtained from the MD simulation. These forces are of the same order of magnitude as those for the original isolated glycine zwitterion, indicating again that the reduction in the forces is due to the interactions of the glycine with the solvent and not simply due to the use of a more correct geometry of the solvated zwitterion. The values of the bond lengths, bond angles, and torsion angles of glycine determined from the averages of the values of these quantities over each of the configurations of the solvated glycine molecule are listed in table 1. It is readily apparent that the geometry of the solvated glycine molecule has changed comparatively little from that of the original isolated form. The bond lengths have changed from their original values by a few hundredths of an Angstrom and the bond angles and most of the dihedral angles by a few degrees. The major difference in the geometry involves the COO group which is no longer in the plane defined by the CCN skeleton, but is now at an angle of approximately 20”. As we mentioned above, this torsional degree of freedom is associated with the normal mode in the isolated glycine zwitterion that involves the imaginary fundamental vibration frequency. We note that the crystalline geometry shows a similar dihedral angle, indicating that the intermolecular forces that stabilize the zwitterion in the solution and in the crystal may be similar. The vibrational frequencies were computed using a vibrational Born-Oppenheimer approximation. We assume that in the time equal to the period of the lowest frequency vibrational mode of the glycine molecule, the position and orientation of the glycine and each of the water molecules remain relatively constant. This assumption is satisfied since the intermolecular forces act to maintain the configuration of the system. The libration frequency in the water dimer is of the order of 150 cm-’ [ 81. The collective motions of the much larger number of water molecules surrounding the glycine zwitterion will be characterized by considerably lower frequencies. We therefore can divide the time-ordered sets of molecular dynamics instantaneous configurations into

59

groups such that the time of final configurations minus the time of the initial configuration of each group is equal to the lowest vibrational period of the glytine, namely z 1.4X lo-l3 s. We begin with the first group and calculate a set of average force constants from the force constants for each configuration of the group. We then calculate a set of fundamental frequencies for this set of force constants. This procedure is repeated for each successive group of configurations. Corresponding frequencies are then averaged to obtain the set of fundamental frequencies characteristic of glycine in aqueous solution. The vibrational calculation was entirely standard with the exception of the fact that the B matrix must be modified according to eq. ( 1) in order to include the effect of the nonzero forces. The B matrix was generated using a set of symmetrized internal coordinates. We chose the internal coordinates given by Destrade et al. [ 16 ] except that the definition of COO torsion used was that of Machida et al. [ 17 1. The internal force constant matrix and fundamental frequencies were calculated with the vibrational program of McIntosh and Peterson [ 18 ] suitably modified [ 19 ] to include the effect of the forces. We also performed the vibrational analysis assuming the forces were zero. We found that in no case did the values of a fundamental frequency calculated using the modified and unmodified B matrix differ by more than 0.5 cm-’ thus showing that the inclusion of the solvent in the molecular dynamics calculation has rendered the linear forces on the glycine atoms totally negligible. Moreover, all the fundamental frequencies are now real. The calculated fundamental frequencies for 4 isotopomers of glycine; undeuterated glycine, carbondeuterated glycine (NH: CD2 COO- ), undeuterated glycine in DzO (ND: CH,COO- ), and carbondeuterated glycine in D20 (ND: CD2COO- ); are listed in reverse numerical order in tables 3-6. We have found that for undeuterated glycine, the fret quencies obtained from the set of force constants found by averaging the instantaneous force constants over all the molecular dynamic configurations are in close agreement with the frequencies obtained by the method described above. The only exceptions are those high frequency stretches for which there is no experimental solution data. Because the amount of

J.S. Alper et al. / Mbratronal structure of the solvated gi_vclnezwttenon

60

Table 3 Calculated

and observed

Frequency

(cm-

frequencies

of NH: CH,COO-

’)

calculated scaled

3645 3619 3535 3144 3084 1969 1924 1828 1775 1663 1567 1545 1483 1352 1286 1110 1091 959 840 754 711 612 439 252

3390 3366 3288 2924 2868 1831 1790 1700 1651 1547 1457 1437 1379 1257 1196 1032 1015 892 781 701 661 569 408 234

aqueous

and observed

Frequency

(cm-’ )

frequencies

calculated

observed

ab imtio

Table 4 Calculated

a1

3011 2968 1615

1440 1410 1327 1118 1027 896 665 577 502

‘) Fromrefs. [21] and [22]. b, From table 7 of ref. [ 201 whrch is based on ref.

of ND: CH,COO-

observed

crystal b,

ab mitio

scaled

aqueous

overlap 3152 3058 3009 2973 1634 1618 1597 1533 1446 1417 1337 1315 1138 1119 1037 919 897 704 613 545 507 364 238

3145 3084 2683 2659 2530 1832 1661 1531 1483 1463 1423 1391 1351 1199 1122 1067 949 918 728 713 605 591 415 240

2925 2869 2495 2473 2353 1704 1545 1424 1379 1360 1323 1294 1257 1115 1043 992 882 854 677 663 563 550 386 224

3011 2969

a.bJ See footnotes

2201 1623 1438 1409 1320 1273 1200

1000 964 840 640 595 500

‘)

crystal b’ 3008 2972 overlap 2341 2203 1590 1440 1408 1322 1268 1193 1182 1168 1047 1004 966 823 772 674 598 493 395 337 232

for table 3.

[ 3 1.

computational time required for the repeated calculations of sets of vibrational frequencies for each of the isotopomers is extremely large compared with the time required using the one set of force constants obtained by averaging over all configurations, we report the results from the latter method. A comparison of the calculated with the experimental results is somewhat hampered because in aqueous solution, the absorption due to the water molecules, especially in the mid-infrared, interferes with that of the glycine. As a result, most of the spectroscopic data, especially those reported in the more recent studies, are based on the a-crystalline form of glycine [ 16,17,20]. In addition, the solution data we have found [ 2 1,22 ] are 25-30 years old and predates the advent of modern high-resolution laser Raman

and Fourier Transform IR spectrometers. Although the intermolecular glycine-water forces in solution certainly differ from the glycine-glycine interactions in the a-crystal, the fundamental frequencies found in solution and in the crystalline phase are remarkably similar. With few exceptions, corresponding frequencies differ by only a few wavenumbers so that it is useful to use both sets of results in our comparison. Because our focus is on the solvated glycine, our primary comparison is with the solution studies and we use the crystalline studies as a supplement. We report the experimental fundamental frequencies for the four isotopomers of glycine for both solution and crystalline phase in tables 3-6. At this stage of our work, we are unable to give a quantitative explanation for the fact that the experi-

61

J.S. Alper et al. / Vibrationalstructure of the solvatedglycine zwitterion Table 5 Calculated

and observed

Frequency

(cm-’

frequencies

) observed

calculated ab initio

scaled

3645 3619 3535 2341 2246 1969 1924 1824 1775 1504 1454 1418 1241 1179 1030 1019 938 929 835 729 624 599 434 233

3390 3365 3288 2178 2089 1831 1789 1696 1650 1399 1352 1319 1154 1096 958 948 872 864 777 678 580 557 404 217

a,b) See footnotes

of NH: CD2COO-

aqueous

2260 2161 1617

1404 1199 1107 1041 912 872 824 659 499

Table 6 Calculated

and observed

Frequency

(cm-‘)

frequencies

observed

calculated a)

crystal b)

ab initio

scaled

overlap 3151 not obs. 2257 2164 1644 1620 1587 1525 1407 1216 1188 1113 1044 928 911 873 803 675 567 530 472 359 234

2683 2659 2530 2340 2245 1827 1496 1426 1396 1352 1276 1232 1197 1079 1029 956 913 851 708 621 602 579 411 224

2495 2473 2353 2177 2088 1699 1391 1326 1298 1258 1187 1146 1113 1003 957 889 849 792 658 578 561 539 382 208

for table 3.

mental solution and a-crystal frequencies are in such good agreement with each other or for the fact that our calculated geometry for the solvated glycine is so close to the experimental geometry found for the cycrystal. Unless these results are purely coincidental, which seems unlikely, there must be some underlying similarity in the environments created by the solvent and by the crystal. These effects seem worthy of further investigation. Vibrational frequencies calculated for gas phase molecules using SCF methods are generally lo- 15% too high [ 23 1. The calculated ab initio frequencies for the glycine zwitterion in solution also show this degree of divergence from the experimental aqueous glycine data. With just a few minor exceptions, the calculated vibrational assignments agree with the experimental ones. A more detailed discussion of the assignments will be given in a paper now in progress.

a.b) See footnotes

of ND: CD2COO-

aqueous

2266 2179 1610 1402 1200

1126 1073 950 867 812 716 600 480

a1

crystal b’ 2349 2340 2263 overlap 2177 1580 1405 1197 1176 1169 1135 1077 1011 957 921 872 804 711 650 534 468 393 335 222

for table 3.

The average absolute per cent deviation of the calculated from the 56 aqueous experimental frequencies is 14.2Oh. If the force constant matrix is multiplied by a factor of 0.865, corresponding to scaling the frequencies by a factor of 0.930, the agreement with experiment is considerably improved. The average absolute percent deviation of the calculated from the aqueous frequencies becomes 7.3%. Similar agreement (7.4%) was found when the calculated frequencies were compared with the experimental crystal frequencies.

5. Conclusions Ab initio calculations of the gas phase glycine zwitterion show the presence of an unstable vibrational mode. The NH: and COO- groups twist toward each

62

J.S. Alper et al. / ?‘lbratlonal structure of the solvated gl.vcme zwlttenon

other in an attempt to form hydrogen bonds. This motion is reflected in the forces on the various glytine atoms that give rise to an imaginary fundamental frequency corresponding to the unstable mode. In aqueous solution, the water-glycine intermolecular interactions stabilize the zwitterion. The resulting structure of the zwitterion is a staggered conformer with the plane of the COO- group rotated by an angle of approximately 20” from the plane formed by the nitrogen and 2 carbon atoms of the molecule. This torsional angle is also seen in experimental studies of the a-crystalline form of the zwitterion. The molecular dynamics calculations show that after equilibrium is reached, the average force on each glycine atom due to the other glycine atoms and the surrounding water molecules is negligible. All the vibrational modes of glycine are now stable and the fundamental frequencies are all real. The agreement between the calculated and experimental aqueous frequencies is of the same order as that typically obtained by ab initio calculations of the vibrational spectra of gas phase molecules.

References [I ] E.B. Wilson Jr., J.C. Decius and P.C. Cross, Molecular Vibrations (Dover, New York, 1980). [2] P. Pulay, in: Modem Theorettcal Chemistry, ed. H.F. Schaefer III (Plenum Press, New York, 1977 ) [ 31 P.G. Jonsson and A. Kvick, Acta Cryst. B 28 ( 1972) 1827.

[ 41 J.S. Binkley, M.J. Fnsch, D.J. DeFrees, K. Raghavachari, R.A. Whiteside, H.B. Schlegel, E.M. Fluder and J.A. Pople, GAUSSIANIZ, Revision H Version (Carnegie-Mellon, Pittsburgh, 1986). [ 51L. Carozzo, G. Corongiu, C. Petrongolo and E. Clementi, J. Chem. Phys. 68 (1978) 787. [6] J.R. Reamers and R.O. Watts, Mol. Phys. 52 (1984) 357. [ 71 D.F. Coker, R.E. Miller and R.O. Watts, J. Chem. Phys. 82 (1985) 3554. [8] D.F. CokerandR.0. Watts, J. Phys. Chem. 91 (1987) 2513. [ 9 ] J.R. Reimers. R.O. Watts and M.L. Klein, Chem. Phys. 64 (1982) 95. [lo] J.R. Reimers and R.O. Watts, Chem. Phys. 91 ( 1984) 201. [ 111 M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids (Oxford, Umv. Press, Oxford, 1987). [ 12 ] R.O. Watts, Mol. Phys. 28 ( 1974) 1069. [ 13 ] S. Roman0 and E. Clementt, Intern. J. Quantum Chem. 19 (1978) 839. [ 14 ] M. Mezei, P.K. Mehrotra and D.L. Beveridge, J. Biomol. Struct. Dyn. 2 ( 1984) 1. [ 151 H. Goldstem, Classical Mechamcs, 2nd Ed. (AddisonWesley, Reading, 1980). [ 161 C. Destrade, C. Garrigou-Lagrange and M.-T. Forel, J. Mol. Struct. THEOCHEM 10 ( 197 1) 203. [ 17 ] K. Machida. A. Kagayama, Y. Satto, Y. Kuroda and T. Uno, Spectrochim. Acta A 33 ( 1977) 569. [ 181 D.F. McIntosh and M.R. Peterson, QCPE 11 (1977) 342. [ 191 M.A. Lowe, J.S. Alper, R. Kawiecki and P.J. Stephens, J. Phys. Chem. 90 (1986) 41. [20] M. Kakihana, M. Akiyama, T. Nagumo and M. Okamoto. Z. Naturforsch. 43a ( 1988) 774. [ 211 M. Takeda. R.E.S. Iavazzo, D. Garfinkel, I.H. Scheinberg and J.T. Edsall, J. Am. Chem. Sot. 80 ( 1958) 3813. [22] S.A.S. Ghazanfar, D.V. Myersand J.T. Edsall, J. Am. Chem. Sot. 86 (1964) 3439. [ 231 G. Fogarasi and P. Pulay, Ann. Rev. Phys. Chem. 35 ( 1984) 191.