- Email: [email protected]
H-Bond in methanol: a molecular dynamics study
Journal of Molecular Structure, 250 (1991) 147-170 Elsevier Science Publishers B.V., Amsterdam
147
H-Bond in methanol: a molecular dynamics study* J. Alonso**‘, F.J. Bermejob, M. Garcia-Hern&rdez”*l, J.L. Martinez” and W.S. Howells” “Rutherford Appleton Laboratory, Chilton, Didcot, Oxon, OX1 1 OQX (UK) *Znstituto de Estructura de la Materia, Serrano 119, E-28006 Madrid (Spain) ‘Znstitut Laue Langeoin, 156X, F-38042 Grenoble-Cedex (France) (Received 15 October 1990)
Abstract A set of molecular dynamics (MD) simulations of methanol-d, at three temperatures in the liquid range (200,250 and 300 K) have been carried out. The equations of motion of 256 molecules interacting through a potential model due to Haughney et al. [J. Phys. Chem., 91 (1987) 49341 were solved using the velocity version of the Verlet algorithm. This rather large number of molecules was required for studying the behaviour of the system at momentum transfers as low as 0.25 A-‘. It was found that the system experiences long period fluctuations, and therefore very long MD runs (of the order of 100 ps) are necessary in order to obtain accurate statistical averages. Computed static properties are in good agreement with those reported by Haughney et al. and the neutron weighted g (r) and the static structure factor compare favourably with available neutron diffraction data. The study of time-dependent properties through centre-of-mass autocorrelation functions (VACF, F. (Q, t ) and F( Q, t ) ) and their memory functions reveals features unknown in simple liquids and very similar to those found in liquid water. A close agreement between centreof-mass single-particle autocorrelation functions and the translational part of QENS data is also observed. The dynamic structure factor for the centres of mass show distinctive side peaks in the same region of the (Q, w) plane where recently collective excitations have been studied using coherent neutron scattering thus establishing the presence of propagating short wavelength modes.
INTRODUCTION
The microscopic dynamical properties of hydrogen-bonded liquids have attracted a great deal of attention in recent times, with special emphasis on the investigation of single-particle and collective dynamics [ 11. Some of the recent efforts have been devoted to the experimental investigation of such phenomena by means of neutron spectroscopy [ 2,3]. The quan*Presented at the International Symposium on Hydrogen Bond Physics held at 11Ciocco, Barga, Italy, 11-14 September 1990. ‘Permanent address: Instituto de Estructura de la Materia, Serrano 119, E-28006 Madrid, Spain.
0022-2860/91/$03.50
0 1991 Elsevier Science Publishers B.V. All rights reserved.
148
tities measured in such experiments can be related to the single-particle (i.e. low frequency rotational and translational excitations) and collective (phononlike) motions. However, difficulties in both the experimental measurements and in the interpretation of the measured spectra which have to be done in terms of fully microscopic models hinder future progress unless some additional information from other sources is available. As an example, let us recall the fact that most of the available literature concerning experimental neutron scattering data is based upon some crude models to analyze the data, mostly based on jump or diffusion-like approximations [ 41. The atomistic computer simulation methods are nowadays firmly established as a very useful tool to complement the experimental observations as well as to provide information on those quantities not amenable to experimental investigation [ 41. The purpose of the present paper is to contrast the results arising from recent neutron scattering experiments [ 2,3] with those derived from molecular dynamics (MD) simulation for a relatively simple hydrogen-bonded liquid, methanol. The single-particle dynamics has been studied by means of time-of-flight spectroscopy [ 21 and the data has been analyzed on the basis of a recently proposed model [5] where the different association states are considered explicitly. Also, the collective dynamics has been investigated by means of tripleaxis spectroscopy [ 31. In this latter case the relatively low resolution achieved in the measurement precluded the determination of direct estimates of the lineshape of the underlying dynamical structure factor; and, therefore, one of the purposes of this paper is to compare the computer simulated structure factors with those obtained by means of deconvolution of the experimental intensities. The outline of the paper is as follows. Details about the MD simulation and the molecular model are presented in the next section. All results regarding hydrogen bonds are then summarized in the section HB network analysis, whilst the following section provides results on static properties. The penultimate section is devoted to the study of centre-of-mass dynamical correlation functions. The results obtained are finally discussed in the Conclusions section. THE MD SIMULATIONS
A computer program for simulating molecular liquids was written following the hints given by Allen and Tildesley [ 61. The program generates the trajectory of a system of molecules according to Newtonian classical mechanics ( NVE-P ensemble ) and using cubic periodic boundary conditions to simulate an infinite system. The molecules are described as a collection of interaction
149
sites whose Cartesian equations of motion are integrated using the velocity version of the Verlet algorithm [ 6,7]. Because of our intention to produce long simulation runs a rigid molecular model [8] was employed to keep the computations within reasonable limits. The RATTLE algorithm [6,9] was used to implement the holonomic constraints required to keep all intramolecular distances fixed with a tolerance of 10w6. The combination of these two algorithms made it possible to use a time step as large as lo-l4 s with good energy conservation (2 x lo-l4 over periods of 30 ps). A noticeable, though small, drift in the total energy was detected for long runs and we chose to rescale particle velocities every time step in order to ensure we were really sampling the microcanonical ensemble. All results presented here were obtained using the integration scheme mentioned above with a step of 10-14 s. Because of the excellent overall agreement with experimental data we have chosen model Hl of Haughney et al. [8] for our simulations. We have introduced just two slight modifications to the original model. On the one hand, we have substituted the mass of the hydrogen atoms by that of deuterium in order to simulate the fully deuterated compound. This choice has been motivated by the availability of recent neutron scattering data for the deuterated compound [ 2,3,10] with which to compare the results. On the other hand, instead of using the Ewald method for long range dipole-dipole interactions, a switch function [ 111 was introduced to smoothly switch off all interactions between molecular pairs whose centres-of-mass separation were greater than a given cut-off distance. The interaction energy between two molecules cyand /? then becomes
(1) where Rap is the distance between centres of mass for both molecules; S(X) is the unique fifth-order polynomial satisfying S(r)=1 =o
for x
(2)
for x>R&
and having continuous
first and second derivatives at the endpoints of the
TABLE 1 Int.errnolecularpotential parameters [81 Site, i
Oii/A
(d&)(K)
Carbon Oxygen Deuterium Deuterium
3.861 3.083 0.000 0.000
91.15 87.94 0.00 0.00
(hydroxyl ) (methyl)
qil
Ie I
0.297 - 0.700 - (qc+qo) 0.000
150 TABLE 2 Details of the MD runs. Values of internal energy and pressure do not include long range corrections Quantity
Unite
Density Elapsed time
g cme3 picoseconds
Temperature IJinte’
Kelvin kJ mol-’ kJ mol-’ atmospheres atmospheres 10m6 cm2 s-i ps-2 cm-’
GE Pressure PressureLRc Diffusion coefficient Einstein frequency, C$ Qll
Run 1 0.9430
177.00 202.3 f 0.3 - 39.88 -0.37 513.0 - 193.2 0.144 463.7 114.3
Run 2 0.9400
138.95 248.1 f 0.9 -37.78 -0.37 1285.6 - 192.0 0.54 443.6 111.8
Run 3 0.8850 100.10 297.7 f 1.0
-34.79 -0.35 1639.7 - 170.2 1.69 402.2 106.4
interval [ 111. In the present calculations we took RL= 12.25 A and RU= 12.65 A for all simulations. The function V”,, (rz ,r$ ) is the potential function referred to as model Hl by Haughney et al. [ 81 and has the form of a sum of sitesite interactions (3) where Uij( 1r z - r,$ 1) is a site-site interaction which only depends on the distance, rij, between site i in molecule LXand site i in molecule j? (4) Uij=
Eij
=
(Oii+Ojj)/2
(5)
Jtii
(6)
Cjj
All required potential parameters are summarized in Table 1. Three simulations of a 256 molecule system were carried out at temperatures near 200,250 and 300 K. Values of densities and run lengths are shown in Table 2. The initial configurations for each run were generated from final configurations of previous runs after an equilibration period no shorter than 50 ps. HB NETWORK ANALYSIS
Though the debate on hydrogen bond definition is still open, it seems there are two preferred criteria: energetic [ 121 and geometric [8]. Both of them sometimes being combined with a time analysis of the hydrogen bond [ 13-151. In this work, for the sake of simplicity, we shall use a purely geometric criteria
151
in order to define the hydrogen bond. Following Haughney et al. [8] we consider two molecules as bonded if r (0 *..O) 53.5 A, r(O**.D,) 12.6 A, and the angle ~D,O~~~0130”. Percentages of molecules in different bonding states are shown in Table 3. We observe a small fraction of molecules in a “fourfold” bonded state that could possibly be regarded as molecules involved in multiple collisions and not really as being hydrogen bonded. Fractions for the most traditional bonding states agree very well with those results previously reported for similar molecular models [ 8,12,15,16]. The average number of hydrogen bonds per molecule, nne, is also included in the Table. Standard deviations for average values were estimated by computing the statistical inefficiency (SI ) associated with each quantity [ 61, namely (7)
4 (g) ) = wN3bs)1’2 h!)
where s is the statistical inefficiency, Nabsis the number of observations, and a(g) is the standard deviation of g. We have found SI values as large as lo3 steps for some of the bonding fractions; this being a consequence of the long period oscillations (of the order of 30 ps) experienced by these molar fractions. The ultimate cause of these oscillations would be the slowness of relaxation phenomena in liquid methanol due both to the existence of the hydrogen bond network and to the low temperature range studied. Cross correlations between populations of different bonded states makes it possible for the SI for nnn to be about 30 steps for all temperatures. It is also interesting to notice the high TABLE 3 Hydrogen bond statistics. All fractions expressed in percentages and hydrogen bond lifetimes in picoseconds Run 1 0.052 + 0.012 6.60 f0.17 k
f4 nHB
GiBe rHBb ?HB=
Run 2 0.44 12.2
Run 3 kO.03 f 0.3
1.84 19.8
kO.05 ?I0.3
87.8 5.51 ff0.19 0.4 0.0027 + 0.0010 1.9881+ 0.0005
79.6 7.7 f 0.5 + 0.3 0.018 f 0.004 1.9471 z!z0.0009
69.8 8.5 f0.2 kO.5 0.056 + 0.010 1.8510 I?I0.0018
9.ld 20.4d 442.4d
2.6 5.5 74.7d
1.2 2.2 2l.Od
*From the time integral of the H(t) function. bFrom the time integral of the Cc (t ) function. “From the time integral of the C,(t) function. dExtrapolated value using a multiexponential fitting.
152
5
0
10
t/ps
15
20
0.9-
0.80.706s 0.5v- 0.40.3-
0
I, 0
, IO
, 20
I
, 30
I,
40
I,, 50 60 t IPS
70
60
Fig. 1.Hydrogen-bond autocorrelation functions at different temperatures. Note the different time periods covered by each plot.
153
correlation that exists between fluctuations of temperature and nnn about their average values which at 300 K is close to 37%. In order to characterize the hydrogen-bond dynamics we have used an approach based on bond autocorrelation functions. For that purpose we have computed three different functions. 1. H(t), the fraction of bonds unbroken after time t. 2. Cc(t), the autocorrelation function for a “continuous” definition of the hydrogen bond. 3. Ci (t), the autocorrelation function for an “intermittent” definition of the hydrogen bond. Note that the first two functions require the bond to be uninterrupted. If a bond breaks and reforms again, the reformed bond is treated as a new bond; while for the third function the new bond will be recognized as an old one. Details on the computation of these functions are given elsewhere [ 8,141. Hydrogen-bond lifetimes can be computed from rnn =
I
X(t) dt
(3)
0
where X(t) is one of the three above-mentioned autocorrelation functions. Results for the bond lifetimes are shown in Table 3. In those cases when the function X(t) has not yet decayed to zero a multiexponential fitting was integrated analytically in order to obtain the bond lifetime. HB autocorrelation functions are plotted for all three temperatures in Fig. 1. Differences between H( t ) and Cc (t ) arise from the greater weight given to short-lived bonds in the former function. Because most short-lived bonds are simply binary collisions or unstable bonds that break and reform in rapid sequence, we think function Cc(t) should be preferred. However, multiexponential fittings show that both functions have similar relaxation times at large t. A recent MD study of HB in liquid methanol [ 151 has shown that Cr( t) gives the best description of the HB dynamics because time stability of the bonds is implicitly required, and therefore C,(t) depends to a lesser degree than H(t) or Cc (t ) on the exact form in which the hydrogen bond is defined. . This function shows a new picture of liquid methanol in which the relaxation of the HB network has a time scale at least an order of magnitude bigger than previously thought, and bond lifetimes are almost in the nanosecond scale. STATIC PROPERTIES
Most of the thermodynamic properties we computed and monitored in real time (i.e. at the same time as the trajectories). The average values of temperature, intermolecular potential energy and pressure for the three runs are summarized in Table 2. The standard deviation of average temperatures was esti-
154
mated by computing the statistical inefficiency (see eqn. (7) ). Long range corrections (LRC) to the intermolecular energy and pressure were estimated by assuming total absence of spatial correlations beyond the cut-off distance [6], and are included in Table 2 for completeness. Because the standard deviation of average pressure and internal energy are smaller than the corresponding long range corrections they have not been included in the Table. Typical uncertainties in the average value of the internal energy can be roughly estimated as less than 1%. However, for the pressure averages these are far less reliable because the LRC amount to 20%, and even to 60% of the actual value are 200 K. Instantaneous values of dynamical variables were observed to experience long period (30-50 ps) oscillations which also manifest themselves in great statistical inefficiencies. Owing to its large amplitude, the influence of these oscillations on calculated average values can be noticeable unless the simulation is carried out long enough to average them. The only transport coefficient we determined was the diffusion coefficient. It was computed via the Einstein relation D=ilim$<
]ZZa(t+s)-&(a)
t-w
12)
(9)
where 2,(t) is the position of the centre of mass of molecule cy at time t; and the average is taken over molecules, (Y,and initial conditions, s. Alternatively the diffusion coefficient can also be computed from the time integral of the velocity autocorrelation function (IO) We have observed that both methods give indistinguishable results within the statistical uncertainties. However, it is necessary to obtain the velocity autocorrelation function up to times of about 10 ps in order to obtain a good agreement between both methods. Diffusion constants are given in Table 2 at the three temperatures. A comparison of the MD values with those obtained experimentally [ 21 show that the simulation value is about 25% smaller. Finally, we have calculated the Einstein frequency which is obtained from the second time derivative of the normalized velocity autocorrelation function at zero time
sL:yp
a2
+ >I (
$~t+~)-ti,(~))
(11)
t=o
where the calculated frequencies are also shown in Table 2 for reference. The static spatial correlations have been investigated via pair radial distri-
155
bution functions (RDFs). Centres-of-mass RDFs, g,, (r), for methanol at 200 and 300 K are plotted in Fig. 2. The main change with temperature is an enhancement of the structure at low temperature. The narrow peak at 3.5 A corresponds to hydrogen-bonded molecules, and then the first shell of non-hydrogen-bonded molecules appears as a second peak at 4.5 A. At longer distances g,, (r ) closely resembles the RDFs found for simple liquids. In Fig. 3 the six atom-atom RDFs involving only atoms C, 0 and Do are shown for the simulation at 200 K. These partial RDFs are very similar to those reported by Haughney et al. [8] at room temperature, except for the peak heights which have increased at 200 K. Again, the main change observed on reducing the temperature is the increase in peak heights. The existence of a
-E 1.0 UY
0.0
I--5
6
7
0
9
r/B,
(b)
0.6-
0,,/1
2,,,I,,,,,, 3
4
5
6
~,,,,,, 7
0
9
r/A
Fig. 2. (a) Centres-of-mass radial distribution functions for liquid methanol at 200 K (line) and 300 K (dote); (b) Neutron-weighted RDF at 200 K.
156
oc 23456789
r/h
Fig. 3. Set of selected atom-atom radial distribution functions for liquid methanol at 200 K.
hydrogen bond network in the liquid produces very strong correlations at short distances, as can be seen from the figure; the position of the first peak in the partial RDFs being determined by the HB geometry. To enable comparison with neutron diffraction data [lo] we have computed the neutron-weighted RDF, gN (r ), which is a combination of atom-atom RDFs given by
(12) where gij (r) is an atom-atom RDF and bi is the coherent scattering length of
157
nucleus type i. The function gN(r) (also shown in Fig. 2) exhibits a great number of maxima and shoulders corresponding to the maxima of the atomatom RDFs. It was known that this molecular model produces X-ray and neutron diffraction patterns close to those measured experimentally [ 81. The function gN(r) computed from MD has the same phase as its experimental counterpart, though some discrepancies concerning peak heights exist. CENTRE-OF-MASS
CORRELATION FUNCTIONS
The detailed microscopic information provided by MD simulations enables separate investigation of translational and rotational motions, and their contributions to the dynamic structure factor, S (&,a). In this work we have limited ourselves to study the molecular translation (centres-of-mass motion), leaving rotation for a forthcoming paper. However, the centre-of-mass correlation functions constitute a bench-mark to test the usual approximations employed in the experimental data analysis. Single-particle dynamics
The study of the centre-of-mass velocity autocorrelation function (VACF) has revealed an unsuspected complexity. From Fig. 4 we can see the influence of temperature on the VACF. As expected, the cage effect (i.e. the so-called backscattering) is weaker as temperature increases; the depth of the main minimum at 300 K being about one half the depth at 200 K. Thus indicating a lost of rigidity in the response of the liquid at frequencies of some terahertz. An interesting feature of the VACF at 200 K is the small local maximum at about
0.4
0.6
0.6
t/p5
Fig. 4. Normalized centre-of-mass velocity autocorrelation functions computed at 200 and 300 K.
158
0.25 ps which is absent at 300 K where only a shoulder is noticeable. At very short times the VACF is almost independent of the temperature for the three states studied here; this can be understood by making a short-time expansion of the normalized VACF [ 171 to obtain W(t)=+
L
+...
(13)
Values of Einstein frequencies are quoted in Table 2, where their small sensitivity to temperature variations can be noticed. Fourier transforms of VACFs, V/(U), are shown in Fig. 5. At 200 K a broad shoulder is apparent at frequencies close to the Einstein frequency; this being almost non-existent at 300 K. These temperature variations are enhanced in the spectra of the normalized force autocorrelation function (FACF) which is related to that of the VACF through the equation (14) It is very interesting to see the way in which the fine structure of the band centred at 200 cm-l, present at 200 and 250 K, disappears in the 300 K simulation making clear a qualitative change in the response of the liquid at intermediate frequencies. Assuming the dynamics of liquid methanol at frequencies of a few terahertz is dominated by the presence of an extensive hydrogen bond network, these results seem to indicate the existence of some kind of dynamical transition (possibly the percolation of the HB network) somewhere between 250 and 300 K. However, results of another recent MD simulation [ 151 suggest that simulated liquid methanol is, almost certainly, percolated at 200 K while it is not at 300 K. The memory function (MF) of the VACF, m (t ), was computed by solving the memory function equation [ 17,181 t
m(t)=C(t)-Sdsm(t-s)W(s)
(15)
0
where k(t) and 9(t) are the first and the negative of the second derivative of the normalized VACF, w(t), respectively. The memory function at 200 K is shown in Fig. 6 along with the VACF and the FACF. The most striking characteristic of the MF is its oscillatory behaviour inherited, without any doubt, from the FACF to which is closely related, at least at short times. Memory functions reported in the literature for simple liquids used to be simpler functions than the corresponding ACFs and very different from the one reported here. In order to investigate the quasielastic region of the (4,~) space we have
, 150
;, 200
004-
I
,
0
0
100
I 300
Wavenumber
I 200
Wavenumber
/ cm-’
I 400
/cm-’
I
, 300
,
, 350
,
, 400
560
600
0
0.01
: 0.03 i? 0.02
t; 0.05 i!i m 0.04
E I 0.06 L
0.07
0
“;“I
-
-
0
,
,
100
, 50
Fig. 5. Fourier transforms of normalized VACFs and FACFs for 200 and 300 K simulations.
0
2003iiT 002-
E
E 2 0.06* g! 0.05-
0.07-
0.04-
, 250 0.
,
0.02 -
, 100
0.
,
0.02-:
, 50
< 0.06> 0.04 -
0.06-
;
I
O.lO0.08-
4 O.lO-
0 2
K
0.12-
2
:
200
0.16 E o.l4-
L O.OS-
5 0.14;; 0.12-
0.16-
I,
I 100
4
, 150
,
,
.I,
, 200
300 Wavenumber
200
Wavenumber
,
,
,
,
/ cm-’
K
I
K
,
I,
, 350
500
spy
, 300
300
,
300
?? 400
/cm-’
, 250
I
,
6bO
, 400
160
ooooo00 CM
Velocity
I--.
CM
Force
+++++
CM
Velocity
ACF ACF memory
function
0.6
0.4
0.8
t I PS
Fig. 6. Normalized centre-of-mass velocity, force autocorrelation functions and memory function of the VACF at 200 K.
0
5
10
15
20
25
30
t IPS
Fig. 7. Constant Q plots of the self part of$he centres-of-mass intermediate scattering function for liquid methanol at 250 K. Q values (in A-‘) as displayed on the graph.
calculated the self part of the centres-of-mass intermediate scattering function, defined as F,(Q,t) = (exp{&
[k(t+s)
-k(s)
I}>
(16)
where ( - - * ) denotes an average over initial conditions, s, molecules, (x, and wave-vectors 4 of magnitude Q. Several constant Q-plots of this function computed from data of the 250 K simulation are shown in Fig. 7. From the figure it can be seen that no spurious features seem to appear at times longer than the simulation recurrence time (about 1.5 ps) [ 191; so we will improve our
161
frequency resolution computing the intermediate scattering function for times beyond this limit. The t-dependence of F,( Q,t) is qualitatively the same in the explored Qrange: a fast decay during the first picosecond followed by a much slower decay at long times. Though this function is not measurable in the laboratory, it can be compared with the translational part of the incoherent scattering function obtainable by quasielastic neutron scattering (QENS ) [ 21. Because of the different time scales (at least two) present in F, (Q,t) and for comparison purposes we have carried out a series of fittings using a biexponential model function of the form F Yde’(Q,t)=C1exp(-t/r1)+C,exp(-t/~2)
(17)
where the coefficients Ci are constrained to be positive numbers with sum equal to one. Though the model of eqn. (17) does violate several known moments of F, (Q,t) we have, so far, been unable to find a simple function able to provide a good fit to the MD data in the whole Q-range studied. Half widths at half maximum (HWHM) for the two components can be obtained by taking the inverse of zi and are displayed in Fig. 8 as a function of Q2, the Fickian regime DMDQ2is also shown for reference. The narrower component follows the Fickian line at low Q and progressively deviates from it as Q increases. For Q greater than 1.0 the wider component width has an approximate linear dependence with Q” and seems to have a non-null width in the limit Q+O; thus indicating a non-diffusive nature for this component. However, it is extremely difficult to extract any definitive conclusion because of the ill-conditioning of the fitting for low Q-values (where the wider component becomes quite unstable due to its small amplitude) and the deficiencies of the model function mentioned above. Amplitudes Ci are plotted in the second part of Fig. 8, and again we see the anomalous character of the wider component. Instead of losing intensity with increasing Q like the diffusional component, the wider one gains intensity monotonically starting from zero at Q= 0. For Q > 1.0 a triexponential model can be employed, though due to ill-conditioning the Q-dependences of the parameters worsen. The biggest change with respect to the biexponential model is that the narrower component has split into two. Finally, in order to overcome the deficiencies of the multiexponential models and to investigate to some extent the model-dependence of the results we have used one more model of the form F ,““de’(Q9t)=i$lGenp
(-$(
tz+;)““-3)
where {Ci,Zi,y}are adjustable parameters subject to the constraints
(18)
162
Ci>O
(19)
i=l
T/he functions appearing in eqn. (18) were first introduced by Egelstaff and ‘&hofield [ 201 for overcoming the deficiencies of the simple diffusion model.
0
2
1
3
4
5
6
8
7
9
10
Q*/A-*
I
I
I,,
,
1,
I,
I
I
I
I1
I
I,
-%
0.9- “, *
0.800
0.7 G 2 0.6.c a 0.55
* 0 0
0.4 -
*
0.3 -
*
** 0 x0
*
*
0
00
Cl 0
**
0
0.20.1 --E
*
arcs , 0
1
1,
I, 2
I, 3
I 4
I h’6
I
1 f
, 8
I
1 9
I 10
Q’/K’
Fig. 8. Widths and amplitudes for fitting to F,(Q,t) using the biexponential model. (top) Half widths at half maximum versus Q2. The narrower component is shown with greater detail in the inset. The straight line represents the Fickian limit DMDQ2,with DMD taken from Table 2; (bottom) amplitude versus Q2.
-*
0
-0,
-al
-b *
-CD -
N
ok
-a,\
*
c3
*
-V
-m -N
+
164
The new model is better behaved from a mathematical point of view than the multiexponential, having the right time derivative at t = 0. Results using this model are plotted in Fig. 9 where it can be seen that amplitudes are much better determined than widths or the friction coefficient, y. Again, we can see that the narrower component reduces to the Fickian at low Q. The wider component behaves qualitatively very much like the multiexponential models though we can notice how just one Egelstaff-Schofield-type function is enough to give a good description of F,(Q,t) at Q beyond 5.0 A-‘. It is known that in the Q+oo limit Egelstaff-type functions reduce to that of the free particle. Single-particle motion in liquid methanol evidences a complexity not present in simple liquids [ 17,19,20]. As we have seen there are good descriptions in the low (hydrodynamics) and high (free particle) Q limits; but for intermediate values a kinetic model is urgently needed in order to get some insight on the dynamics of HB liquids in the picosecond time scale. Collective dynamics From the centres-of-mass trajectories we have computed the spatial Fourier components of the microscopic number density according to [ 191 pd(t)=
F exp[ -i&&(t)]
CY=l
(20)
where N,,, is the number of molecules in the simulation box, and ii,(t) is the position vector of the centre of mass of molecule (Y at time t. Then the intermediate scattering function is calculated as the autocorrelation function of the Fourier components p4 ( t ) [ 191 $‘(Q,t)
+&++.Ldb))
(21)
m
where the average is taken over initial conditions, s, and wave-vectors Q of magnitude Q. The intermediate scattering function has been drawn in Fig. 10 for the four smallest Q values accessible with the 200 K simulation. The oscillations at low t are the unmistakable signature of collective excitations which become rapidly overdamped as Q increases. Some noise, a consequence of the finite simulation length, is noticeable for times longer than 3 ps. Another consequence of the finite simulation length is the observed anisotropy of the system as shown in Fig. 11, though due to the extremely long relaxation times observed in the simulated system one should perform simulations of nanosecond length in order to overcome this difficulty. However, the simulated system behaves much more isotropically when temperature is increased. The dynamical structure factor can be obtained by Fourier transforming the intermediate scattering function
165
=-&5dtF(Q,t)
S(Q,o)
exp( -iot)
(22)
--oo and it has been plotted (Fig. 12) for the lowest Q values and at 200 K. Because no windowing of F( Q,t) has been employed and F( Q,t ) is a very slow_decaying function ripples appear superimposed on the spectra. The Q=0.25 A-’ spec-
Fig. 10. Constant Q plots of the centres-of-mass intermediate scattering function for simulated liquid methanol at 200 K. Q values as shown in figure.
0.040
0 0
2
4
t/P
6
8
10
Fig. 11. Centzes-of-mass intermediate scattering function for simulated liquid methanol at 200 K and Q= 0.25 A-‘. Due to the finite simulation time length, functions calculated for different wavevectors with same magnitude are different. Simple lines, from top to bottom, Q,, &, and &; crossed line, average.
166
,Ol
x (rn’D)S w3
167 1.8
,
I
I
,
*
1.6;
1.4-
*
**o:*+ o+
lI:
+ 1.2-
+
*: *0
k : l.O0 ?
0 * + Cl 0
o"cO+
* + to
*t
0.8-
+
0.6 -
+0 0 + +o +o 0
:: +
0.4 I 0
0.5
1.0
+1J
1.5
2.0
Q/a-'
4.5-
z
4.0 3.5-
J k
3.02.5-
*** * +
*
0.5
+
0
*
0+
$0
.s-2.01.5IO-
3
0 to"
$f
+o
+ 00
+o+o
+g
8 + 0
I 0.5
I 1.0
1 1.5
I
I 2.0
Q/A_'
Fig. 13. Temperature and momentum transfer’dependence of the frequency moments for the viscoelastic component obtained fitting the dynamic structure factor with eqn. (23) for three temperatures: 200 K (empty dots); 250 K (asterisks); and 300 K (crosses).
trum exhibits a distinctive side peak at frequencies of about 1 THz that migrates to higher frequencies when Q is increased. Recent inelastic neutron scattering experiments by Bermejo et al. [ 31 have established the existence of D collective excitations in liquid methanol for momentum transfers up to 0.6 A-l. In order to compare MD results with experimental data we have performed preliminary fittings of our MD S (Q,o) data with a variant of the viscoelastic model used by Bermejo et al. [3]; the major difference being the substitution of the experimental quasielastic response by two adjustable lorentzians. Our model function then becomes
168
+iilCii rm2+~~,2,21 I
Smodel(Q,O)=CgSVisc(Q,W)
(23)
where Svisc(Q,o) corresponds to the viscoelastic component and it is given by
1211 S”“(Q,o)
=
o&C -o&z [o.r(oY -o_$)]“+
(24)
with the relaxation time defined by z -‘=2J(of
-o$)/n
(25)
Values of the frequency moments, o_+,and Ok,for the viscoelastic component are given in Fig. 13. The overall agreement with the results of Bermejo et al. is surprisingly good though experimental frequency moments are displaced to somewhat higher frequencies. Among the possible causes for the quantitative disagreement are the inadequacies of the molecular model or the quasielastic component of the fitting function is not good enough (see discussion on single particle subsection), but mainly, we are comparing MD centres of mass dynamic structure factors with the full dynamic structure factor measured with inelastic neutron scattering. These results are very encouraging and at the present moment we are working in two different directions trying to improve the agreement with the experimental data. Firstly, we are investigating other functions as models for the quasielastic component and, secondly, we are looking forward to obtaining the first spectra for the full dynamic structure factor computed from the MD trajectories. CONCLUSIONS
We have presented the results of a thorough study of centre-of-mass dynamical properties of one molecular model proposed by Haughney et al. [ 81. As expected our MD results agree with these whenever a comparison is possible (i.e. the 300 K run). A good agreement with simulations by Jorgensen [ 121 and Matsumoto and Gubbins [ 151 employing a similar molecular model is, in general, also observed. The comparison of MD predictions with available neutron scattering data of Bermejo et al. [2,3] shows that the model is able to closely reproduce the main dynamical features of real liquid methanol. The F, (Q,t) functions are reminiscent of those obtained for water. At least two time scales are present in this function: a fast decaying component dominates the dynamics during the first picosecond, and progressively a much slower component takes over at long times. Even so, the F, (Q,t) seems to indicate the existence of more than one single process which could be attributed to different
169
states of bonding. However, it is difficult to decompose the total F,(Q,t ) into different substates [ 51 since the HB admits different definitions. The memory function of the velocity autocorrelation function is too complex, and the known approximations employed for simpler liquids (Rb for example) fail since there is clearly more than one single component. Further theoretical developments and a better knowledge of the rotational dynamics (not just translational motion) will be necessary in order to gain a better understanding of the memory function. The computed dynamic structure factors are in good agreement with inelastic neutron scattering data, and this gives further support to the methodology employed to analyse the experimental intensities. A detailed analysis of the lineshapes of S (8,~) has to be performed in order to improve the deficiencies of the viscoelastic approximation. In particular the non-propagating part(central peak) has to be re-expressed in order to take into account some of the contributions that the simple exponential form for the memory function discards. Though similar to water when studied at room temperature, relaxation phenomena in liquid methanol become so slow that relaxation times are almost macroscopic (nanosecond time scale) when studied at temperatures close to the melting point (i.e. 175.4 K). This feature makes the simulation of liquid methanol increasingly difficult when reducing the temperature because of the very long runs required to obtain good statistical averages. ACKNOWLEDGEMENTS
All calculations were carried out at the Neutron Division in RAL. Thanks are given for the provision of the necessary computational resources. Members of the Computing Group are also thanked for their excellent technical assistance whenever it was required. The work was supported in part by CICYT grant No PB86-0617. REFERENCES 1 2
3
4 5 6
S.H. Chen and J. Teixeira, Adv. Chem. Phys., 64 (1986) 1. F.J. Bermejo, F. Batallin, E. Enciso, R. White, A.J. Dianoux and W.S. Howells, J. Phys., Condensed Matter, 2 (1990) 1301. F.J. Bermejo, F. Batallkn, W.S. Howells, C.J. Carlile, E. Enciso, M. Garcia-Hernandez, M. Alvarez and J. Alonso, J. Phys., Condensed Matter, 2 (1990) 5005. F.J. Bermejo, F. Bat&b, E. Enciso, M. Garcia-Hernarrdez, J. Alonso and J.L. Martinez, Europhys. Lett., 12 (1990) 129. F.J. Bermejo, F. Batall& J.L. Martinez, M. Garcia-Herndndez and E. Enciso, J. Phys., Condensed Matter, 2 (1990) 6659. G. Ciccotti and W.G. Hoover (Eds.), Molecular Dynamics simulation of Statistical Mechanical systems, North Holland, 1986. D. Bertolini, M. Cassettari, P. Grigolini, G. Salvetti and A. Tani, J. Chem. Phys., 91 (1989) 1179. M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987.
170 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
WC. Swope, H.C. Andersen, P.H. Berens and K.R. Wilson, J. Chem. Phys., 76 (1982) 637. M. Haughney, M. Ferrario and I.R. McDonald, J. Phys. Chem., 91 (1987) 4934. H.C. Andersen, J. Comput. Phys., 52 (1983) 24. D.G. Montague, J.C. Dore and S. Cummings, Mol. Phys., 53 (1984) 1049. F.J. Bermejo, and W.S. Howells, CDBOD neutron diffraction data, unpublished work, 1989. T.A. Andrea, W.C. Swope and H.C. Andersen, J. Chem. Phys., 79 (1983) 4576. W.L. Jorgensen, J. Phys. Chem., 90 (1986) 1276, and references cited therein. F. Sciortino and S.L. Fornili, J. Chem. Phys., 90 (1989) 2786. D.C. Rapaport, Mol. Phys., 50 (1983) 1151. M. Mataumoto and K.E. Gubbins, J. Chem. Phys., 93 (1990) 1981. G. Pdlinkk, E. Hawlicka and K. Heinzinger, J. Phys. Chem., 91 (1987) 4334. J.P. Boon and S. Yip, Molecular Hydrodynamics, McGraw-Hill, New York, 1982. C. Hoheisel, Comput. Phys. Rep., 12 (1990) 29. J.P. Hansen and I.R. McDonald, Theory of Simple Liquids, Academic Press, London, 1986. P.A. Egelstaff. An Introduction to the Liquid State, Academic Press, New York, 1967, p. 129. S.W. Lovesey, Theory of Neutron Scattering from Condensed Matter, Oxford University Press, Oxford, 1986.
147
H-Bond in methanol: a molecular dynamics study* J. Alonso**‘, F.J. Bermejob, M. Garcia-Hern&rdez”*l, J.L. Martinez” and W.S. Howells” “Rutherford Appleton Laboratory, Chilton, Didcot, Oxon, OX1 1 OQX (UK) *Znstituto de Estructura de la Materia, Serrano 119, E-28006 Madrid (Spain) ‘Znstitut Laue Langeoin, 156X, F-38042 Grenoble-Cedex (France) (Received 15 October 1990)
Abstract A set of molecular dynamics (MD) simulations of methanol-d, at three temperatures in the liquid range (200,250 and 300 K) have been carried out. The equations of motion of 256 molecules interacting through a potential model due to Haughney et al. [J. Phys. Chem., 91 (1987) 49341 were solved using the velocity version of the Verlet algorithm. This rather large number of molecules was required for studying the behaviour of the system at momentum transfers as low as 0.25 A-‘. It was found that the system experiences long period fluctuations, and therefore very long MD runs (of the order of 100 ps) are necessary in order to obtain accurate statistical averages. Computed static properties are in good agreement with those reported by Haughney et al. and the neutron weighted g (r) and the static structure factor compare favourably with available neutron diffraction data. The study of time-dependent properties through centre-of-mass autocorrelation functions (VACF, F. (Q, t ) and F( Q, t ) ) and their memory functions reveals features unknown in simple liquids and very similar to those found in liquid water. A close agreement between centreof-mass single-particle autocorrelation functions and the translational part of QENS data is also observed. The dynamic structure factor for the centres of mass show distinctive side peaks in the same region of the (Q, w) plane where recently collective excitations have been studied using coherent neutron scattering thus establishing the presence of propagating short wavelength modes.
INTRODUCTION
The microscopic dynamical properties of hydrogen-bonded liquids have attracted a great deal of attention in recent times, with special emphasis on the investigation of single-particle and collective dynamics [ 11. Some of the recent efforts have been devoted to the experimental investigation of such phenomena by means of neutron spectroscopy [ 2,3]. The quan*Presented at the International Symposium on Hydrogen Bond Physics held at 11Ciocco, Barga, Italy, 11-14 September 1990. ‘Permanent address: Instituto de Estructura de la Materia, Serrano 119, E-28006 Madrid, Spain.
0022-2860/91/$03.50
0 1991 Elsevier Science Publishers B.V. All rights reserved.
148
tities measured in such experiments can be related to the single-particle (i.e. low frequency rotational and translational excitations) and collective (phononlike) motions. However, difficulties in both the experimental measurements and in the interpretation of the measured spectra which have to be done in terms of fully microscopic models hinder future progress unless some additional information from other sources is available. As an example, let us recall the fact that most of the available literature concerning experimental neutron scattering data is based upon some crude models to analyze the data, mostly based on jump or diffusion-like approximations [ 41. The atomistic computer simulation methods are nowadays firmly established as a very useful tool to complement the experimental observations as well as to provide information on those quantities not amenable to experimental investigation [ 41. The purpose of the present paper is to contrast the results arising from recent neutron scattering experiments [ 2,3] with those derived from molecular dynamics (MD) simulation for a relatively simple hydrogen-bonded liquid, methanol. The single-particle dynamics has been studied by means of time-of-flight spectroscopy [ 21 and the data has been analyzed on the basis of a recently proposed model [5] where the different association states are considered explicitly. Also, the collective dynamics has been investigated by means of tripleaxis spectroscopy [ 31. In this latter case the relatively low resolution achieved in the measurement precluded the determination of direct estimates of the lineshape of the underlying dynamical structure factor; and, therefore, one of the purposes of this paper is to compare the computer simulated structure factors with those obtained by means of deconvolution of the experimental intensities. The outline of the paper is as follows. Details about the MD simulation and the molecular model are presented in the next section. All results regarding hydrogen bonds are then summarized in the section HB network analysis, whilst the following section provides results on static properties. The penultimate section is devoted to the study of centre-of-mass dynamical correlation functions. The results obtained are finally discussed in the Conclusions section. THE MD SIMULATIONS
A computer program for simulating molecular liquids was written following the hints given by Allen and Tildesley [ 61. The program generates the trajectory of a system of molecules according to Newtonian classical mechanics ( NVE-P ensemble ) and using cubic periodic boundary conditions to simulate an infinite system. The molecules are described as a collection of interaction
149
sites whose Cartesian equations of motion are integrated using the velocity version of the Verlet algorithm [ 6,7]. Because of our intention to produce long simulation runs a rigid molecular model [8] was employed to keep the computations within reasonable limits. The RATTLE algorithm [6,9] was used to implement the holonomic constraints required to keep all intramolecular distances fixed with a tolerance of 10w6. The combination of these two algorithms made it possible to use a time step as large as lo-l4 s with good energy conservation (2 x lo-l4 over periods of 30 ps). A noticeable, though small, drift in the total energy was detected for long runs and we chose to rescale particle velocities every time step in order to ensure we were really sampling the microcanonical ensemble. All results presented here were obtained using the integration scheme mentioned above with a step of 10-14 s. Because of the excellent overall agreement with experimental data we have chosen model Hl of Haughney et al. [8] for our simulations. We have introduced just two slight modifications to the original model. On the one hand, we have substituted the mass of the hydrogen atoms by that of deuterium in order to simulate the fully deuterated compound. This choice has been motivated by the availability of recent neutron scattering data for the deuterated compound [ 2,3,10] with which to compare the results. On the other hand, instead of using the Ewald method for long range dipole-dipole interactions, a switch function [ 111 was introduced to smoothly switch off all interactions between molecular pairs whose centres-of-mass separation were greater than a given cut-off distance. The interaction energy between two molecules cyand /? then becomes
(1) where Rap is the distance between centres of mass for both molecules; S(X) is the unique fifth-order polynomial satisfying S(r)=1 =o
for x
(2)
for x>R&
and having continuous
first and second derivatives at the endpoints of the
TABLE 1 Int.errnolecularpotential parameters [81 Site, i
Oii/A
(d&)(K)
Carbon Oxygen Deuterium Deuterium
3.861 3.083 0.000 0.000
91.15 87.94 0.00 0.00
(hydroxyl ) (methyl)
qil
Ie I
0.297 - 0.700 - (qc+qo) 0.000
150 TABLE 2 Details of the MD runs. Values of internal energy and pressure do not include long range corrections Quantity
Unite
Density Elapsed time
g cme3 picoseconds
Temperature IJinte’
Kelvin kJ mol-’ kJ mol-’ atmospheres atmospheres 10m6 cm2 s-i ps-2 cm-’
GE Pressure PressureLRc Diffusion coefficient Einstein frequency, C$ Qll
Run 1 0.9430
177.00 202.3 f 0.3 - 39.88 -0.37 513.0 - 193.2 0.144 463.7 114.3
Run 2 0.9400
138.95 248.1 f 0.9 -37.78 -0.37 1285.6 - 192.0 0.54 443.6 111.8
Run 3 0.8850 100.10 297.7 f 1.0
-34.79 -0.35 1639.7 - 170.2 1.69 402.2 106.4
interval [ 111. In the present calculations we took RL= 12.25 A and RU= 12.65 A for all simulations. The function V”,, (rz ,r$ ) is the potential function referred to as model Hl by Haughney et al. [ 81 and has the form of a sum of sitesite interactions (3) where Uij( 1r z - r,$ 1) is a site-site interaction which only depends on the distance, rij, between site i in molecule LXand site i in molecule j? (4) Uij=
Eij
=
(Oii+Ojj)/2
(5)
Jtii
(6)
Cjj
All required potential parameters are summarized in Table 1. Three simulations of a 256 molecule system were carried out at temperatures near 200,250 and 300 K. Values of densities and run lengths are shown in Table 2. The initial configurations for each run were generated from final configurations of previous runs after an equilibration period no shorter than 50 ps. HB NETWORK ANALYSIS
Though the debate on hydrogen bond definition is still open, it seems there are two preferred criteria: energetic [ 121 and geometric [8]. Both of them sometimes being combined with a time analysis of the hydrogen bond [ 13-151. In this work, for the sake of simplicity, we shall use a purely geometric criteria
151
in order to define the hydrogen bond. Following Haughney et al. [8] we consider two molecules as bonded if r (0 *..O) 53.5 A, r(O**.D,) 12.6 A, and the angle ~D,O~~~0130”. Percentages of molecules in different bonding states are shown in Table 3. We observe a small fraction of molecules in a “fourfold” bonded state that could possibly be regarded as molecules involved in multiple collisions and not really as being hydrogen bonded. Fractions for the most traditional bonding states agree very well with those results previously reported for similar molecular models [ 8,12,15,16]. The average number of hydrogen bonds per molecule, nne, is also included in the Table. Standard deviations for average values were estimated by computing the statistical inefficiency (SI ) associated with each quantity [ 61, namely (7)
4 (g) ) = wN3bs)1’2 h!)
where s is the statistical inefficiency, Nabsis the number of observations, and a(g) is the standard deviation of g. We have found SI values as large as lo3 steps for some of the bonding fractions; this being a consequence of the long period oscillations (of the order of 30 ps) experienced by these molar fractions. The ultimate cause of these oscillations would be the slowness of relaxation phenomena in liquid methanol due both to the existence of the hydrogen bond network and to the low temperature range studied. Cross correlations between populations of different bonded states makes it possible for the SI for nnn to be about 30 steps for all temperatures. It is also interesting to notice the high TABLE 3 Hydrogen bond statistics. All fractions expressed in percentages and hydrogen bond lifetimes in picoseconds Run 1 0.052 + 0.012 6.60 f0.17 k
f4 nHB
GiBe rHBb ?HB=
Run 2 0.44 12.2
Run 3 kO.03 f 0.3
1.84 19.8
kO.05 ?I0.3
87.8 5.51 ff0.19 0.4 0.0027 + 0.0010 1.9881+ 0.0005
79.6 7.7 f 0.5 + 0.3 0.018 f 0.004 1.9471 z!z0.0009
69.8 8.5 f0.2 kO.5 0.056 + 0.010 1.8510 I?I0.0018
9.ld 20.4d 442.4d
2.6 5.5 74.7d
1.2 2.2 2l.Od
*From the time integral of the H(t) function. bFrom the time integral of the Cc (t ) function. “From the time integral of the C,(t) function. dExtrapolated value using a multiexponential fitting.
152
5
0
10
t/ps
15
20
0.9-
0.80.706s 0.5v- 0.40.3-
0
I, 0
, IO
, 20
I
, 30
I,
40
I,, 50 60 t IPS
70
60
Fig. 1.Hydrogen-bond autocorrelation functions at different temperatures. Note the different time periods covered by each plot.
153
correlation that exists between fluctuations of temperature and nnn about their average values which at 300 K is close to 37%. In order to characterize the hydrogen-bond dynamics we have used an approach based on bond autocorrelation functions. For that purpose we have computed three different functions. 1. H(t), the fraction of bonds unbroken after time t. 2. Cc(t), the autocorrelation function for a “continuous” definition of the hydrogen bond. 3. Ci (t), the autocorrelation function for an “intermittent” definition of the hydrogen bond. Note that the first two functions require the bond to be uninterrupted. If a bond breaks and reforms again, the reformed bond is treated as a new bond; while for the third function the new bond will be recognized as an old one. Details on the computation of these functions are given elsewhere [ 8,141. Hydrogen-bond lifetimes can be computed from rnn =
I
X(t) dt
(3)
0
where X(t) is one of the three above-mentioned autocorrelation functions. Results for the bond lifetimes are shown in Table 3. In those cases when the function X(t) has not yet decayed to zero a multiexponential fitting was integrated analytically in order to obtain the bond lifetime. HB autocorrelation functions are plotted for all three temperatures in Fig. 1. Differences between H( t ) and Cc (t ) arise from the greater weight given to short-lived bonds in the former function. Because most short-lived bonds are simply binary collisions or unstable bonds that break and reform in rapid sequence, we think function Cc(t) should be preferred. However, multiexponential fittings show that both functions have similar relaxation times at large t. A recent MD study of HB in liquid methanol [ 151 has shown that Cr( t) gives the best description of the HB dynamics because time stability of the bonds is implicitly required, and therefore C,(t) depends to a lesser degree than H(t) or Cc (t ) on the exact form in which the hydrogen bond is defined. . This function shows a new picture of liquid methanol in which the relaxation of the HB network has a time scale at least an order of magnitude bigger than previously thought, and bond lifetimes are almost in the nanosecond scale. STATIC PROPERTIES
Most of the thermodynamic properties we computed and monitored in real time (i.e. at the same time as the trajectories). The average values of temperature, intermolecular potential energy and pressure for the three runs are summarized in Table 2. The standard deviation of average temperatures was esti-
154
mated by computing the statistical inefficiency (see eqn. (7) ). Long range corrections (LRC) to the intermolecular energy and pressure were estimated by assuming total absence of spatial correlations beyond the cut-off distance [6], and are included in Table 2 for completeness. Because the standard deviation of average pressure and internal energy are smaller than the corresponding long range corrections they have not been included in the Table. Typical uncertainties in the average value of the internal energy can be roughly estimated as less than 1%. However, for the pressure averages these are far less reliable because the LRC amount to 20%, and even to 60% of the actual value are 200 K. Instantaneous values of dynamical variables were observed to experience long period (30-50 ps) oscillations which also manifest themselves in great statistical inefficiencies. Owing to its large amplitude, the influence of these oscillations on calculated average values can be noticeable unless the simulation is carried out long enough to average them. The only transport coefficient we determined was the diffusion coefficient. It was computed via the Einstein relation D=ilim$<
]ZZa(t+s)-&(a)
t-w
12)
(9)
where 2,(t) is the position of the centre of mass of molecule cy at time t; and the average is taken over molecules, (Y,and initial conditions, s. Alternatively the diffusion coefficient can also be computed from the time integral of the velocity autocorrelation function (IO) We have observed that both methods give indistinguishable results within the statistical uncertainties. However, it is necessary to obtain the velocity autocorrelation function up to times of about 10 ps in order to obtain a good agreement between both methods. Diffusion constants are given in Table 2 at the three temperatures. A comparison of the MD values with those obtained experimentally [ 21 show that the simulation value is about 25% smaller. Finally, we have calculated the Einstein frequency which is obtained from the second time derivative of the normalized velocity autocorrelation function at zero time
sL:yp
a2
+ >I (
$~t+~)-ti,(~))
(11)
t=o
where the calculated frequencies are also shown in Table 2 for reference. The static spatial correlations have been investigated via pair radial distri-
155
bution functions (RDFs). Centres-of-mass RDFs, g,, (r), for methanol at 200 and 300 K are plotted in Fig. 2. The main change with temperature is an enhancement of the structure at low temperature. The narrow peak at 3.5 A corresponds to hydrogen-bonded molecules, and then the first shell of non-hydrogen-bonded molecules appears as a second peak at 4.5 A. At longer distances g,, (r ) closely resembles the RDFs found for simple liquids. In Fig. 3 the six atom-atom RDFs involving only atoms C, 0 and Do are shown for the simulation at 200 K. These partial RDFs are very similar to those reported by Haughney et al. [8] at room temperature, except for the peak heights which have increased at 200 K. Again, the main change observed on reducing the temperature is the increase in peak heights. The existence of a
-E 1.0 UY
0.0
I--5
6
7
0
9
r/B,
(b)
0.6-
0,,/1
2,,,I,,,,,, 3
4
5
6
~,,,,,, 7
0
9
r/A
Fig. 2. (a) Centres-of-mass radial distribution functions for liquid methanol at 200 K (line) and 300 K (dote); (b) Neutron-weighted RDF at 200 K.
156
oc 23456789
r/h
Fig. 3. Set of selected atom-atom radial distribution functions for liquid methanol at 200 K.
hydrogen bond network in the liquid produces very strong correlations at short distances, as can be seen from the figure; the position of the first peak in the partial RDFs being determined by the HB geometry. To enable comparison with neutron diffraction data [lo] we have computed the neutron-weighted RDF, gN (r ), which is a combination of atom-atom RDFs given by
(12) where gij (r) is an atom-atom RDF and bi is the coherent scattering length of
157
nucleus type i. The function gN(r) (also shown in Fig. 2) exhibits a great number of maxima and shoulders corresponding to the maxima of the atomatom RDFs. It was known that this molecular model produces X-ray and neutron diffraction patterns close to those measured experimentally [ 81. The function gN(r) computed from MD has the same phase as its experimental counterpart, though some discrepancies concerning peak heights exist. CENTRE-OF-MASS
CORRELATION FUNCTIONS
The detailed microscopic information provided by MD simulations enables separate investigation of translational and rotational motions, and their contributions to the dynamic structure factor, S (&,a). In this work we have limited ourselves to study the molecular translation (centres-of-mass motion), leaving rotation for a forthcoming paper. However, the centre-of-mass correlation functions constitute a bench-mark to test the usual approximations employed in the experimental data analysis. Single-particle dynamics
The study of the centre-of-mass velocity autocorrelation function (VACF) has revealed an unsuspected complexity. From Fig. 4 we can see the influence of temperature on the VACF. As expected, the cage effect (i.e. the so-called backscattering) is weaker as temperature increases; the depth of the main minimum at 300 K being about one half the depth at 200 K. Thus indicating a lost of rigidity in the response of the liquid at frequencies of some terahertz. An interesting feature of the VACF at 200 K is the small local maximum at about
0.4
0.6
0.6
t/p5
Fig. 4. Normalized centre-of-mass velocity autocorrelation functions computed at 200 and 300 K.
158
0.25 ps which is absent at 300 K where only a shoulder is noticeable. At very short times the VACF is almost independent of the temperature for the three states studied here; this can be understood by making a short-time expansion of the normalized VACF [ 171 to obtain W(t)=+
L
+...
(13)
Values of Einstein frequencies are quoted in Table 2, where their small sensitivity to temperature variations can be noticed. Fourier transforms of VACFs, V/(U), are shown in Fig. 5. At 200 K a broad shoulder is apparent at frequencies close to the Einstein frequency; this being almost non-existent at 300 K. These temperature variations are enhanced in the spectra of the normalized force autocorrelation function (FACF) which is related to that of the VACF through the equation (14) It is very interesting to see the way in which the fine structure of the band centred at 200 cm-l, present at 200 and 250 K, disappears in the 300 K simulation making clear a qualitative change in the response of the liquid at intermediate frequencies. Assuming the dynamics of liquid methanol at frequencies of a few terahertz is dominated by the presence of an extensive hydrogen bond network, these results seem to indicate the existence of some kind of dynamical transition (possibly the percolation of the HB network) somewhere between 250 and 300 K. However, results of another recent MD simulation [ 151 suggest that simulated liquid methanol is, almost certainly, percolated at 200 K while it is not at 300 K. The memory function (MF) of the VACF, m (t ), was computed by solving the memory function equation [ 17,181 t
m(t)=C(t)-Sdsm(t-s)W(s)
(15)
0
where k(t) and 9(t) are the first and the negative of the second derivative of the normalized VACF, w(t), respectively. The memory function at 200 K is shown in Fig. 6 along with the VACF and the FACF. The most striking characteristic of the MF is its oscillatory behaviour inherited, without any doubt, from the FACF to which is closely related, at least at short times. Memory functions reported in the literature for simple liquids used to be simpler functions than the corresponding ACFs and very different from the one reported here. In order to investigate the quasielastic region of the (4,~) space we have
, 150
;, 200
004-
I
,
0
0
100
I 300
Wavenumber
I 200
Wavenumber
/ cm-’
I 400
/cm-’
I
, 300
,
, 350
,
, 400
560
600
0
0.01
: 0.03 i? 0.02
t; 0.05 i!i m 0.04
E I 0.06 L
0.07
0
“;“I
-
-
0
,
,
100
, 50
Fig. 5. Fourier transforms of normalized VACFs and FACFs for 200 and 300 K simulations.
0
2003iiT 002-
E
E 2 0.06* g! 0.05-
0.07-
0.04-
, 250 0.
,
0.02 -
, 100
0.
,
0.02-:
, 50
< 0.06> 0.04 -
0.06-
;
I
O.lO0.08-
4 O.lO-
0 2
K
0.12-
2
:
200
0.16 E o.l4-
L O.OS-
5 0.14;; 0.12-
0.16-
I,
I 100
4
, 150
,
,
.I,
, 200
300 Wavenumber
200
Wavenumber
,
,
,
,
/ cm-’
K
I
K
,
I,
, 350
500
spy
, 300
300
,
300
?? 400
/cm-’
, 250
I
,
6bO
, 400
160
ooooo00 CM
Velocity
I--.
CM
Force
+++++
CM
Velocity
ACF ACF memory
function
0.6
0.4
0.8
t I PS
Fig. 6. Normalized centre-of-mass velocity, force autocorrelation functions and memory function of the VACF at 200 K.
0
5
10
15
20
25
30
t IPS
Fig. 7. Constant Q plots of the self part of$he centres-of-mass intermediate scattering function for liquid methanol at 250 K. Q values (in A-‘) as displayed on the graph.
calculated the self part of the centres-of-mass intermediate scattering function, defined as F,(Q,t) = (exp{&
[k(t+s)
-k(s)
I}>
(16)
where ( - - * ) denotes an average over initial conditions, s, molecules, (x, and wave-vectors 4 of magnitude Q. Several constant Q-plots of this function computed from data of the 250 K simulation are shown in Fig. 7. From the figure it can be seen that no spurious features seem to appear at times longer than the simulation recurrence time (about 1.5 ps) [ 191; so we will improve our
161
frequency resolution computing the intermediate scattering function for times beyond this limit. The t-dependence of F,( Q,t) is qualitatively the same in the explored Qrange: a fast decay during the first picosecond followed by a much slower decay at long times. Though this function is not measurable in the laboratory, it can be compared with the translational part of the incoherent scattering function obtainable by quasielastic neutron scattering (QENS ) [ 21. Because of the different time scales (at least two) present in F, (Q,t) and for comparison purposes we have carried out a series of fittings using a biexponential model function of the form F Yde’(Q,t)=C1exp(-t/r1)+C,exp(-t/~2)
(17)
where the coefficients Ci are constrained to be positive numbers with sum equal to one. Though the model of eqn. (17) does violate several known moments of F, (Q,t) we have, so far, been unable to find a simple function able to provide a good fit to the MD data in the whole Q-range studied. Half widths at half maximum (HWHM) for the two components can be obtained by taking the inverse of zi and are displayed in Fig. 8 as a function of Q2, the Fickian regime DMDQ2is also shown for reference. The narrower component follows the Fickian line at low Q and progressively deviates from it as Q increases. For Q greater than 1.0 the wider component width has an approximate linear dependence with Q” and seems to have a non-null width in the limit Q+O; thus indicating a non-diffusive nature for this component. However, it is extremely difficult to extract any definitive conclusion because of the ill-conditioning of the fitting for low Q-values (where the wider component becomes quite unstable due to its small amplitude) and the deficiencies of the model function mentioned above. Amplitudes Ci are plotted in the second part of Fig. 8, and again we see the anomalous character of the wider component. Instead of losing intensity with increasing Q like the diffusional component, the wider one gains intensity monotonically starting from zero at Q= 0. For Q > 1.0 a triexponential model can be employed, though due to ill-conditioning the Q-dependences of the parameters worsen. The biggest change with respect to the biexponential model is that the narrower component has split into two. Finally, in order to overcome the deficiencies of the multiexponential models and to investigate to some extent the model-dependence of the results we have used one more model of the form F ,““de’(Q9t)=i$lGenp
(-$(
tz+;)““-3)
where {Ci,Zi,y}are adjustable parameters subject to the constraints
(18)
162
Ci>O
(19)
i=l
T/he functions appearing in eqn. (18) were first introduced by Egelstaff and ‘&hofield [ 201 for overcoming the deficiencies of the simple diffusion model.
0
2
1
3
4
5
6
8
7
9
10
Q*/A-*
I
I
I,,
,
1,
I,
I
I
I
I1
I
I,
-%
0.9- “, *
0.800
0.7 G 2 0.6.c a 0.55
* 0 0
0.4 -
*
0.3 -
*
** 0 x0
*
*
0
00
Cl 0
**
0
0.20.1 --E
*
arcs , 0
1
1,
I, 2
I, 3
I 4
I h’6
I
1 f
, 8
I
1 9
I 10
Q’/K’
Fig. 8. Widths and amplitudes for fitting to F,(Q,t) using the biexponential model. (top) Half widths at half maximum versus Q2. The narrower component is shown with greater detail in the inset. The straight line represents the Fickian limit DMDQ2,with DMD taken from Table 2; (bottom) amplitude versus Q2.
-*
0
-0,
-al
-b *
-CD -
N
ok
-a,\
*
c3
*
-V
-m -N
+
164
The new model is better behaved from a mathematical point of view than the multiexponential, having the right time derivative at t = 0. Results using this model are plotted in Fig. 9 where it can be seen that amplitudes are much better determined than widths or the friction coefficient, y. Again, we can see that the narrower component reduces to the Fickian at low Q. The wider component behaves qualitatively very much like the multiexponential models though we can notice how just one Egelstaff-Schofield-type function is enough to give a good description of F,(Q,t) at Q beyond 5.0 A-‘. It is known that in the Q+oo limit Egelstaff-type functions reduce to that of the free particle. Single-particle motion in liquid methanol evidences a complexity not present in simple liquids [ 17,19,20]. As we have seen there are good descriptions in the low (hydrodynamics) and high (free particle) Q limits; but for intermediate values a kinetic model is urgently needed in order to get some insight on the dynamics of HB liquids in the picosecond time scale. Collective dynamics From the centres-of-mass trajectories we have computed the spatial Fourier components of the microscopic number density according to [ 191 pd(t)=
F exp[ -i&&(t)]
CY=l
(20)
where N,,, is the number of molecules in the simulation box, and ii,(t) is the position vector of the centre of mass of molecule (Y at time t. Then the intermediate scattering function is calculated as the autocorrelation function of the Fourier components p4 ( t ) [ 191 $‘(Q,t)
+&++.Ldb))
(21)
m
where the average is taken over initial conditions, s, and wave-vectors Q of magnitude Q. The intermediate scattering function has been drawn in Fig. 10 for the four smallest Q values accessible with the 200 K simulation. The oscillations at low t are the unmistakable signature of collective excitations which become rapidly overdamped as Q increases. Some noise, a consequence of the finite simulation length, is noticeable for times longer than 3 ps. Another consequence of the finite simulation length is the observed anisotropy of the system as shown in Fig. 11, though due to the extremely long relaxation times observed in the simulated system one should perform simulations of nanosecond length in order to overcome this difficulty. However, the simulated system behaves much more isotropically when temperature is increased. The dynamical structure factor can be obtained by Fourier transforming the intermediate scattering function
165
=-&5dtF(Q,t)
S(Q,o)
exp( -iot)
(22)
--oo and it has been plotted (Fig. 12) for the lowest Q values and at 200 K. Because no windowing of F( Q,t) has been employed and F( Q,t ) is a very slow_decaying function ripples appear superimposed on the spectra. The Q=0.25 A-’ spec-
Fig. 10. Constant Q plots of the centres-of-mass intermediate scattering function for simulated liquid methanol at 200 K. Q values as shown in figure.
0.040
0 0
2
4
t/P
6
8
10
Fig. 11. Centzes-of-mass intermediate scattering function for simulated liquid methanol at 200 K and Q= 0.25 A-‘. Due to the finite simulation time length, functions calculated for different wavevectors with same magnitude are different. Simple lines, from top to bottom, Q,, &, and &; crossed line, average.
166
,Ol
x (rn’D)S w3
167 1.8
,
I
I
,
*
1.6;
1.4-
*
**o:*+ o+
lI:
+ 1.2-
+
*: *0
k : l.O0 ?
0 * + Cl 0
o"cO+
* + to
*t
0.8-
+
0.6 -
+0 0 + +o +o 0
:: +
0.4 I 0
0.5
1.0
+1J
1.5
2.0
Q/a-'
4.5-
z
4.0 3.5-
J k
3.02.5-
*** * +
*
0.5
+
0
*
0+
$0
.s-2.01.5IO-
3
0 to"
$f
+o
+ 00
+o+o
+g
8 + 0
I 0.5
I 1.0
1 1.5
I
I 2.0
Q/A_'
Fig. 13. Temperature and momentum transfer’dependence of the frequency moments for the viscoelastic component obtained fitting the dynamic structure factor with eqn. (23) for three temperatures: 200 K (empty dots); 250 K (asterisks); and 300 K (crosses).
trum exhibits a distinctive side peak at frequencies of about 1 THz that migrates to higher frequencies when Q is increased. Recent inelastic neutron scattering experiments by Bermejo et al. [ 31 have established the existence of D collective excitations in liquid methanol for momentum transfers up to 0.6 A-l. In order to compare MD results with experimental data we have performed preliminary fittings of our MD S (Q,o) data with a variant of the viscoelastic model used by Bermejo et al. [3]; the major difference being the substitution of the experimental quasielastic response by two adjustable lorentzians. Our model function then becomes
168
+iilCii rm2+~~,2,21 I
Smodel(Q,O)=CgSVisc(Q,W)
(23)
where Svisc(Q,o) corresponds to the viscoelastic component and it is given by
1211 S”“(Q,o)
=
o&C -o&z [o.r(oY -o_$)]“+
(24)
with the relaxation time defined by z -‘=2J(of
-o$)/n
(25)
Values of the frequency moments, o_+,and Ok,for the viscoelastic component are given in Fig. 13. The overall agreement with the results of Bermejo et al. is surprisingly good though experimental frequency moments are displaced to somewhat higher frequencies. Among the possible causes for the quantitative disagreement are the inadequacies of the molecular model or the quasielastic component of the fitting function is not good enough (see discussion on single particle subsection), but mainly, we are comparing MD centres of mass dynamic structure factors with the full dynamic structure factor measured with inelastic neutron scattering. These results are very encouraging and at the present moment we are working in two different directions trying to improve the agreement with the experimental data. Firstly, we are investigating other functions as models for the quasielastic component and, secondly, we are looking forward to obtaining the first spectra for the full dynamic structure factor computed from the MD trajectories. CONCLUSIONS
We have presented the results of a thorough study of centre-of-mass dynamical properties of one molecular model proposed by Haughney et al. [ 81. As expected our MD results agree with these whenever a comparison is possible (i.e. the 300 K run). A good agreement with simulations by Jorgensen [ 121 and Matsumoto and Gubbins [ 151 employing a similar molecular model is, in general, also observed. The comparison of MD predictions with available neutron scattering data of Bermejo et al. [2,3] shows that the model is able to closely reproduce the main dynamical features of real liquid methanol. The F, (Q,t) functions are reminiscent of those obtained for water. At least two time scales are present in this function: a fast decaying component dominates the dynamics during the first picosecond, and progressively a much slower component takes over at long times. Even so, the F, (Q,t) seems to indicate the existence of more than one single process which could be attributed to different
169
states of bonding. However, it is difficult to decompose the total F,(Q,t ) into different substates [ 51 since the HB admits different definitions. The memory function of the velocity autocorrelation function is too complex, and the known approximations employed for simpler liquids (Rb for example) fail since there is clearly more than one single component. Further theoretical developments and a better knowledge of the rotational dynamics (not just translational motion) will be necessary in order to gain a better understanding of the memory function. The computed dynamic structure factors are in good agreement with inelastic neutron scattering data, and this gives further support to the methodology employed to analyse the experimental intensities. A detailed analysis of the lineshapes of S (8,~) has to be performed in order to improve the deficiencies of the viscoelastic approximation. In particular the non-propagating part(central peak) has to be re-expressed in order to take into account some of the contributions that the simple exponential form for the memory function discards. Though similar to water when studied at room temperature, relaxation phenomena in liquid methanol become so slow that relaxation times are almost macroscopic (nanosecond time scale) when studied at temperatures close to the melting point (i.e. 175.4 K). This feature makes the simulation of liquid methanol increasingly difficult when reducing the temperature because of the very long runs required to obtain good statistical averages. ACKNOWLEDGEMENTS
All calculations were carried out at the Neutron Division in RAL. Thanks are given for the provision of the necessary computational resources. Members of the Computing Group are also thanked for their excellent technical assistance whenever it was required. The work was supported in part by CICYT grant No PB86-0617. REFERENCES 1 2
3
4 5 6
S.H. Chen and J. Teixeira, Adv. Chem. Phys., 64 (1986) 1. F.J. Bermejo, F. Batallin, E. Enciso, R. White, A.J. Dianoux and W.S. Howells, J. Phys., Condensed Matter, 2 (1990) 1301. F.J. Bermejo, F. Batallkn, W.S. Howells, C.J. Carlile, E. Enciso, M. Garcia-Hernandez, M. Alvarez and J. Alonso, J. Phys., Condensed Matter, 2 (1990) 5005. F.J. Bermejo, F. Bat&b, E. Enciso, M. Garcia-Hernarrdez, J. Alonso and J.L. Martinez, Europhys. Lett., 12 (1990) 129. F.J. Bermejo, F. Batall& J.L. Martinez, M. Garcia-Herndndez and E. Enciso, J. Phys., Condensed Matter, 2 (1990) 6659. G. Ciccotti and W.G. Hoover (Eds.), Molecular Dynamics simulation of Statistical Mechanical systems, North Holland, 1986. D. Bertolini, M. Cassettari, P. Grigolini, G. Salvetti and A. Tani, J. Chem. Phys., 91 (1989) 1179. M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987.
170 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
WC. Swope, H.C. Andersen, P.H. Berens and K.R. Wilson, J. Chem. Phys., 76 (1982) 637. M. Haughney, M. Ferrario and I.R. McDonald, J. Phys. Chem., 91 (1987) 4934. H.C. Andersen, J. Comput. Phys., 52 (1983) 24. D.G. Montague, J.C. Dore and S. Cummings, Mol. Phys., 53 (1984) 1049. F.J. Bermejo, and W.S. Howells, CDBOD neutron diffraction data, unpublished work, 1989. T.A. Andrea, W.C. Swope and H.C. Andersen, J. Chem. Phys., 79 (1983) 4576. W.L. Jorgensen, J. Phys. Chem., 90 (1986) 1276, and references cited therein. F. Sciortino and S.L. Fornili, J. Chem. Phys., 90 (1989) 2786. D.C. Rapaport, Mol. Phys., 50 (1983) 1151. M. Mataumoto and K.E. Gubbins, J. Chem. Phys., 93 (1990) 1981. G. Pdlinkk, E. Hawlicka and K. Heinzinger, J. Phys. Chem., 91 (1987) 4334. J.P. Boon and S. Yip, Molecular Hydrodynamics, McGraw-Hill, New York, 1982. C. Hoheisel, Comput. Phys. Rep., 12 (1990) 29. J.P. Hansen and I.R. McDonald, Theory of Simple Liquids, Academic Press, London, 1986. P.A. Egelstaff. An Introduction to the Liquid State, Academic Press, New York, 1967, p. 129. S.W. Lovesey, Theory of Neutron Scattering from Condensed Matter, Oxford University Press, Oxford, 1986.