Methyl group dynamics in paracetamol and acetanilide: probing the static properties of intermolecular hydrogen bonds formed by peptide groups

Chemical Physics 244 Ž1999. 49–66 www.elsevier.nlrlocaterchemphys

Methyl group dynamics in paracetamol and acetanilide: probing the static properties of intermolecular hydrogen bonds formed by peptide groups M.R. Johnson

a,)

, M. Prager b, H. Grimm b, M.A. Neumann C.C. Wilson d

a,c

, G.J. Kearley a ,

a

Institut Laue-LangeÕin (ILL), B.P. 156, 38042 Grenoble Cedex 9, France Institut fur des Forschungszentrums Julich, D-52425 Julich, Germany ¨ Festkorperforschung ¨ ¨ ¨ Laboratoire de Spectrometrie ´ Physique, UniÕersite´ J. Fourier Grenoble, CNRS (UMR 5588), B.P. 87, 38402 Saint-Martin-d’Heres ` Cx, France d ISIS, Rutherford Appleton Laboratory, Chilton, Didcot, OX11 0QX, UK b

c

Received 16 November 1998

Abstract Measurements of tunnelling and librational excitations for the methyl group in paracetamol and tunnelling excitations for the methyl group in acetanilide are reported. In both cases, results are compared with molecular mechanics calculations, based on the measured low temperature crystal structures, which follow an established recipe. Agreement between calculated and measured methyl group observables is not as good as expected and this is attributed to the presence of comprehensive hydrogen bond networks formed by the peptide groups. Good agreement is obtained with a periodic quantum chemistry calculation which uses density functional methods, these calculations confirming the validity of the one-dimensional rotational model used and the crystal structures. A correction to the Coulomb contribution to the rotational potential in the established recipe using semi-emipircal quantum chemistry methods, which accommodates the modified charge distribution due to the hydrogen bonds, is investigated. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Methyl group dynamics; Intermolecular hydrogen bonds; Peptide groups

1. Introduction Methyl groups provide the most widespread example of coherent proton tunnelling in the solid state due to the indistinguishability of the trio of protons

)

Corresponding author. Tel.: q33-76207139; Fax: q3376483906; E-mail: [email protected]

and the small moment of inertia of the molecular rotor w1,2x. The approximately exponential dependence of the tunnel effect on the height of the potential barrier makes the methyl group a highly sensitive probe of its chemical environment. In this context, hundreds of tunnelling measurements have been accumulated, empirical rules being established by manipulating the chemical composition of samples. For example, by studying series of alkanes w3x,

0301-0104r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 1 - 0 1 0 4 Ž 9 9 . 0 0 1 3 8 - X

50

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

ketones w4x and carboxylic acids w5x, variations in tunnel frequencies have been interpreted in terms of the crystal packing of different length chains, while in trihalogenated mesitylenes w6x, increasing the size of the halogen atoms from chlorine through to iodine results in a weakening of the rotational potentials. More recently however, the advent of powerful desktop computers and the proliferation of force-field and quantum chemistry based molecular dynamics software packages have facilitated the numerical investigation of tunnelling results based on atom–atom interactions. Early calculations have been repeated w7–12x and new systems studied w13–18x to give a series of compounds for which a correlation is sought between observed and calculated methyl group observables. By using high level ab initio calculations on single molecules combined with Van der Waals ŽVDW. and Coulomb interactions to describe intermolecular interactions, the rotational potential has been calculated to a precision of ; 80%, tunnel frequencies are correct to within a factor of three, for a limited series of compounds w19x. Thus, in practice, the methyl group is now a quantitative probe of its environment. Since the molecular mechanics ŽMM. calculations are based on observed crystal structures, the limit to the number of systems Ž; 10. studied is due directly to the lack of crystallographic data at liquid helium temperatures: whereas there have been hundreds of tunnelling measurements the corresponding structural data is counted in tens. ŽThe number of systems which has now been included in the correlation, including the samples studied here, is ; 20.. The sample which works least well in terms of the calculated amplitude of the rotational potential, a value of ; 50 meV instead of ; 20 meV being obtained, is acetamide ŽCH 3 CONH 2 . w7x, one of the simplest molecules containing the peptide group –CONH–. Since peptide groups are responsible for hydrogen bond networking, the approximation used in our MM calculations, that a crystal can be decomposed into single molecules connected by weak intermolecular interactions, is less accurate when peptide groups are present. In assessing the extent to which the intermolecular hydrogen bonding compromises the accuracy of the MM calculation, it is essential to study other molecules containing the peptide link, namely in this article paracetamol Ž p-

hydroxy acetanilide: CH 3 CONHC 6 H 4 OH. and acetanilide ŽCH 3 CONHC 6 H 5 .. The molecules studied here are rare examples of systems where the low temperature crystal structures are known and the methyl group tunnelling dynamics have not been measured previously. Paracetamol has been the subject of recent single crystal diffraction work w20x and measurements of tunnelling and librational excitations by inelastic neutron scattering ŽINS. are reported here. The structure of acetanilide is known at 15 K w21x, the vibrational spectrum has been measured and force field based molecular dynamic simulations have been performed w22x, here we report the tunnelling measurement which is at the current limit of INS resolution. Small molecules containing peptide groups have often been studied as model systems for polypeptides and a-helix proteins, notably acetanilide w23x and n-methyl acetamide w24x. However, in both these systems, the amide vibrations seem to show anomalous behaviour. In acetanilide, the temperature dependence of some of the amide vibrations shows strong anharmonicity and unexpected sidebands, while in n-methyl acetamide, anomalous intensity in INS vibrational spectra at 20 K has resulted in the N–H stretch frequency being assigned at 1575 cmy1 rather than ; 3000 cmy1 . In the latter case, an ionic model of the hydrogen bond ŽN dy . . . Hq . . . O dy . is invoked instead of a covalent representation ŽNH . . . O. since the assignment of the stretching frequency implies a significant weakening of the N–H covalent bond. This in turn results in a lengthening of the N–H separation and the peptide structure is intermediate between amide-like ŽOCNH. and imidol-like ŽHOCN.. Since in the samples studied here the methyl group is bonded to the peptide group, the tunnelling dynamics of the methyl rotor are influenced by the electronic wavefunction in the hydrogen bond. It is the purpose of this article to show how the methyl group is a quantitative reporter of the hydrogen bonds. The article is set out as follows: Section 2 establishes the theoretical background for analysing methyl group dynamics, Section 3 describes the spectroscopic measurements on paracetamol and acetanilide, Section 4 brings together the essential crystallographic information for paracetamol and acetanilide, and in Section 5 calculations on both systems are

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

described. For the remainder of the article, paracetamol will be referred to as PA and acetanilide as AN.

2. The methyl group Hamiltonian, experimental observables and determination of the rotational potential The Hamiltonian describing single particle motion ŽSPM. of a rigid methyl group hindered by a potential V is: "2 E2 Hsy

2 I Ew 2

q V SPM

Ž 1.

where I is the moment of inertia of the rotor and w is the rotation angle. " 2r2 I is the rotational constant B. The potential V can be represented as a Fourier expansion: V SPM Ž w . s Ý b k e i3 k w .

Ž 2.

k

The potential has at least three-fold symmetry due to the indistinguishability of protons, that is the rotation R of a rigid, symmetric rotor through 1208 generates an energetically equivalent molecular configuration. It is normally sufficient to retain only Fourier coefficients b k up to second order Žy2 -s k -s 2., the first order terms generally dominate and an initial interpretation of spectroscopic data can be based upon a simple potential of the form V3 cosŽ3 w .. The eigenstates c of H have either A or E symmetry depending upon their transformation under R as follows: R c A s c A , R c E a s ´c Ea , R c E b s ´ Uc E b with ´ s e i2 p r3

Ž 3.

H of Eq. Ž1. can be developed in a basis of free rotor states, diagonalisation of the Hamiltonian matrix giving wavefunctions c which are linear sums of free rotor states:

c X s Ý c k e iŽ3 kqs. w

Ž 4.

k

The value of s depends on the symmetry X of the wavefunction. s is equal to 0, 1, y1 for symmetry A, Ea and E b , respectively. The energy levels consti-

51

tute a set of librational levels, each of which is composed of an A singlet and an E doublet, the energy difference between these components being the tunnel splitting hn i of that level, where the ground state tunnel splitting is hn 0 . The average energy difference between rotational levels in the ground librational state and those in a higher librational state defines the librational energy Ei . Spectroscopic observables of the methyl group are therefore the ground state tunnel splitting hn 0 and librational excitations. The equilibrium orientation of the rotor wmean , as defined by the position of the minimum of V in Eq. Ž2., is an observable to be compared with diffraction measurements. When V does not have a simple form, the equilibrium orientation is obtained from the wavefunction, where the integral ŽEq. Ž5.. extends between successive minima of V. min 1

Hmin 2 w P < c Ž w . <

wmean s 3

2

dw.

Ž 5.

Within the SPM approximation, the authors have established a numerical method which reproduces observed tunnel splittings from low temperature crystal structures w19x. The essential details of the calculation of V SPM Ž w . are as follows. Ž1. The molecular and crystal geometry from crystallographic measurements are retained, only the triangle of methyl protons is symmetrised to ensure C3 symmetry in the SPM approximation. When the structure determination is by X-rays and the C–H bond lengths and H–C–C angles are unrealistic, average values from neutron measurements are imposed. Ž2. The internal rotational potential is determined using ab initio Hartree–Fock methods as the variation in the ground state electronic energy with methyl group orientation. The 6-31GUU basis set and second order perturbation corrections ŽMP2. in the GAMESS w25x program are used. Ž3. Coulomb interactions with neighbouring molecules are calculated using the electrostatic potential derived ŽESPD. charges from the series of ab initio calculations, which depend on the methyl orientation. Ž4. VDW interactions have the physically reasonable exponential form for the repulsive term, the parameters are taken from the Universal Force Field w26x.

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

52

Ž5. Cut-offs for Coulomb and VDW interactions are taken big enough to ensure convergence. The SPM approach, including the calculation of the rotational potential as outlined in Ž1. – Ž5., will be referred to as the standard calculation for the remainder of the article. It is the determination of potential energy curves Žor surfaces. by displacing atoms along defined trajectories which we refer to here as MM. Each crystallographically distinct methyl group can give rise to a single tunnel splitting within the SPM approximation. When more tunnelling excitations than methyl groups are observed, the SPM approach must be developed to include coupling to other degrees of freedom. For methyl acetate w14x and durene w15x, the SPM Hamiltonian ŽEq. Ž1.. has been extended to describe an ensemble of coupled methyl groups by incorporating pairwise coupling potentials in the total Hamiltonian, higher order coupling, for example between triplets of rotors, is assumed to be negligible. E2 H s y Ý Bi i

Ew i2

q V Ž w 1 , . . . , wn . .

Ž 6.

The potential function V can be separated into one and two dimensional terms. n

V Ž w 1 , . . . , wn . s

n

Ý Vi s Ž wi . q Ý is1

Vi ,pj Ž w i , w j .

‘true’ two dimensional terms. The static potentials V s now contain all of the one-dimensional terms and are related to the single particle potentials V SPM by: Vi SPM Ž w i . s Vi s Ž w i . q

Ý Vi ,pj Ž wi , wj,SPM . .

Ž 9.

j/i

Here the angles w j,SPM denote the orientations of the space fixed neighbouring methyl groups in the SPM. In practice, the static potentials V s are calculated from V SPM and V p. Intermolecular, Coulomb and VDW contributions for V p are determined as described in Ž1. – Ž5. above for the SPM. For lithium acetate w27,28x, a model of a coupled pair of methyl groups was investigated using the Hamiltonian ŽEq. Ž6.. and potential energy calculations as described above, and was found to be inappropriate. An alternative model of rotation–translation coupling, in which the centre of mass of the rotor precesses during rotation of the methyl group, was elaborated w16x, and gives excellent agreement with a wealth of crystallographic and spectroscopic data. The model involves a two-dimensional harmonic oscillator coupled to a one-dimensional rotor and the parameters in the corresponding Hamiltonian are derived by fitting the potential energy surface calculated from the crystal structure Ž1–5 above.. Thus we are equipped with three models for investigating methyl group rotation: SPM, coupling between methyl groups and rotation–translation motion.

i , js1,i-j

Ž 7. Vi ,pj Ž w i , w j . s Ý k

Ý bk ,l ,i , j e i3 k w e i3 l w . i

j

Ž 8.

3. Spectroscopic measurements and preliminary interpretation

l

The pair potentials V p of two methyl groups i and j is obtained by calculating the potential energy of the pair of rotors for orientations w i s 108h1 and w j s 108h 2 with h1 , h 2 s 0, . . . ,11. Fourier coefficients b k ,l ,i , j are adjusted to reproduce the calculated potential energy. As for SPM potentials, only Fourier coefficients up to second order Žy2 -s k,l -s 2. are retained. The potential energy terms corresponding to the coefficients b 0, l ,i , j and b k,0,i , j are one-dimensional and are absorbed in the static potentials for clarity, that means b 0, l ,i , j s b k,0,i , j s 0 is set in V p without loss of generality and V p contains only

3.1. Paracetamol Commercially available PA Ž4-acetamido-phenol or p-hydroxy AN. from Sigma was used without further purification. The material was filled into a flat 1-mm thick sample can and put into a variable temperature cryostat. The sample was oriented at 458 to the beam. The backscattering spectrometer BSS1 at the DIDO research reactor of the research centre FZJ in Julich, Germany w29x, was used to measure the tun¨ nel spectrum. The neutron wavelength of this instru-

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

˚ . is determined by the siliconŽ111. ment Ž l s 6.27 A monochromator and analysers, the energy resolution is 1 meV ŽFWHM.. Spectra were taken within the energy range "13 meV at five sample temperatures between 5.4 K ŽFig. 1. and 35 K. Fig. 1 clearly shows the tunnelling transitions at "3.10 meV. The spectrum is described by the scattering function of a tunnelling methyl group, a purely elastic contribution from the other atoms of the molecule, especially its protons, and a flat background. This scattering function is convoluted numerically with the resolution function of the instrument, the solid line of Fig. 1 is the result of such a fit. The intensity of the tunnel transition relative to the elastic line is 1:11 and agrees well with 1:10.5 expected for a unique methyl group on the elastic background of the remaining six protons of the molecule at the largest momentum ˚ y1 w30x. transfer Q s 1.88 A The spectra obtained for higher sample temperatures have been fitted by constraining the relative intensities of the elastic and inelastic lines to be unchanged. With increasing temperature, however, a quasielastic line emerges from the elastic line with a width not predicted by theoretical models, but which is smaller than that of the tunnel transition w31x. This will change the relative elastic and inelastic intensity at high temperature when the classical limit is approached and thus influence the fit parameters. Furthermore, due to the Lorentzian wings of the resolution function, the tunnel line is not well-separated from the elastic peak and the derived fit parameters, especially at higher temperatures, have large error

Fig. 1. Tunnel spectrum of paracetamol measured at 5 K Ž Qs1.6 ˚ y1 . at BSS1, Julich. A ¨

53

Fig. 2. Temperature dependence of the position and width of the tunnel transitions in paracetamol.

bars. The position hn and width G of the tunnelling lines obtained this way is shown in an Arrhenius plot in Fig. 2. From the slope of the linewidth, an activation energy of 11 meV is derived, the shift shows a similar slope. Once the structure of the tunnel spectrum has disappeared, incoherent reorientation of the methyl group gives rise to a quasielastic signal which quickly becomes too broad to measure on a backscattering spectrometer. The decay of the elastic intensity was however recorded in the temperature range 50–200 K. A monotonic decrease m s d IrdT in elastic intensity is observed due to a Debye–Waller effect dominated by the three methyl protons. These measurements were completed in the regime of lattice modes by spectra from the spectrometer SV22, a thermal time-of-flight instrument at the w29x. Density of states measame reactor in Julich ¨ surements were taken at two sample temperatures, T s 20 K and T s 90 K. Fig. 3 shows the sum over all detectors in the range of scattering angles 708 2 u - 1308 after subtraction of a background and transformation from time-of-flight to energy scale. The spectra show two strong peaks at about 11.5 meV and 17.5 meV and a weaker peak at 26 meV. These features are well-defined only at low temperature, they have similar temperature dependences and are strongly damped at higher temperature. Interpretation of the spectroscopic data is based initially on the measured tunnelling transitions since

54

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

Fig. 3. Density of states spectra of paracetamol in the - 30 meV energy transfer range measured at SV22, Julich. ¨

these excitations occur at energies where coupling to other excitations is uncommon and a SPM approach is often sufficient. For a rotational constant of 0.672 meV Žderived from the diffraction data., the tunnel splitting of 3.10 meV determines a three-fold rotational potential of height Ž V3 s. 42 meV. According to theory w32x, the temperature-dependent tunnelling linewidth should show an activation energy E01 equal to the libration energy due to resonant coupling to phonons. The measured activation energy of 11 meV is in good agreement with the first peak in the librational spectrum. However, the libration energy calculated from V3 is 14 meV, this is approximately the average of the first two peaks in the librational spectrum. The first overtone in the V3 potential occurs at Ž E02 s. 26 meV which coincides with the weak peak in the Fig. 3 and shows the potential wells to be approximately harmonic Ž2 E01 –E02 .. Thus, the spectral features between 10 meV and 20 meV cannot be explained by an SPM model of methyl rotation and a ‘reasonable’ rotational potential. Structure could arise in the librational spectrum due to the presence of methyl–methyl coupling w15x, although no evidence of such coupling is seen in the tunnelling spectrum, or coupling to a lattice phonon of similar symmetry. If as a result of coupling the librational mode shows dispersion, the form of the librational band in a density of states measurement depends on the shape of the dispersion curve, for example, two peaks Žvan Hove anomalies. could occur corresponding to flat regions ŽEwrEk ; 0. at the ends of the dispersion curve. A quantum mechan-

ical model has been presented for an analogous system, 2-butine w33x. Single crystal measurements of the dispersion of the librational excitation andror reliable lattice dynamics calculations are required to test this idea. Alternatively, there may simply be other molecular vibrations in the energy range under consideration, which may or may not be coupled to the methyl group libration. The amide group is not ‘rigid’ and soft vibrational modes are commonly observed: in n-methyl acetamide w24x ŽCH 3 CONHCH 3 . two methyl librations and a torsional mode Žrotation about the C–N bond. are assigned at 21.9 meV, 20.5 meV and 19.3 meV, respectively, while n-alanine w34x ŽCH 3 CONHCHŽCH 3 .CONHCH 3 . shows a continuous distribution of peaks up to 60 meV. However, while AN w22x, which is a molecule closely related to PA, does have peaks in the vicinity of the libration, they are of much lower intensity. In addition, thermal factors for the methyl group carbon atom, derived from the single crystal diffraction measurements, do not show any anomalies which are indicative of large amplitude soft modes. The linear decrease of the elastic intensity, slope m, is related to a linear increase, slope S s d² u 2 :ŽT .rdT of the mean-square displacement of the methyl protons with temperature. In a harmonic oscillator model, which is reasonable for the proposed three-fold potential, the slope S and the mean-square displacement at zero temperature are simply related ² u 2 : 0 s Ž " vr2 k . S. With the mea˚ 2rdeg one obtains sured slope S ; 3.2. 10y3 A 2 ˚ 6² u : 0 s 0.20 A corresponding to an angular extension of 118 HWHM. This agrees perfectly with the zero point mean square amplitude ² w : s 118 HWHM derived from diffraction w20x. Spectroscopy gives a consistent result. Tunnelling and spectroscopy are described by a three-fold rotational potential. At T s 0 K, the proton density is derived from the calculated ground state wavefunction as cc U , giving 11.58 HWHM. 3.2. Acetanilide The sample was purchased in powder form from Aldrich. A flat sample, thickness 1 mm, aligned with the scattering angle of 458, was maintained at 2 K in a standard orange ŽILL. cryostat. The tunnelling

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

55

4. Essential structural information and hydrogen bond networks

Fig. 4. Comparison of the tunnel spectrum Žcrosses. of acetanilide and the instrument resolution function Žsolid curve., T s 2 K, ˚ y1 ŽIN16,ILL.. The difference spectrum Ždashed curve. Qs1.3 A shows the tunnel peaks, the tunnel splitting is 0.30"0.1 meV.

spectrum was measured on the backscattering spectrometer IN16 at the ILL, Grenoble w35x using polished siliconŽ111. monochromator and analysers Ž l ˚ . which constitutes the highest resolution s 6.27 A set-up ŽFWHMs 0.4 meV.. It is also the lowest flux set-up, the incident beam is approximately seven times less intense than in the standard set-up using an unpolished monochromator. Fig. 4 shows the measured spectrum Žcrosses. and the resolution function Žsolid curve. with the same peak heights. Excess intensity is seen in the flanks of the elastic peak and the difference spectrum shows clearly the tunnel peaks. The tunnel splitting, 0.3 " 0.1 meV, is consistent with the libration frequency of 18.0 meV, taken from the published vibrational spectrum, and a three-fold rotational potential of amplitude 65 meV Žrotational constant 0.664 meV as derived from the diffraction data.. The ratio of inelastic to elastic intensity is 1:12, smaller than for PA due to the ˚ y1 . and the lower momentum transfer Ž Qave s 1.3 A presence of more Bragg peaks in the detectors Žthe unit cell of AN is twice as big as that of PA.. The resolution function is a distorted Lorentzian since the monochromator, composed of thick strain-free crystals has a higher resolution than the analysers. Due to this and the fact that the tunnel peaks are not resolved from the elastic peak, fitting of the spectrum by numerical convolution of an elastic and two tunnel peaks with the resolution function was unstable and did not improve the above analysis based on the difference spectrum.

Crystallographic information for both samples is assembled in Table 1. In particular we note that the unit cell of PA is about half the size of that of AN, there being four and eight molecules, respectively, in the unit cells. Both structures have been measured using single crystal neutron diffraction techniques at low temperature ŽT s 15 K. w21,20x. Fig. 5 shows the smallest clusters of molecules for each system which include all the hydrogen bonds formed by a central molecule. PA forms hydrogen bonds with four neighbours, the hydrogen bonds linking the peptide groups with hydroxy groups. On the other hand, AN forms hydrogen bonds with only two neighbours but these bonds link peptide groups together giving rise to poly-peptide chains.

5. Calculation results 5.1. The standard calculation Fig. 6 shows the total rotational potentials, together with the intramolecular component and the intermolecular VDW and Coulomb components, for both samples, derived using the standard calculation. The origin of the rotation coordinate corresponds to a methyl group C–H bond being parallel to the peptide C5O double bond, positive angles denote clockwise rotation when looking from the methyl group to the molecule. In Table 2, the corresponding calculated observables are compared with Table 1 Crystallographic parameters for PA and AN

Temperature of measurement ŽK. Space group Z ˚. a ŽA

PA

AN

15 P 21 r c 4

15 Pbca 8

12.667

19.466

˚. b ŽA ˚. c ŽA a b g

9.166

9.332

7.073 90 115.51 90

7.735 90 90 90

˚ 3. Volume ŽA

741.2

1405.1

56

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

Fig. 5. Molecular clusters of Ža. PA and Žb. AN extracted from the measured crystal structures to show the intermolecular hydrogen bonds Ždashed lines. in these systems.

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

57

Fig. 6. Rotational potentials from the standard calculation for Ža. PA and Žb. AN.

the measured values, the rotational constant B is 0.672 meV for PA and 0.664 meV for AN, as determined from the crystal structures. Three-fold rotational potentials have been derived from the Table 2 Results of the standard calculation

measured tunnel frequency ŽmeV. standard calculation tunnel frequency ŽmeV. measured libration frequency ŽmeV. standard calculation libration frequency ŽmeV. measured equilibrium orientation standard calculation equilibrium orientation experimental V3 potential ŽmeV. standard calculation barrier height ŽmeV. a

PA

AN

3.10 0.998 13.8 a 18.4 75.88 88.68 41 48

0.30 1.15 18.0 14.5 54.78 46.58 65 56

This frequency does not correspond a unique spectral peak but is the average of two frequencies Žsee text Section 3.1..

measured tunnel frequencies and the amplitudes of the calculated rotational potentials are given to quantify the error in the calculated potential. In both cases, there are significant six-fold contributions to the standard calculation rotational potentials, so fitting with V3 potentials is unreasonable and the quoted potential amplitudes serve only for rough comparison. Nevertheless, the level of agreement is generally very poor. For PA, the amplitude of the barrier is overestimated by ; 7 meV while for AN the barrier is underestimated by ; 9 meV. In these systems, the equilibrium orientation, as well as the tunnel frequency, is a sensitive probe of the intermolecular interactions, since the intramolecular potential is small Ž- 7 meV. compared to the other contributions Ž) 20 meV.. For an isolated molecule, the electrostatically favoured orientation of the methyl group corresponds to 08 and deviations from this in

58

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

the crystal, due to intermolecular interactions ŽVDW, Coulomb and hydrogen bonds., are large Žy75.88 for PA and 54.78 for AN.. The errors in these orientations are significant in both cases, q12.88 for PA and y8.28 for AN. 5.2. Sources of error There are three possible sources of error in any such calculation of the rotational potential: the crystal structure may be incorrect, the rotational trajectory imposed may not be appropriate or the calculational techniques employed may not be sufficiently accurate. In both cases, the crystal structures are beyond reproach, they have been measured recently by neutron diffraction at 15 K and the protons are well-localised at temperatures comparable to the tunnelling measurements. The rotational trajectory imposed in this study is the simplest, corresponding to SPM. The need to consider complex rotational dynamics is estimated from potential energy calculations. Methyl–methyl coupling potentials have been calculated for both systems, the biggest coupling potential Ž V pŽm ax. ŽPA. s 5.73 meV, V pŽm ax. ŽAN. 1.10 meV. varies from ; 2% ŽAN. to ; 11% ŽPA. of the static potential Ž V s ŽPA. f 53 meV, V s ŽAN. f 50 meV.. The coupling potential in PA causes a splitting of the librational band of 0.7 meV, it should be an order of magnitude greater in order to explain the observed splitting of 6 meV. This result is consistent with coupling studies in other systems; the strongest coupling potential must be comparable to the static potential if spectral structure is to be observed w15x. In the two systems studied here methyl–methyl coupling effects can be neglected. Rotation–translation coupling was alluded to in the original calculations on acetamide as a way to attenuate the rotational potential w7x. However, the common feature of the three systems under investigation is the hydrogen bond network formed by the peptide groups which has a tendency to stabilise the crystal structures. Any precession of the centre of mass of the methyl group would be accompanied by stretching of the hydrogen bonds and this is energetically unfavourable. In addition, rotation–translation motion generally manifests itself in anomalous thermal factors for the methyl group atoms but reasonable anisotropic thermal parameters have been refined in both cases.

The calculational techniques employed have been shown to be valid for about 20 systems which do not include intermolecular hydrogen bonds. It is precisely the existence of such hydrogen bonds which compromises the distinction here between inter- and intramolecular interactions and thus the treatment only of single molecules with ab initio techniques. The contributions to the rotational potential as shown in Fig. 6 reveal that the Coulomb term is comparable to the VDW term while the intramolecular term is the least important Ž Vintra - 7 meV.. In view of this and the fact that the methyl rotation is independent of the hydrogen bond dynamics in these MM calculations, that is the interaction along the hydrogen bond is not solicited, solutions to be discussed in the next section focus on the electronic Ži.e., charge. distribution in the hydrogen bonds. 5.3. Extending ab initio calculations to more than one molecule 5.3.1. Density functional calculations on a periodic system Molecules with about 20 atoms are already at the current computational limit for the Hartree–Fock methods implemented in GAMESS-UK using UNIX workstations and extending this method to dimers or clusters is not tractable. High level quantum chemistry calculations are possible on periodic systems employing either Hartree–Fock methods w36x or density functional techniques ŽDFT.. It is the latter, as implemented in CASTEP w37x in the Cerius 2 1 program suite and run on a Silicon Graphics workstation and a parallel version of the same code ŽCETEP. run on a Cray T3E, which has been used here. In this application of density functional methods, the electron density, which describes the ground state of the crystal, is derived by the variational method from a basis set of plane waves which fill the volume of the unit cell. The computational cost therefore depends on the volume of the unit cell as well as the number of atoms therein. On a Silicon Graphics Indigo2 ŽR10000 processor., it was only possible to perform calculations on PA using the ‘medium’ basis set Žplane wave energy cut-off 600 eV.. Even using the

1

Cerius 2 is a product of Molecular Simulations.

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

coarse basis set Ženergy cut-off 200 eV., AN is intractable. Both systems were therefore studied using CETEP on a Cray T3E where it was possible in the CPU time available to investigate basis set convergence. The generalised gradient approximation ŽGGA. was used for the exchange-correlation term in the energy functional, which is expressed in terms of electron density rather than a many-electron wavefunction, GGA being more accurate than the commonly used local density approximation ŽLDA. for regions of low electron density. Single k-point sampling for integration over the Brillouin zone was used, the default for unit cells of this size. More k-points have only been required for much smaller systems, such as nitromethane Žunit cell volume ˚ 3 . and methyl iodide Ž V s 325 A˚ 3 ., in V s 275 A ˚ The which the shortest unit cell dimension is F 5 A. results are assembled in Table 3. The potential parameters clearly converge with increasing basis set so that, with an 800 eV cut-off, the barrier heights are accurate to within 1 meV. This stability is obtained despite the total crystal energy still evolving significantly as a function of the basis set cut-off; from y10087 eV at 500 eV cut-off to y10101 eV at 800 eV cut-off for PA and from y16733 eV at 500 eV cut-off to 16753 eV at 600 eV cut-off. For both systems, the calculated potentials are three-fold and slightly underestimated. Nevertheless, the V3 potentials relate accurately to the experimental observables: for PA, the tunnel frequency is within a factor of 2, for AN, it is within a factor of 3, the libration frequencies are within 2 meV and, signifi-

59

cantly, the equilibrium orientations are correct to within 18. The virtue of the DFT calculation is that it confirms the crystal structure and SPM rotational trajectory here for PA and AN. On the other hand, the computational cost Ž; 200 h CPU time per rotational potential. is a limiting factor. A computationally more efficient solution is desirable and insight into the deficiencies of the standard calculation is sought from semi-empirical quantum chemistry calculations.

5.3.2. Semi-empirical methods Up to 200 atoms can be handled by the semi-empirical methods embodied in MOPAC6 w38,39x. With the possibility of calculating an ensemble of up to nine PA molecules, two solutions can be envisaged. The first solution is to calculate the rotational potential for the central methyl group in the biggest cluster possible, knowing that these will be smaller than the ˚ radius clusters used in the standard calcula20 A tion, where the radius ensures convergence of the variation in the interactions with rotation angle. The second solution is to calculate charges for a central molecule taking into account the neighbouring molecules connected to the probe by hydrogen bonds, then using these charges instead of single molecule charges in the standard calculation. 5.3.2.1. Semi-empirical cluster calculations. MOPAC AM1 calculations on all systems give intramolecular rotational potentials which are stronger than the

Table 3 Results of DFT calculations for PA and AN as a function of the basis set cut-off. The potential parameters relate to the expression Ž V3 r2.Ž1 q cosŽ3Ž w q w 0 ... 500 eV cut-off

600 eV cut-off

800 eV cut-off

Observed

PA Potential Ž V3 , w 0 . ŽmeV, 8. Tunnel splitting ŽmeV. Libration frequency ŽmeV. Equilibrium orientation Ž8.

38.0, y14.6 4.53 13.4 74.6

36.8, y15.3 5.26 13.2 75.3

36.4, y15.0 5.42 13.1 75.0

3.10 13.8 a 75.8

AN Potential Ž V3 , w 0 . ŽmeV, 8. Tunnel splitting ŽmeV. Libration frequency ŽmeV. Equilibrium orientation Ž8.

55.3, 6.60 0.729 16.5 53.4

53.7, 6.44 0.850 16.3 53.6

a

This frequency does not correspond a unique spectral peak but is the average of two frequencies Žsee Section 3.1..

0.3 18.0 54.7

60

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

GAMESS potentials ŽFig. 7a.. Furthermore, the minimum of the MOPAC AM1 intramolecular potential occurs at 08 in each case, whereas in the GAMESS calculations, the minima occur at 1008 for PA and 308 for AN. The MOPAC potentials are less sensitive than the GAMESS potentials to the differences in molecular geometry between the samples. The variation in ground state electronic energy as a function of the methyl group orientation can be calculated with the AM1 method for several molecules. Fig. 7b shows three such curves for one, five and nine molecules of PA. A cluster of nine molecules

can be regarded as a reasonable approximation to the total rotational potential, only Coulomb interactions to more distant molecules have to be added. However, although the potential has a reasonable amplitude Ž V3 s 49 meV., the equilibrium orientation remains incorrect, being dominated by the single molecule potential which is strong and incorrect. For this reason, we reject the possibility of using MOPAC to calculate the rotational potential from large clusters of molecules. The PM3 Hamiltonian in MOPAC gives weaker intramolecular rotational potentials Žsee Fig. 7a for PA. but they are again relatively insensi-

Fig. 7. MOPAC derived rotational potentials Ža. for a single molecule: PA-AM1 Hamiltonian Ždashed curve, solid circles., AN-AM1 Hamiltonian Žsolid curve, triangles. and PA-PM3 Hamiltonian Ždashed curve, open circles. and Žb. for larger clusters of PA: one molecule Ždashed curve, circles., five molecules Ždashed curve, squares. and nine molecules Žsolid curve, triangles..

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

tive to the molecular geometry and have the minimum at 08. 5.3.2.2. Semi-empirical charges Õs. ab initio charges. Since the Coulomb rotational potential has been targeted as a likely source of error, MOPAC atomic charges for single molecules must first be compared with the GAMESS charges ŽFig. 8.. Despite the fact that the MOPAC charges used here are derived from the calculated electron density ŽAM1 Hamiltonian, Coulson charges., whereas the GAMESS charges are derived to reproduce the electrostatic potential around the molecule, the correspondence between the two

61

methods in terms of the variation of charge from atom-to-atom is reasonable. The MOPAC charges are significantly underestimated but, with an overall scale factor of 2.32 derived from a least squares fit, the root-mean-square-difference between GAMESS and MOPAC charges per atom can be reduced from 0.0420 to 0.0239. Fig. 9 shows the corresponding Coulomb rotational potential for each system, they have approximately the same phase as those derived in the GAMESS calculation. The amplitude of these Coulomb potentials is underestimated by factors of 3 and 1.5, respectively, and applying the charge scale factor of Ž2.32 2 s. 5.38 results in strongly overesti-

Fig. 8. Atomic charges for Ža. PA and Žb. AN. The dashed curves with circles show ESPD GAMESS charges while MOPAC AM1 charges are shown by a dashed curve with crosses for a single molecule and by a solid curve when the effect of the hydrogen bonds is incorporated. Atoms 1–8 are carbon and 9 is nitrogen in both cases. For PA 10–11 are oxygen and 12–20 are hydrogen, and for AN 10 is oxygen and 11–19 are hydrogen. Arrows indicate oxygen and hydrogen atoms involved in hydrogen bonds.

62

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

Fig. 9. Coulomb rotational potentials for Ža. PA and Žb. AN. The dashed curves with circles are the result of the GAMESS calculations, MOPAC AM1 calculations on single molecules are shown by dashed curves with triangles and on hydrogen bond clusters by the solid curves.

mated potentials. Charges derived with the PM3 Hamiltonian are in poorer agreement than the AM1 charges, notably for the carbon atom in the benzene ring to which the peptide group is bonded in PA which is assigned a positive, instead of a negative charge. 5.3.2.3. Charge modifications due to hydrogen bonding and transforming the Coulomb rotational potential. MOPAC AM1 is therefore used to investigate the effect of intermolecular hydrogen bonds on the charge distribution in a single molecule and the consequences for the Coulomb potential. The

MOPAC AM1 calculations are repeated for the smallest clusters of molecules which include all the hydrogen bonds for a central probe molecule ŽFig. 5.. The Coulson charges calculated in this way for the probe molecule are then used instead of single molecule charges to calculate the intermolecular Coulomb interactions. Oxygen and hydrogen atoms in hydrogen bonds experience the biggest charge modifications. For example, in PA, the hydrogencharge in the peptide group changes from 0.23e to 0.25e, the oxygen-charge from y0.32e to y0.36e Ž; 10%. while for the hydroxy group, in which the oxygen forms two hydrogen bonds, the respective

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

63

Fig. 10. Coulomb rotational potentials as a function of polarisation p for Ža. PA and Žb. AN. The smallest amplitude curves correpond to no polarisation. The amplitudes grow monotonically with polarisation, each curve correponds to 0.1e being subtracted successively from the peptide oxygen atom.

changes are 0.23e to 0.28e and y0.26e to y0.32e Ž; 20%.. Fig. 9 shows the Coulomb potentials calculated with the hydrogen bond cluster-derived charges. The results are interesting showing the desired tendencies. For PA the new charges give rise to a 78 shift in the Coulomb potential to lower angle, with a small decrease in amplitude. 9.81 2 ™

Ž 1 q cos Ž 3 Ž w y 42.7 . . . 9.55 2

Ž 1 q cos Ž 3 Ž w y 36.0 . . . .

Ž 10 .

For AN, the Coulomb potential increases in amplitude by ; 25%. 12.0 2 ™

Ž 1 q cos Ž 3 Ž w y 4.30 . . . 15.1 2

Ž 1 q cos Ž 3 Ž w y 1.12 . . . .

Ž 11 .

Incorporating the modified charges gives a stronger effect for AN than for PA, despite the fact that even stronger charge deformations are encountered in PA in the hydrogen bonds formed by the

64

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

hydroxy groups, because the methyl group in PA is located between three hydroxy oxygen atoms which are symmetrically displaced about the rotor. Total potentials have not been shown for MOPAC AM1 calculations since the charges are underestimated and the intramolecular potentials are considered to be too strong. One solution is to ‘map’ the transformation of Coulomb rotational potentials, Eqs. Ž10. and Ž11., onto the GAMESS potentials. For PA, the transformation ŽEq. Ž10.. is essentially the change in phase, D w s y78. When combined with the VDW and intramolecular potentials the final phase shift is smaller, D w s y48, and the total potential increases in amplitude by 2 meV. For AN, the transformation ŽEq. Ž11.. is dominated by the increase in amplitude, the transformed total potential increases in magnitude by 5.2 meV Ž56.2 meV ™ 61.4 meV., the corresponding transformed observables are 0.693 meV for the tunnel frequency and 15.4 meV for the libration frequency. 5.3.2.4. Applying the semi-empirical charge transformation to the ab initio charges. An alternative solution is to apply the charge transformation induced by the intermolecular hydrogen bonds, as calculated by MOPAC AM1, to the GAMESS charges. In this way, charge rescaling is avoided. The charge transformation is taken as absolute changes in charges rather than percentage changes as the latter leads to a bigger non-zero molecular charge which must be corrected. Since the charge differences are small in MOPAC due to the underestimation of the charges, the charge transformation can be applied progressively to the GAMESS charges. Fig. 10 shows the results for both samples, increasing polarisation, as defined by the charge subtracted from the oxygen atom in the peptide group, increases the amplitude of the Coulomb potential. The amplitude changes are indeed dominated by the polarisation of the hydrogen bond oxygen and hydrogen atoms, that is, where the biggest charge modifications occur. However, the phase change of the potential in PA ŽFig. 10a. arises only when the minor changes in all other atoms are included which indicates the sensitivity of the Coulomb potential to the precise distribution of the charges. This sensitivity is the reason why the potential amplitude increases in PA, whereas, when the potential is transformed di-

Fig. 11. Tunnel frequency Žsolid curve. and libration frequency Ždashed curve. as a function of polarisation for AN. Experimental values are indicated by the horizontal lines, solid for the tunnel frequency, dashed for the libration frequency.

rectly ŽEq. Ž10.., there is a slight decrease. The evolution of the total rotational potentials with polarisation are simple to derive from Figs. 6 and 9. For AN, for which the transformed Coulomb potentials in Figs. 9 and 10 are consistent, the tunnel splittings and libration frequencies as a function of polarisation are shown in Fig. 11. The first point, polarisation 0.1e, already corresponds to a slightly stronger polarisation than is predicted by MOPAC. A polarisation of 0.2e would give the correct tunnel splitting, 0.3e the correct libration frequency! Fig. 11 has the virtue of demonstrating how the methyl group probes the charge distribution in neighbouring hydrogen bonds and may, for example, be used to test experimentally the ionicity of hydrogen bonds proposed for n-methyl acetamide w24x.

6. Conclusion The results of the standard calculation on this class of molecules, namely those containing the peptide group, are poor. Methyl–methyl coupling is not responsible for the disagreement as coupling potentials have been calculated to be negligibly small. Rotation–translation coupling has been considered but the systems here are not likely candidates, there being no evidence of such motion in the thermal parameters from the diffraction measurements. Thus, the error in the standard calculation is not attributed to the simplicity of the SPM model or the crystal structures and this is confirmed by the excellent

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

rotational potentials obtained from the CETEP calculations. That the otherwise successful standard calculation fails where the DFT calculation succeeds, this due to the omission of hydrogen bonds in the former, indicates that the methyl group is a sensitive probe of intermolecular hydrogen bonds. Proven DFT calculations can therefore be used to investigate the anomalous properties of hydrogen bonds in polypeptides. Polarisable atoms like oxygen and nitrogen have been encountered in a number of systems Žacetic acid, acetophenone, diacetyl., the only notable problem has arisen when calculating the equilibrium orientation of the methyl group in nitromethane w19x. Since the calculated rotational potentials are essentially determined by the VDW and Coulomb terms, and VDW parameters appear to be satisfactory for a large number of systems, the latter have been targeted as the most likely source of error. MOPAC cluster calculations reveal that the charges of oxygen atoms involved in theses hydrogen bonds are underestimated by 10–20%. Including this polarisation effect from MOPAC in the GAMESS Coulomb rotational potentials shows the required response, but the effect is small, due to the ry1 dependence of the Coulomb interactions, and is sensitive to the details of the charge distribution. The amplitude of the rotational potential for AN calculated in this way Ž61.4 meV. is significantly bigger than that determined by CETEP Ž53.7 meV., the spectroscopic observables from each calculation are comparable, but the equilibrium orientation of the methyl group remains a problem in the modified standard calculation. Applying a charge transformation through a variable polarisation, rather than a transformation directly on the Coulomb rotational potential, gives similar effects, although they can be bigger if the polarisation is allowed to exceed that predicted by the MOPAC cluster calculations. While the results of the cluster calculations are qualitatively interesting, particularly for AN, ultimately they do not allow a reliable improvement of the standard calculation. The Coulomb correction depends sensitively on the charge distribution and it may be that ab initio cluster calculations would give satisfactory results. However, such calculations are intractable with large basis sets for several molecules each with ; 20 atoms.

65

Acetamide has been referred to in this article as the original molecule of interest containing a peptide group. All of the calculations presented in this paper have been applied to acetamide and the results are highly inaccurate, the amplitude of the rotational potential is overestimated by a factor ranging from 2.5 to 10, the upper limit being given by the DFT calculation which took over 1000 h of CPU time. Such results lead us to believe that the crystal structure may be inaccurate w40x. We plan to pursue this idea with new diffraction measurements.

Acknowledgements The authors are grateful to Eric Sandre, ´ Mike Payne and Julian Gale for their assistance with initial CASTEPrCETEP calculations.

References w1x W. Press, Single-particle rotations in molecular crystals, Springer Tracts in Modern Physics 92 Ž1981. . w2x M. Prager, A. Heidemann, Chem. Rev. 97 Ž1997. 2933. w3x K.J. Abed, S. Clough, Chem. Phys. Lett. 142 Ž1987. 209. w4x K.J. Abed, S. Clough, A.J. Horsewill, M.A. Mohammed, Chem. Phys. Lett. 147 Ž1988. 624. w5x B. Black, G. Majer, A. Pines, Chem. Phys. Lett. 201 Ž1993. 550. w6x J. Meinnel, W. Hausler, M. Mani, M. Tazi, M. Nusimovici, ¨ M. Sanquer, B. Wyncke, A. Heidemann, C.J. Carlile, J. Tomkinson, B. Hennion, Physica B 180–181 Ž1992. 711. w7x M. Prager, N. Wakabayashi, M. Monkenbusch, Z. Phys. B 94 Ž1994. 69. w8x M. Prager, M. Monkenbusch, R.M. Ibberson, W.I.F. David, D. Cavagnat, J. Chem. Phys. 98 Ž1993. 5653. w9x M. Prager, W.I.F. David, R.M. Ibberson, J. Chem. Phys. 95 Ž1991. 2473. w10x A.M. Alsanoosi, A.J. Horsewill, S. Clough, J. Phys.: Condens. Matter 1 Ž1989. 643. w11x A.M. Alsanoosi, A.J. Horsewill, Chem. Phys. 160 Ž1992. 25. w12x B.M. Rice, S.F. Trevino, J. Chem. Phys. 94 Ž1991. 7478. w13x M.R. Johnson, M. Neumann, B. Nicolai, P. Smith, G.J. Kearley, Chem. Phys. 215 Ž1997. 343. w14x M.A. Neumann, M.R. Johnson, A. Aibout, A.J. Horsewill, Chem. Phys. 229 Ž1998. 245. w15x M.A. Neumann, M.R. Johnson, P. Radaelli, H.P. Trommsdorff, S.F. Parker, J. Chem. Phys. 110 Ž1999. 516. w16x P. Schiebel, G.J. Kearley, M.R. Johnson, J. Chem. Phys. 108 Ž1998. 2375. w17x B. Nicolai, G.J. Kearley, M.R. Johnson, F. Fillaux, E. Suard, J. Chem. Phys. 109 Ž1998. 9062.

66

M.R. Johnson et al.r Chemical Physics 244 (1999) 49–66

w18x B. Nicolai, E. Kaiser, G.J. Kearley, F. Fillaux, A. Cousson, W. Paulus, Chem. Phys. 226 Ž1997. 1. w19x M. Neumann, M.R. Johnson, J. Chem. Phys. 107 Ž1997. 1725. w20x C.C. Wilson, Chem. Phys. Lett. 280 Ž1997. 531. w21x S.W. Johnson, J. Eckert, M. Barthes, R.K. McMullan, M. Muller, J. Phys. Chem. 99 Ž1995. 16253. w22x R.L. Hayward, H.D. Middendorf, U. Wanderlingh, J.C. Smith, J. Chem. Phys. 102 Ž1995. 5525. w23x M. Barthes, in: S. Cusack, H. Buttner, M. Ferrand, P. ¨ Langan, P. Timmins ŽEds.., Proceedings of a Workshop on Inelastic and Quasielastic Neutron Scattering in Biology, Adenine Press, 1997, p. 35 and references therein. w24x F. Fillaux, J.P. Fontaine, M.-H. Baron, G.J. Kearley, J. Tomkinson, Chem. Phys. 176 Ž1993. 249. w25x M.F. Guest, J. Kendrick, J.H. van Lenthe, K. Schoeffel, P. Sherwood, GAMESS: Generalised atomic and molecular electronic structure system, Daresbury Lab. w26x A.K. Rappe, C.J. Casewit, K.S. Colwell, W.A. Goddard, W.M. Skiff, J. Am. Chem. Soc. 114 Ž1992. 10024. w27x S. Clough, A. Heidemann, A.J. Horsewill, M.N.J. Paley, Z. Phys. B: Condens. Matter 55 Ž1984. 1. w28x W. Hausler, A. Huller, Z. Phys. B: Condens. Matter 59 ¨ ¨ Ž1985. 177.

w29x Neutronenstreuexperiments am FRJ-2 in Julich, available ¨ from the Secretary ŽNeutron Scattering., IFF, Forchungszentrum Julich, D-52425 Julich. ¨ ¨ w30x M. Matuschek, A. Huller, Can. J. Chem. 66 Ž1988. 495. ¨ w31x A. Wurger, A. Huller, Z. Phys. B 78 Ž1990. 479. ¨ ¨ w32x A.C. Hewson, J. Phys. C15, 1985, pp. 3841 and 3855. w33x K.J. Tilli, B. Alefeld, Mol. Phys. 36 Ž1978. 287. w34x J. Baudry, R.L. Hayward, H.D. Middendorf, J.C. Smith, in: Cusack, S., Buttner, H., Ferrand, M., Langan, P., Timmins, ¨ P. ŽEds.., Proceedings of a Workshop on Inelastic and Quasielastic Neutron Scattering in Biology, Adenine Press, 1997, p. 49. w35x http:rrwww.ill.frrYellowBookr. w36x CRYSTAL—Electronic Structure of Periodic Systems Žsee http:rrwserv1.dl.ac.ukrTCSrSoftwarerCRYSTAL.. w37x M.C. Payne, M.P. Teter, D.C. Allan, T.A. Arias, J.D. Joannopoulos, Rev. Mod. Phys. 64 Ž1992. 1045. w38x MOPAC: A Generalised Molecular Orbital Package, J.J.P. Stewart, Frank J. Seiler Research Lab., U.S. Air Force Academy, CO. w39x W. Thiel, in: I. Progogine, A. Rice ŽEds.., Advances in Chemical Physics, Volume XCIII, Wiley, 1996, p. 703. w40x G.A. Jeffrey, J.R. Ruble, R.K. McMullan, D.J. Defrees, J.S. Binkley, J.A. Pople, Acta Cryst. B 36 Ž1980. 2292.