Tunnelling paths for proton transfer reactions in solution. A discussion over a simplified model

THEO CHEM Journal of Molecular Structure (Theochem) 371 (1996) 161-169

Tunnelling paths for proton transfer reactions in solution. A discussion over a simplified model’ Ricard Gelabert, Miquel Moreno, JosC M. Lluch” Departament

de Quimica, Universitat Autbnoma de Barcelona,

08193 Bellaterra

Barcelona Spain

Received 28 December 1995; accepted 21 March 1996

Abstract Tunnelling of proton transfer reactions in solution is analysed through a molecular dynamics study based on a classical trajectories plus semiclassical tunnelling methodology. The potential energy surface is modelled as a proton moving on a symmetric double well interacting with a couple of dipoles, the whole system being kept collinear. Three different tunnelling paths are considered: (a) only motion of the proton through a symmetric potential, (b) only motion of the proton through an asymmetric potential, and (c) motion of both the proton and the solvent through a symmetric potential. Strictly speaking, motions of type (a) are never found. By comparing the tunnelling rates obtained through paths (b) and (c) it is found that at low energies (c) is dominant whereas (b) rate constants are slightly greater at high energies. It is predicted that for a more realistic model, with more solvent molecules included, paths of type (c) will dramatically lower their probability so that tunnelling will be dominated by paths of type (b). Finally, analysis of the (b) tunnelling paths reveals that the total tunnelling rate is dominated by the few paths which have the smallest asymmetry. This result agrees with the extended idea that proton transfer rates in solution depend on the difficulty in reaching a fluctuation of the solvent which yields a symmetric energy profile for the motion of the proton. Keywords:

Proton transfer reactions; Tunnelling

Paths; Solvent effect; Semiclassical

1. Introduction Because

proton

transfer

fulfils

a central

role

in

many processes,

understanding and modelling of proton-transfer rates is an important challenge [l-8] (the entire issue of Ref. [6] is devoted to proton-transfer reactions). Theoretical understanding of these reactions is difficult because of three main reasons: in the first place, protons are very light particles and consequently can only be treated properly via quantum mechanics; second, the proton motion is, or could be, strongly coupled to the motion of the rest of the

* Corresponding author. ’ Dedicated to Professor birthday.

J. Bertrkn on the occasion

of his 65th

atoms in the molecule; and third, if the reaction takes place in a polar solvent then there is also a strong intermolecular coupling between the proton’s and the solvent molecules’ motions. There has been considerable interest recently in the microscopic study of such reactions in solution[9-221. Most of the theoretical methodologies used to this goal have been based on molecular dynamics methods. The most prominent quantum effect which one can expect from proton-transfer reactions is the so-called tunnelling effect. Tunnelling stems from the propagation through a potential energy barrier of a wave packet linked to a given energy state. The semiclassical interpretation of this propagation is that the energy is not to vary along the tunnelling path. Bearing in mind that the proton-transfer system is surrounded by

0166-1280/96/$15.00 0 1996 Elsevier Science B.V. All rights reserved PII SO166-1280(96)04635-O

tunnelling model

162

R. Gelabert et al.lJournal of Molecular Structure (Theochem) 371 (1996) 161-169

molecules of solvent, the possible tunnelling paths could be classified into four different groups, according to the kind of motion which takes place along the path: (a) There is only motion of the proton through a symmetric potential (keeping the solvent structure frozen) (b) Only proton motion occurs through an asymmetric potential (also keeping the solvent structure frozen) (c) Motion of both proton and solvent molecules takes place through a symmetric potential (d) Motion includes both proton and solvent molecules, this time through an asymmetric potential It has to be noted that proton motion implies motion of the proton coupled to the internal molecular vibrations. The analysis of the kind of motions described above shows that only for the one detailed in (a) is a solvent reorganization previous to the proton-transfer compulsory. In the other cases this previous solvent reorganization is possible but not necessary. As it is well known, tunnelling processes are favoured when the energy profile through the barrier is symmetric [23]. The tunnelling amplitude falls down dramatically even when this symmetry is only slightly distorted [24]. Furthermore, the more mass the tunnelling system has, the lower the tunnelling amplitude becomes. Taking these two points into account, the tunnelling paths of type (d) above can safely be disregarded. On the other hand, Hynes and co-workers [ 171 have developed a theory for proton transfer reactions in solution which leads to a thermal rate constant which depends on, among other factors, the activation free energy for the solvent (AGi) to reach a symmetric

Fig. 1. Schematic

representation

situation, yielding an energy profile of type a in our terminology. The aim of this paper is to discuss what is the relative weight of the different contributions, derived from each of the different tunnelling paths cited above, to the rate constant at a given energy, and in addition to that, to test the conclusion previously obtained by Hynes and co-workers. To achieve this, we have considered a very simplified model of proton transfer in solution.

2. Calculational details

2.1. Reaction model The proton-transfer has been modelled as a symmetric intramolecular monodimensional protontransfer in the gas phase. The solvent has been modelled as a pair of rigid dipoles which reproduce the dipolar momentum of water. For the sake of simplicity, the whole system (proton and two dipoles) was restricted to move in a collinear fashion. A scheme of this model is shown in Fig. 1. In this model the reaction consists of the transfer of the proton from the left side of the potential to the right, passing through zero. Any other atoms of the reactive system have not been considered, so that their coupling with the proton motion is eliminated. In this way, this paper can focus exclusively on the solvent effects. Each dipole is placed each side of the x-axis and their motion is confined in such a way that they cannot go beyond a certain point XL. This restriction has a two-fold purpose: first, it serves to prevent the collapse of each dipole on the proton (a foreseeable consequence of the Coulombic attraction between them),

of the model used to study proton transfer in solution in this paper.

163

R. Gelaberi et al.Uournal of Molecular Structure (Theochem) 371 (1996) 161-169

and second, it serves to reproduce to some extent the “volume” which the solute has in reality but of which it has been deprived in this model. The explicit formula for the potential energy surface for a model such as this is the following one:

+ VD-&D,7XD2) +

C

i=x,

vrep(xi)

(1)

*DlrXD2

Where: V,(x,)=bx$

--ax;

(2)

Here, a and b are parameters taken from Ref. [25] to reproduce the energy profile of the proton transfer for the hydrogenoxalate anion. xH stands for the coordinate of the proton and xD1 and xD2 stand for the coordinates of both dipoles. VH_D and VD_D represent the Coulombic interaction between a dipole and a point charge (proton), and between a pair of dipoles. As for v rep. . Vrep(Xi)

=

kp

[(&y2+(gg2] C3)

Vrep is provided to prevent both dipoles from collapsing on the proton. The value of XL represents the

abscissa value beyond which no dipole would pass. A value of xL = 3.50 A has been used throughout the calculations. The numerical value for krep has been taken from Ref. [26]. Having a total of three independent variables, it is not possible to represent this potential on a twodimensional graph. Fig. 2 represents a cut from the whole hypersurface with dipole 2 fixed at infinite distance from the origin. This plot represents the potential energy surface for the proton transfer restricted to the neighbourhood of the reactant’s valley.

2.2. Tunnelling model The semiclassical tunnelling model combined with the use of quasiclassical trajectories as described by Makri and Miller [27] has been adopted. This method seems to be especially adequate to extract useful chemical information in order to understand the features that govern proton transfers in chemical and biochemical systems [28]. There are other semiclassical theories which describe tunnelling, such as the instanton mode1[29], but they require the use of classical sophisticated complicated more and mechanics. The main features of the model are summarized

-6.0

Fig. 2. A selected cut of the potential energy hypersurface

for the proton transfer, in which dipole 2 is held at infinite distance

164

R. Gelabert et aLlJournal of Molecular Structure (Theochem) 371 (1996) 161-169

here. The net probability that the system has tunnelled at time t, P&t), is given semiclassically by P&t)

= 1 h(t - t,) eeao n

(4)

where h is Heaviside’s step function: h(t) = 1,

if f > 0,

h(F) =O,

if (f
(5)

t, are the times at which the classical trajectories reach a turning point, a! is a parameter which depends on the symmetry of the energy profile along the tunnelling path:

CY-1,

Symmetric profile

(Y=2,

Asymmetric profile

(6)

Finally, 6 is the classical action integral through the energy barrier, given by 8= iIm

(7) s barrier

where s < and s > are the initial and ending points of the tunnelling trajectories. It is worth noting that the classical action integral, as defined above, has to be calculated by using mass-weighted Cartesian coordinates. The tunnelling rate is directly related to the time averaged net probability, Poet, over the initial phase of the vibrational motion in the well: ksym= j-$&)

k asym =

f-/P&t))

symmetric tunnelling path

asymmetric tunnelling path

(8)

Tunnelling rate constants for tunnelling paths of type a, will be evaluated by means of ksym, as will those of type c, while those of type b will be computed by means of kasym, Here, a, b and c refer to the different types of tunnelling path as classified in the introduction.

3. Results and discussion As previously stated, the weight on the tunnelling rate constant of the different contributions from each different tunnelling path a through c has been

analysed. In order to do so a methodology of classical trajectories has been used. The initial conditions for the trajectories have been chosen from the microcanonical ensemble for each given energy level, and in addition to that, the initial position has been forced to lie in the left well (i.e. xH is negative in sign). The Verlet velocity algorithm [30] has been used throughout the calculations to integrate Newton’s equations of motion. The integration step was taken to be 0.2 fs, and the trajectories obtained preserved the total energy within a 5 0.05% error. In addition to that, the ratio of the root of mean square deviations (henceforth, RMSD) of total energy to the RMSD of kinetic energy for each given trajectory was computed. This ratio never exceeded the value of 0.001 for the trajectories under study, which should be compared to a maximum threshold value of 0.1, as suggested in Ref. [31]. All trajectories were integrated for a total time of 200.0 fs and a total of 1000 trajectories were calculated at each energy. A total of ten different total energy situations were computed. The results which appear in this paper for the magnitudes at each energy are actually mean values over the individual results for each of the corresponding 1000 trajectories. The monodimensional tunnelling model presented in the previous section needs to be generalized in order to be applied to our multidimensional case. The first point to address is the definition of a turning point, that is, the points suitable for starting the tunnelling trajectory have to be defined. It has been assumed, in an analog way to what happens in a monodimensional system, that tunnelling events are likely to happen when the trajectory experiences a classical turning point for the proton (henceforth, CTP), i.e. a point in which the speed component for the tunnelling particle (the proton) changes its sign from positive to negative. When one of these points is found, and in order to quantify the magnitude of the tunnelling effect which it can produce, it is necessary to define what the tunnelling trajectory will be like in the multidimensional energy surface. Three cases have been considered, corresponding to paths of type a, b and c. In the first two situations, the tunnelling trajectory under the barrier is obtained by keeping all the dipoles fixed at the coordinates they had when the CTP was found, and then moving the proton towards positive XH

165

R. Gelabert et al./Journal of Molecular Structure (Theochem) 371 (1996) 161-l 69 0.50

0.45

0.40

I

0.35

I

0.30 5i 3 x

0.25

g i! c1

0.20 1

I

0.15

I

0.10

0.05 -

0.00

i

I

-I

_ 0

20

40

60

80

100

120

140

160

180

200

t (fs) Fig. 3. The net tunnelling

probability

P&(t) for a sample trajectory of 200 fs duration.

values. Due to a certain configuration of the solvent

molecules at that given CTP, it is possible that the exit of this tunnelling trajectory does not exist. This is likely to happen when the structure of the solvent greatly stabilizes the reactant’s structure. As for the case in which there is motion of the solvent and the proton (type c), the exit of the tunnelling trajectory is obtained by changing the sign of xH, and then making the following set of changes: xDI(exit) = - xD2(entry) and x&exit) = - xDI(entry). Obviously, now the exit point always exists, no matter what the energy is. The tunnelling path is defined as the linear path in massweighted Cartesian coordinates which links both structures. Once the tunnelling path has been defined, 8 is evaluated according to Eq. (7). The tunnelling amplitude, P&t), is calculated as shown in Eq. (4), separately for each different tunnelling path considered. An example of what this step-function looks like can be seen in Fig. 3.

Obtaining tunnelling rate constants from net tunnelling amplitudes is straightforward by means of Eq. (8). It is worth mentioning that, as a working hypothesis, it has been assumed that all tunnelling takes place via tunnelling paths of types a, b or c, without mixing them. Hence, a value of k,(E) has been obtained from the values of P,,et coming from tunnelling paths of type a, averaged over 1000 trajectories at energy E. Needless to say, an equivalent procedure has been followed to calculate kb(E) and k,(E). Fig. 4 summarizes these results. At low energies, the tunnelling rate constants derived from a tunnelling path of type c, k,(E), are dominant. In part, this is because at these energies, some CTPs would have no possible exit if tunnelling had to occur via motion of the proton only. When global energy is increased, there is a significant increase in kb(E) values. As expected, k,(E) = 0 for all the energies studied, because no CTP was ever found which yielded a totally symmetric profile

166

R. Gelaberi et aLlJournal of Molecular Structure (Theochem) 371 (1996) 161-169 12.0

11.0

10.0

9.0

/

7.0

/ 6.0

/

/

/

/+

/

/

/

/ ,J

/

J /

5.0

4

4.0

0.0

1.0

2.0

3.0

4.0

5.0

6.0

7.0

8.0

9.0

10.0

Total Energy &cal/mol) Fig. 4. Values of k,,(E) (dashed line) and k,(E) (solid line) in s-‘. Given that k,(E) is zero for all energies represented in this graph.

for the motion of the proton with the dipoles frozen in their positions. At the higher energies studied, both k&3) and k,(E) have a similar value. However, given that type c motion implies motion of both solvent molecules and the proton, if the number of solvent molecules is increased the length of the hmnelling paths of type c will dramatically increase and the same will happen to 0. Consequently, the individual probability for each tunnelling event of type c will fall down drastically once the number of solvent molecules increases, and hence kb(E) is likely to be clearly higher than k,(E). According to this, all the results seem to indicate that tunnelling in solution occurs involving only motion of the proton. Hynes and co-workers [17], in their study of tunnelling reactions, conclude that tunnelling in solution occurs with motion only of the proton through a symmetric potential, i.e. a tunnelling motion of type a in our terminology, and that the rate depends on the

studied, its value cannot be

difficulty in reaching a fluctuation of the solvent which yields a symmetric energy profile for the motion of the proton. In the end, tunnelling rates depend on the value of free energy (AG$ needed in order to reach such fluctuations. In the present work the results so far seem to show that tunnelling in solution has to occur via motion of the proton only, but through an asymmetric profile. However, the values of k,,(E) which have been obtained might comprise some CTPs in which the asymmetry is so slight that a nearly-symmetric profile, that is, a nearly type-u situation, has been achieved. Some effort has been devoted to analyse this point. Tunnelling trajectories encountered during a trajectory in which only the proton moves and which keep the solvent structure frozen have, most likely, an asymmetric energy profile. However, in order to somehow quantify the asymmetry of such a profile, an “asymmetry parameter” 6 has been defined for

R. Gelabert et aLlJournal of Molecular Structure (Theochem) 371 (1996) 161-169

Fig. 5. Evaluation of the 6 asymmetry parameter (see text). s < and s > are the entry and exit points, respectively, x> is the value for which the potential under the barrier reaches its maximum.

each tunnelling path as: 6= I 2

(9) where xi is the value of xH for which the potential reaches a maximum, and w is the shortest of the two segments which link the maximum in energy along the tunnelling path and either the entry point or the exit point of the tunnelling trajectory (see Fig. 5). xD1 and xD2remain frozen along the tunnelling path when computing 6. Both parameters need to be evaluated at each new CTP. It is clearly seen that 6 = 0 only if the energy profile through the barrier is symmetric with respect to the maximum, and for asymmetric profiles the greater the asymmetry, the greater the value of 6. In this way, by evaluating 6 at each CTP, it is possible to classify that CTP with respect to its symmetry. This parameter has been evaluated over all the trajectories and a summary of the results is presented in Table 1,

167

of the tunnelling trajectory.

Table 2 and Table 3. Only three tables of results, corresponding to three selected total energies, are presented here for the sake of conciseness. As stated previously, strictly speaking the profile through the barrier will always be an asymmetric one. Consequently, individual tunnelling probabilities for each CTP have been computed using the following Table 1 Analysis of the symmetry properties of the classical turning points. (Total energy = 1.00 kcal mol.t; total number of effective CTF% = 8599) log ( < P ' 1

Number of effective CTPs in this group

-8.97 -8.84 -8.70 -8.56 -8.42 -8.28 -8.14 -8.00 -7.86 -7.72

52 365 876 1415 1553 1375 1443 1113 406 1

< 6 z

0.091 0.082 0.071 0.068 0.060 0.054 0.045 0.032 0.019 0.0044

168

R. Gelabert et al.IJournal of Molecular Structure (Theochem) 371 (19%) 161-169

Table 2 Analysis of the symmetry properties of the classical turning points. (Total energy - 4.00 kcal mol-‘; total number of effective CTPs = 8723) log ( < P ’ )

Number of effective CTPs in this group

-8.70 -8.26 -7.81 -7.36 -6.91 -6.46 -6.01 -5.56 -5.12 -4.67

226 618 1343 1664 2302 1516 684 303 66 1

< 6 >

0.070 0.064 0.050 0.042 0.032 0.027 0.021 0.012 0.0060 0.0019

formula, pi =e-2*1

(10)

which is the WKB probability for a tunnelling event taking place through an asymmetric barrier. CTPs have been grouped according to the individual probability they yield. The first column in Table 1, Table 2 and Table 3 shows the mean probability produced by a particular group of CTPs. CI’Ps which have no possible exit have been excluded from the tables. Column 2 shows how many CTPs fall into every interval. The third column holds the mean 6 value. Studying Table 1, representative of the situation at low energies, it is easy to see that all CTPs produce fairly low tunnelling probabilities. In addition to that, a comparison of the mean values of asymmetry of Table 3 Analysis of the symmetry properties of the classical turning points. (Total energy = 10.00 kcal mol-‘; total number of effective CI’Ps = 8090) log ( < P ’ )

Number of effective CI’Ps in this group

-8.45 -7.54 -6.63 -5.72 -4.81 -3.90 -2.99 -2.08 -1.17 -0.27

247 694 1052 1650 1831 1414 754 324 123 1

< 6 >

0.15 0.072 0.040 0.019 0.011 0.0069 0.0039 0.0016 0.00025 o.OOOOOOO95

Table 1 to those at higher energies (Table 2 and Table 3) shows that, in general, at lower energies tunnelling events occur through more asymmetric barriers. More interesting details can be derived from a study of the situation at higher energies. To begin with, it is easy to see that now the range of probabilities has increased considerably. In fact, the lower end remains more or less invariable, but remarkably, much higher probabilities can be obtained at high energies. This is not surprising, since having more energy means that the solvent has a wider range of configurations accessible and that fewer CI’Ps will have no exit. In addition to this, and with respect to the asymmetry parameter 6, there is a general trend towards lower values of 6 for the most high-probability yielding CTPS. In the high energy cases (Table 2 and Table 3) it is also worth mentioning that a few CTPs produce relatively high probabilities, in such a way that they alone account for the majority of the hmnelling probability obtained at that energy. Furthermore, these CTPs are also those with a lower asymmetry parameter. Then it can be said that all significant tunnelling in these trajectories comes from a few CTPs which have a fairly symmetric profile. In these cases, there will be no noticeable tunnelling probability until one of these particular CTPs appears, which could take time. In the end, the low number of CTPs which produce such probability is clearly related in some way to the AGj value proposed by Hynes and co-workers [ 171. The fact that, in our model, all significant tunnelling comes from CTPs with a nearly-symmetric profile is compatible with their conclusion that tunnelling in solution has to occur via a symmetric energy profile. In conclusion, our molecular dynamics simulation based on classical trajectories with the tunnelling incorporated through the WKB formalism has confirmed the extended idea that proton transfer tunnelling reactions in solution occur through tunnelling paths which involve only motion of the transferring proton. Moreover, tunnelling rates mainly depend on the probability of reaching a fluctuation of the solvent which yields an (almost) symmetric energy profile for the motion of the proton. In our very simple model, symmetric paths which involve motion of both the proton and the solvent molecules are found to yield

R. Gelabert et al.IJournal of Molecular Structure (Theochem) 371 (1996) 161-169

higher rate constants at low energies and only slightly lower ones at higher energies. However, in a more realistic model, with more solvent molecules included, the tunnelling probability of these paths is likely to fall down to almost zero.

Acknowledgements The authors, on behalf of the Reaction Dynamics Group of the UAB, are very pleased to participate in this Special Issue in honour of Professor Juan Bertran. All of us, who are, from a scientific point of view, “son” (J.M.Ll.), “grandson” (M.M.) and “greatgrandson” (R.G.) of Professor Bertran, wish to acknowledge his very intense and successful work at the Universitat Autbnoma de Barcelona. His enthusiasm, humour, joviality and agile mind have always accompanied those who have the privilege and the joy of working with him.

References [l] R.P. Bell, The Proton in Chemistry, 2nd Edn., Chapman and Hall, London, 1973. [2] E. Caldin and V. Gold (Ed?..), Proton Transfer Reactions, Wiley, New York, 1975. [3] V.K. Babamov and R.A. Marcus, .I. Chem. Phys., 74 (1981) 1790; V.K. Babamov, V. Lopez and R.A. Marcus, J. Chem. Phys., 78 (1983) 5621. [4] R. Stewart, The Proton: Applications to Organic Chemistry, Academic Press, Orlando, FL, 1985. [5] S. Scheiner, P. Redfern and M.M. Szczesniak, J. Phys. Chem., 89 (1985) 262. [6] Chem. Phys., 136 (1989) [entire issue]. [7] D. Li and G.A. Voth, J. Phys. Chem., 95 (1991) 10425. [8] J.L. Herek, S. Pedersen and A.H. Zewail, J. Chem. Phys., 97 (1992) 9046.

169

[9] W. Siebrand, T.A. Wildman and M.Z. Zgierski, J. Am. Chem. Sot., 106 (1984) 4083; 106 (1984) 4089. [lo] R. Meyer and R.R. Ernst, J. Chem. Phys., 86 (1987) 784. [ll] P.P. Schmidt, J. Phys. Chem., 90 (1986) 4945. [12] R.I. Cukier and M. Morillo, J. Chem. Phys., 91 (1989) 857; M. Morillo and RI. Cukier, J. Chem. Phys., 92 (1990) 4833; J. Chem. Phys., 98 (1993) 4548. [13] Y.R. Kim, J.T. Yardley and R.M. Hochstrasser, Chem. Phys., 136 (1989) 311. [14] A. Warshel and J.K. Hwang, Science, 246 (1989) 112; A. Warshel and Z.T. Chu, J. Chem. Phys., 93 (1990) 4003. [15] N. Shida, P.F. Barbara and J. AImlBf, J. Chem. Phys., 94 (1991) 3633. [16] A. Suarez and R. Silbey, J. Chem. Phys., 94 (1991) 4809. [17] D. Borgis and J.T. Hynes, J. Chem. Phys.. 94 (1991) 3619; Chem. Phys., 170 (1993) 315. [18] D. Laria, G. Ciccotti, M. Ferrario and R. Kapral, J. Chem. Phys., 97 (1992) 378. [19] T.N. Truong, J.A. McCammon, D.J. Kouri and D.K. Hoffman, J. Chem. Phys., 96 (1992) 8136; P. Bala, B. Lesyng and J.A. McCammon, Chem. Phys., 180 (1994) 271. [20] S.G. Christov, Chem. Phys., 168 (1992) 327. [21] H.J.C. Berendsen and J. Mavri, J. Phys. Chem., 97 (1993) 13464; J. Mavri, H.J.C. Berendsen and W.F. Van Gunsteren, J. Phys. Chem., 97 (1993) 13469. [22] Z. Smerdachina, V. Enchev and L. Lavtchieva, J. Phys. Chem., 98 (1994) 4218. [23] C. Cohen-Tannoudji, B. Diu and F. Laloe, Quantum Mechanics, Wiley, Paris, 1977. [24] R. Gelabert, M. Moreno and J.M. Lluch, J. Comput. Chem., 15 (1994) 125. [25] E. Bosch, M. Moreno, J.M. LIuch and J. Bertrdn, J. Chem. Phys., 93 (1990) 5685. [26] V. Perez, J.M. Lluch and J. Bertran, J. Comput. Chem., 13 (1992) 1057. [27] N. Makri and W.H. Miller, J. Chem. Phys., 91 (1989) 4026. [28] E. Bosch, M. Moreno and J.M. Lluch, Chem. Phys., 159 (1992) 99. [29] V.A. Benderskii, D.E. Makarov and C.A. Wight, Adv. Chem. Phys., 88 (1994) 1. [30] M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. [31] H.J.C. Berendsen and W.F. van Gunsteren, in A.J. Barnes et al. (Eds.), Molecular Liquids-Dynamics and Interactions, Reidel, Reading, MA, 1984, p. 475.