- Email: [email protected]
Calculation of elastic coherent neutron scattering spectra from molecular dynamics data: The NaCN plastic crystal
8 August 1997
CHEMICAL PHYSICS LETTERS ELSEVIER
Chemical PhysicsLetters 274 (1997) 335-340
Calculation of elastic coherent neutron scattering spectra from molecular dynamics data: the NaCN plastic crystal Riccardo Chelli ~,b, Gianni Cardini a,],b, Salvatore Califano a,b a Laboratorio di Spettroscopia Molecolare, Dipartimento di Chimica, Via Gino Capponi 9, 50121 Florence, Italy b European Laborato~for Nonlinear Spectroscopy (LENS), Largo E. Fermi 2, 50125 Florence, Italy
Received 14 March 1997;in final form 6 June 1997
Abstract A recipe for the calculation of elastic coherent neutron scattering spectra from molecular dynamics simulation data is proposed and applied to the case of the plastic phase of NaCN. A good agreement with the experimental results is obtained. © 1997 Published by Elsevier Science B.V.
1. Introduction Sodium cyanide is characterized by a complex phase diagram [1]. At room pressure and temperature the stable phase (I) has sodium chloride structure [2,3], space group Fm3m and site group O h, with a dynamical orientational disorder of the C N - ions characteristic of a plastic phase. The orientational distribution is not completely random, the direction (100> being more favoured with respect to a secondary maximum along the (111> direction and to the direction (110), which is the absolute minimum [4]. By decreasing the temperature to 288 K, a transition occurs to an orthorhombic phase (II) where the C N - ions show a head-to-tail disorder [2,3] and are aligned along the former (110> direction of the cubic cell (now being the b crystallographic axis). Several experimental results are available on phase ! and the existing light [5,6], neutron scattering
1E-maih [email protected].
[4,7,8], NMR [9,10] and Raman [11] data can be interpreted only assuming a strong roto-translational coupling [12-16]. This system has been the subject of numerous theoretical investigations and molecular dynamics simulations [17-19], aimed at finding a molecular description of the dynamics of the crystal able to reproduce the available experimental data. Rowe and Susman [8] have measured the elastic coherent statical structure factor by diffuse neutron scattering to obtain information on the nature of short-range correlations in the plastic phase of the NaCN crystal. In this letter we report the calculation of the same observable by molecular dynamics simulation. We shall show that unreliable results are obtained by performing the calculation as time average unless a long simulation time is used, as shown recently by Ryckaert et al. [20] for incoherent neutron scattering, for which they adopted an improved method in the time domain. Instead good agreement with experimental data is obtained by performing the calculation in the frequency domain using window functions.
0009-2614/97/$17.00 © 1997 Published by Elsevier Science B.V. All rights reserved. PI1 S0009-26 14(97)00680-5
R. Chelli et al. / Chemical Physics Letters 274 (1997) 335-340
336
Table 1 A t o m - a t o m parameters. Potential parameter for the Buckingham type potential V~,t~(R)= A~¢e-b~ R - C ~ R -6 between atoms of type o~ and /3, at a distance R
2. The model
2.1. The interaction potential Several models [17-191 have been proposed to represent the interaction potential for NaCN which differ in the choice of both the atom-atom and the charge distribution on the C N - ion. A constant feature of these models is that the C N - quadrupole is always chosen to be smaller than the one computed using ab-initio methods [21-23]. Ferrario et al. have proposed a seven-charge distribution model [24] in good agreement with the multipolar distribution of the C N - ion computed by density functional theory by Le Sat and Gordon [21]. This potential model works well for KCN but not for NaCN [19,24], probably due to the incorrect balance of the electrostatic interactions with the atom-atom part. The potential model giving the best compromise is the simplified potential model proposed by Klein et al. [19] which, though giving a worse representation of the electrostatics with respect to the seven center charge distribution model, is able to reproduce [19] with good accuracy a large part of the experimental data available on the plastic phase. This potential model includes two distinct parts: an atom-atom contribution and an electrostatic interaction between the ions. The atom-atom potential is chosen in the form of a semi-empirical Buckingham type potential whose parameters are collected in Table 1, while the electrostatic interaction is described using a threecharge model for the CN ion as shown in Fig. 1 and a full charge for the Na + ion.
t~ - [3
A,~t~ (MJ/mol)
b,~0 (,~ - I )
C,~t~ ( k J / m o l ,~6)
Na-Na Na-C Na-N C-C C-N N-N
40.9 84.4 88.6 174.3 182.9 192.0
3.155 3.278 3.378 3.4 3.5 3.6
101.0 509.4 416.6 2569.0 2100.8 1718.0
sponding to 6 unit cells (with cell parameter of 5.875 A) along each box axis. The electrostatic interaction has been computed with the Ewald method [25] using a convergence parameter a = 0.225 ,~- 1 and a cutoff of 1.4175 A,-~ in reciprocal space and 14 A in direct space. The equations of motion were integrated in the NVE ensemble with a time step of 6 fs using the Veflet method plus constraints as implemented by Ciccotti et al. [26]. The system was prepared increasing the temperature by a combination of stochastic collisions [27] and uniform scaling[28,29] of the atomic velocities in a period of 90 ps. Before starting the run, the system was left free to evolve for a further 40 ps. The equidistribution of the energy in the system was checked and the average temperature was found to be 294.6 K. The run was of 40960 steps ( ~ 246 ps) and the configurations obtained were saved every 30 fs for subsequent analysis.
2.2. The simulation 3. Results and discussion
In the present calculations, we have chosen a cubic simulation box with 864 NaCN units corre0.Be ~
- 1.0 e ,,,,,~.
r
0.8837
M/ -0.8e
-'
Fig. 1. Partial charge distribution of the C N - ion. The distances from the center of molecular mass are expressed in ,~,ngstrom and the charges in fractions of an electron. The CN bond length is assumed to be equal to 1.17 A.
3.1. Coherent elastic neutron scattering: a brief summary of the theory An expression for the diffuse scattering can be obtained from the Van Hove density operator [30] p(K,t), for a generic wave-vector K at time t: M
p(K,t) = E bmexp[-iK'Rm(t)],
(1)
m=l
where Rm(t) is the position at time t of the m-th atom and b,, is the neutron scattering length. The
R. Chelli et al./ Chemical Physics Letters 274 (1997) 335-340
337
coherent neutron scattering cross section is proportional to the dynamic structure factor:
I I
1
S(K,to) = ~ f _
(p(K,O);p(-K,t))ei~'dt (2)
where h ~o and hK are the energy and the momentum transferred from the sample to a neutron during the scattering event, respectively, and N is a normalization factor. As usual, the symbol ( •. • ) indicates a correlation function. The static structure factor, S(K), is obtained integrating the dynamic structure factor, S(K, to), over the frequency spectrum:
S(K) =
S(K,oJ) do~ =
#'t
•i
i /
I
J i
I
,,
i
0
T,(II p(K)ll2).
~
l
2
3
4
Q
(3) Usually, S(K) is the sum of two contributions: ace(K), due to the elastic scattering:
Fig. 2. Comparison between the experimental [8] (dashed line) and the calculated elastic coherent static structure factor from Eq. (6) (full line). The dashed line is a fit to the experimental data. The data are reported along the 21r/a(Q 10) line, a being the cell parameter.
1
ace(K) = ~11< p(K))II 2
(4)
and Sca(K), due to inelastic processes: aca(K ) =S(K)
(5)
--ace(K).
3.2. Calculationfrom MD data In this Letter the attention is focused on the calculation of the elastic contribution (Eq. (4)) to the static structure factor (Eq. (3)). For a comparison with the experimental data, we will consider existing data on the sodium cyanide crystal [8]. To this purpose the density operator (Eq. (1)) has been first calculated from MD data and then inserted into Eq. (4). The coherent elastic contribution to S(K) is obtained from direct time averaging over the number T of configurations: Sc~(K) =
1
r
~ E
p(K,t)
12.
(6)
t=l
The results of the application of the direct method to our simulation data, together with the experimental data [8], are shown in Fig. 2. The agreement between the calculated and experimental data is poor. This is disappointing since one expects to be in good statistical conditions. In fact the simulated sample contains 2592 atoms, the system is properly thermalized and the sampling has been extended in time up to
245.76 ps. This is, however, clearly not sufficient for obtaining reliable results in the calculation. The value of See(K) depends strongly on the time length of the simulation over which the average is made, even for times as long as 150 ps. To further confirm the nature of the problem we have performed an average of Sce(K) over all symmetry equivalent K points and obtained that the standard deviation is often as large as the average value calculated, as shown in Fig. 3a. This shows that the time domain calculations are unreliable, at least for the plastic phase of NaCN. To overcome these statistical and numerical problems in the calculation of ace(K), the length of the simulation run should, in principle, be increased until convergence in the running averages is reached. This procedure is too expensive, in terms of computer time, to be pursued. We, therefore, suggest performing the calculation in the frequency domain. The advantage of shifting to the frequency domain is that the data can be treated with window functions to be closer to the experimental setup. The elastic structure factor is experimentally measured in the frequency domain with finite resolution. This means that the observed spectrum results from the convolution of elastic and inelastic (with small aJ value) contributions due to the finite instrumental resolution. The experiment can therefore be simulated by the use of appropriate windows around the oJ = 0 value.
338
R. Chelli et al. / Chemical Physics Letters 274 (1997) 335-340
[ ::,/,::: ':/i'I
a
::
-
.=_-
.a ,m
0
1
2
3
4
5
Q Fig. 3. Average value of the elastic coherent static structure factor calculated by averaging over all the equivalent K points along the (Q 10) direction (full lines). The (a) spectrum is calculated from Eq. (6). The (b) spectrum is calculated from Eq. (9) with a Gaussian window of 5 cm -~ . The dotted lines define the upper and lower limits of the standard deviation. Let us start observing that, for w = 0, the dynamic structure factor of Eq. (2) is the elastic contribution to the coherent scattering. The elastic coherent structure factor is expressed for w = 0 in Eq. (2) as the limiting value of the time average of the density operator autocorrelation function:
s(i¢,0) = sco ( *( ) 1
shows noisy oscillations of amplitude comparable to those observed for the imaginary part. This is the reason why in the time domain it is not possible to compute a reliable value for the coherent elastic structure factor for a generic K point. It is to be noted that this is not true if the calculation is performed at Bragg points. At these points the limiting value is, in percentage, similar to the initial one and the amplitude of the fluctuations can be neglected with respect to the average value of the density operator. The long-time behaviour of the autocorrelation functions is therefore a consequence of statistical and numerical errors, that are known to be larger at longer correlation times [31]. We propose, therefore, to compute the elastic coherent scattering in the frequency domain, as done experimentally, by computing the full S(K,w) (Eq. (2)) and selecting the zero frequency value corresponding to the elastic coherent contribution. Clearly, the values obtained in this way are still affected by a large error but they can be improved by the use of appropriate window functions in the frequency domain. The results obtained with this procedure are valid for any K point. This approach finds its justification in the use in the experiments of energy analyzers with finite resolution. We have investigated two different kinds of windows and both give results in good agreement with the experimental observation for NaCN crystal [8]: a rectangular one (Eq. (8)) and a Gaussian averaging function (Eq. (9)):
T
1
0tlimr_~- ~ ~ (p(K,O);p(-K,t)).
t~0
=
Nc __E0S ( X , o,b ,
(8)
(7) Since at thermodynamic equilibrium this autocorrelation function is real, we can neglect the imaginary part. This is completely justified since any numerical discrepancy is only due to numerical errors. Two different behaviours are observed for the real part of the autocorrelation function depending on the K value. At the Bragg points the autocorrelation functions are noisily oscillating around the zero correlation time value. On the contrary, at a generic K point they are characterized by a rapid decay, followed by small amplitude oscillations compared to the zero time value. In both cases, the long time behaviour does not converge to a constant value (corresponding to the elastic contribution) but rather
1
s o(K)
Nm
Nm
ES(K,"j)
E e j=0
;=0
e
,
,
(9)
where Nc is the number of channels over which the average is performed, N O is the maximum number of channels resulting from the Fourier transform of the density operator autocorrelation function (Eq. (2)) and vT~ln2 is the half-width at half-maximum of the Gaussian. An estimate of the error can be obtained by comparing the value at equivalent symmetry points. The results are reported graphically in Fig. 3b for a Gaussian window of 5 cm -~. With a rectangular window similar results are obtained.
R. Chelli et al. / Chemical Physics Letters 274 (1997) 335-340
"~n
,
'
--
,
0
1
,
,
t
I
t
,'
,,,;
2
3
,
4
5
Q Fig. 4. Comparison between the experimental [8] (dashed line) and the calculated elastic coherent static structure factor fi'om Eq. (9) (full line). The Gaussian window is 5 cm - i. The data are reported along the 21r/a(Q 10) line, a being the cell parameter.
Fig. 4 shows the comparison between the experimental and the calculated elastic coherent structure factors along the (410) direction. The figure shows clearly the large improvement obtained by performing the calculation in the frequency domain. The improvement in the results is due to the fact that also the contribution from the low frequency inelastic scattering present in the experimental spectra is now
339
taken into account in the calculation. To further confirm the validity of this procedure we have computed the elastic scattering for the K points in the (001) plane compatible with our simulation box. A plot is reported in Fig. 5 where points with close intensity value were connected by smooth continuous lines. The agreement with Fig. 1 of Ref. [8] is impressive; the shape and the position of the maxima are practically the same except for the (440) Bragg region where we obtain a broader feature. In conclusion, a procedure is proposed in this work to bypass numerical difficulties entrenched in the calculation of the coherent elastic structure factor through Eq. (6). By use of the dynamic structure factor at to = 0 and appropriate window functions, a good agreement with experimental data is found.
Acknowledgements
This work was supported by the Italian Ministero dell'Universitfi e della Ricerca Scientifica e Tecnologica (MURST), by the Consiglio Nazionale delle Ricerche (CNR) and by the European Union (contract No. ERBFMGECT950017).
References
[1] [2] [3] [4] [5] [6] [7]
//
Qy
[8] [9] [10] 0
'
'
1
'
Q,
3
4
5
Fig. 5, Neutron scattering intensity contour lines in the (0 0 l) plane. Reciprocal lattice vectors are expressed as 2~r/a(Qx Q v O) where a is the cell parameter. The data were computed with a Gaussian window function with a half-width at half-maximum of 5 cm - I . The intensity on the isolines are reported in arbitrary units.
[11] [12] [13] [14] [15] [16] [17]
W. Dultz, H. Krause, Phys. Rev. B 18 (1978) 394. H.J. Verweel, J.M. Bijvoet, Z. Krist. 100 (1938) 201. J.M. Bijvoet, J.A. Lely, Rec. Trav. Chim. 59 (1940) 908. J.M. Rowe, D.G. Hinks, D.L. Price, S. Susman, J.J. Rush, J. Chem. Phys. 58 (1973) 2039. C.H. Wang, S.K. Satija, J. Chem. Phys. 67 (1977) 851. S.K. Satija, C.H. Wang, J. Chem. Phys. 70 (1979) 4437. J.M. Rowe, J.J. Rush, N. Vagelatos, D.L. Price, D.G. Hinks, S. Susman, J. Chem. Phys. 62 (1975) 4551. J.M. Rowe, S. Susman, Phys. Rev. B 29 (1984) 4727. R.E. Wasylishen, B.A. Pettitt, K.R. Jeffrey, J. Chem. Phys. 74 (1981) 6022. D.E. O'Reilly, E.M. Peterson, C.E. Scheie, P.K. Kadaba, J. Chem. Phys. 58 (1973) 3018. D. Fontaine, R.M. Pick, J. Phys. (Paris)40 (1979) 1105. K.H. Michel, K. Parlinski, Phys. Rev. B 31 (1985) 1823. K.H. Michel, J.M. Rowe, Phys. Rev. B 32 (1985) 5827. K.H. Michel, J. Chem. Phys. 84 (1986) 3451. K.H. Michel, J. Naudts, J. Chem. Phys. 68 (1978) 216. J.P. Ryckaert, A. Bellemans, G. Ciccotti, Molec. Phys. 44 (1981) 979. M.L. Klein, I.R. McDonald, Chem. Phys. Lett. 78 (1981) 383.
340
R. Chelli et al. / Chemical Physics Letters 274 (1997) 335-340
[18] R.M. Lynden-Bell, I.R. McDonald, M.L. Klein, Molcc. Phys. 48 (1983) 1093. [19] M.L. Klein, I.R. McDonald, M. Ferrario, J. Chem. Phys. 84 (1986) 3975. [20] J.P. Ryckaert, M.L. Klein, I.R. McDonald, Molec. Phys. 83 (1994) 439. [21] R. LeSar, R.G. Gordon, J. Chem. Phys. 77 (1982) 3682. [22] R. Bonaccorsi, C. Petrongolo, E. Scrocco, J. Tomasi, Chem. Phys. Lett. 3 (1969) 473. [23] J.E. Gready, G.B. Bacskay, N.S. Hush, Chem. Phys. 31 (1978) 467.
[24] M.L. Klein, I.R. McDonald, J. Chem. Phys. 79 (1983) 2333. [25] S.W. De Leeuw, J.W. Perram, E.R. Smith, Proc. R. Soc. Lond. 373 (1980) 27. [26] G. Ciccotti, M. Ferrario, J.P. Ryckaert, Molec. Phys. 47 (1982) 1253. [27] H.C. Andersen, J. Chem. Phys. 72 (1980) 2384. [28] L.V. Woodcock, Chem. Phys. Lett. 10 (1971) 257. [29] F.F. Abraham, S.W. Kock, R.C. Desai, Phys. Rev. Lett. 49 (1982) 923. [30] L. Van Hove, Phys. Rev. 95 (1954) 249. [31] R. Zwanzig, N.K. Ailawadi, Phys. Rev. 182 (1969) 280.
CHEMICAL PHYSICS LETTERS ELSEVIER
Chemical PhysicsLetters 274 (1997) 335-340
Calculation of elastic coherent neutron scattering spectra from molecular dynamics data: the NaCN plastic crystal Riccardo Chelli ~,b, Gianni Cardini a,],b, Salvatore Califano a,b a Laboratorio di Spettroscopia Molecolare, Dipartimento di Chimica, Via Gino Capponi 9, 50121 Florence, Italy b European Laborato~for Nonlinear Spectroscopy (LENS), Largo E. Fermi 2, 50125 Florence, Italy
Received 14 March 1997;in final form 6 June 1997
Abstract A recipe for the calculation of elastic coherent neutron scattering spectra from molecular dynamics simulation data is proposed and applied to the case of the plastic phase of NaCN. A good agreement with the experimental results is obtained. © 1997 Published by Elsevier Science B.V.
1. Introduction Sodium cyanide is characterized by a complex phase diagram [1]. At room pressure and temperature the stable phase (I) has sodium chloride structure [2,3], space group Fm3m and site group O h, with a dynamical orientational disorder of the C N - ions characteristic of a plastic phase. The orientational distribution is not completely random, the direction (100> being more favoured with respect to a secondary maximum along the (111> direction and to the direction (110), which is the absolute minimum [4]. By decreasing the temperature to 288 K, a transition occurs to an orthorhombic phase (II) where the C N - ions show a head-to-tail disorder [2,3] and are aligned along the former (110> direction of the cubic cell (now being the b crystallographic axis). Several experimental results are available on phase ! and the existing light [5,6], neutron scattering
1E-maih [email protected].
[4,7,8], NMR [9,10] and Raman [11] data can be interpreted only assuming a strong roto-translational coupling [12-16]. This system has been the subject of numerous theoretical investigations and molecular dynamics simulations [17-19], aimed at finding a molecular description of the dynamics of the crystal able to reproduce the available experimental data. Rowe and Susman [8] have measured the elastic coherent statical structure factor by diffuse neutron scattering to obtain information on the nature of short-range correlations in the plastic phase of the NaCN crystal. In this letter we report the calculation of the same observable by molecular dynamics simulation. We shall show that unreliable results are obtained by performing the calculation as time average unless a long simulation time is used, as shown recently by Ryckaert et al. [20] for incoherent neutron scattering, for which they adopted an improved method in the time domain. Instead good agreement with experimental data is obtained by performing the calculation in the frequency domain using window functions.
0009-2614/97/$17.00 © 1997 Published by Elsevier Science B.V. All rights reserved. PI1 S0009-26 14(97)00680-5
R. Chelli et al. / Chemical Physics Letters 274 (1997) 335-340
336
Table 1 A t o m - a t o m parameters. Potential parameter for the Buckingham type potential V~,t~(R)= A~¢e-b~ R - C ~ R -6 between atoms of type o~ and /3, at a distance R
2. The model
2.1. The interaction potential Several models [17-191 have been proposed to represent the interaction potential for NaCN which differ in the choice of both the atom-atom and the charge distribution on the C N - ion. A constant feature of these models is that the C N - quadrupole is always chosen to be smaller than the one computed using ab-initio methods [21-23]. Ferrario et al. have proposed a seven-charge distribution model [24] in good agreement with the multipolar distribution of the C N - ion computed by density functional theory by Le Sat and Gordon [21]. This potential model works well for KCN but not for NaCN [19,24], probably due to the incorrect balance of the electrostatic interactions with the atom-atom part. The potential model giving the best compromise is the simplified potential model proposed by Klein et al. [19] which, though giving a worse representation of the electrostatics with respect to the seven center charge distribution model, is able to reproduce [19] with good accuracy a large part of the experimental data available on the plastic phase. This potential model includes two distinct parts: an atom-atom contribution and an electrostatic interaction between the ions. The atom-atom potential is chosen in the form of a semi-empirical Buckingham type potential whose parameters are collected in Table 1, while the electrostatic interaction is described using a threecharge model for the CN ion as shown in Fig. 1 and a full charge for the Na + ion.
t~ - [3
A,~t~ (MJ/mol)
b,~0 (,~ - I )
C,~t~ ( k J / m o l ,~6)
Na-Na Na-C Na-N C-C C-N N-N
40.9 84.4 88.6 174.3 182.9 192.0
3.155 3.278 3.378 3.4 3.5 3.6
101.0 509.4 416.6 2569.0 2100.8 1718.0
sponding to 6 unit cells (with cell parameter of 5.875 A) along each box axis. The electrostatic interaction has been computed with the Ewald method [25] using a convergence parameter a = 0.225 ,~- 1 and a cutoff of 1.4175 A,-~ in reciprocal space and 14 A in direct space. The equations of motion were integrated in the NVE ensemble with a time step of 6 fs using the Veflet method plus constraints as implemented by Ciccotti et al. [26]. The system was prepared increasing the temperature by a combination of stochastic collisions [27] and uniform scaling[28,29] of the atomic velocities in a period of 90 ps. Before starting the run, the system was left free to evolve for a further 40 ps. The equidistribution of the energy in the system was checked and the average temperature was found to be 294.6 K. The run was of 40960 steps ( ~ 246 ps) and the configurations obtained were saved every 30 fs for subsequent analysis.
2.2. The simulation 3. Results and discussion
In the present calculations, we have chosen a cubic simulation box with 864 NaCN units corre0.Be ~
- 1.0 e ,,,,,~.
r
0.8837
M/ -0.8e
-'
Fig. 1. Partial charge distribution of the C N - ion. The distances from the center of molecular mass are expressed in ,~,ngstrom and the charges in fractions of an electron. The CN bond length is assumed to be equal to 1.17 A.
3.1. Coherent elastic neutron scattering: a brief summary of the theory An expression for the diffuse scattering can be obtained from the Van Hove density operator [30] p(K,t), for a generic wave-vector K at time t: M
p(K,t) = E bmexp[-iK'Rm(t)],
(1)
m=l
where Rm(t) is the position at time t of the m-th atom and b,, is the neutron scattering length. The
R. Chelli et al./ Chemical Physics Letters 274 (1997) 335-340
337
coherent neutron scattering cross section is proportional to the dynamic structure factor:
I I
1
S(K,to) = ~ f _
(p(K,O);p(-K,t))ei~'dt (2)
where h ~o and hK are the energy and the momentum transferred from the sample to a neutron during the scattering event, respectively, and N is a normalization factor. As usual, the symbol ( •. • ) indicates a correlation function. The static structure factor, S(K), is obtained integrating the dynamic structure factor, S(K, to), over the frequency spectrum:
S(K) =
S(K,oJ) do~ =
#'t
•i
i /
I
J i
I
,,
i
0
T,(II p(K)ll2).
~
l
2
3
4
Q
(3) Usually, S(K) is the sum of two contributions: ace(K), due to the elastic scattering:
Fig. 2. Comparison between the experimental [8] (dashed line) and the calculated elastic coherent static structure factor from Eq. (6) (full line). The dashed line is a fit to the experimental data. The data are reported along the 21r/a(Q 10) line, a being the cell parameter.
1
ace(K) = ~11< p(K))II 2
(4)
and Sca(K), due to inelastic processes: aca(K ) =S(K)
(5)
--ace(K).
3.2. Calculationfrom MD data In this Letter the attention is focused on the calculation of the elastic contribution (Eq. (4)) to the static structure factor (Eq. (3)). For a comparison with the experimental data, we will consider existing data on the sodium cyanide crystal [8]. To this purpose the density operator (Eq. (1)) has been first calculated from MD data and then inserted into Eq. (4). The coherent elastic contribution to S(K) is obtained from direct time averaging over the number T of configurations: Sc~(K) =
1
r
~ E
p(K,t)
12.
(6)
t=l
The results of the application of the direct method to our simulation data, together with the experimental data [8], are shown in Fig. 2. The agreement between the calculated and experimental data is poor. This is disappointing since one expects to be in good statistical conditions. In fact the simulated sample contains 2592 atoms, the system is properly thermalized and the sampling has been extended in time up to
245.76 ps. This is, however, clearly not sufficient for obtaining reliable results in the calculation. The value of See(K) depends strongly on the time length of the simulation over which the average is made, even for times as long as 150 ps. To further confirm the nature of the problem we have performed an average of Sce(K) over all symmetry equivalent K points and obtained that the standard deviation is often as large as the average value calculated, as shown in Fig. 3a. This shows that the time domain calculations are unreliable, at least for the plastic phase of NaCN. To overcome these statistical and numerical problems in the calculation of ace(K), the length of the simulation run should, in principle, be increased until convergence in the running averages is reached. This procedure is too expensive, in terms of computer time, to be pursued. We, therefore, suggest performing the calculation in the frequency domain. The advantage of shifting to the frequency domain is that the data can be treated with window functions to be closer to the experimental setup. The elastic structure factor is experimentally measured in the frequency domain with finite resolution. This means that the observed spectrum results from the convolution of elastic and inelastic (with small aJ value) contributions due to the finite instrumental resolution. The experiment can therefore be simulated by the use of appropriate windows around the oJ = 0 value.
338
R. Chelli et al. / Chemical Physics Letters 274 (1997) 335-340
[ ::,/,::: ':/i'I
a
::
-
.=_-
.a ,m
0
1
2
3
4
5
Q Fig. 3. Average value of the elastic coherent static structure factor calculated by averaging over all the equivalent K points along the (Q 10) direction (full lines). The (a) spectrum is calculated from Eq. (6). The (b) spectrum is calculated from Eq. (9) with a Gaussian window of 5 cm -~ . The dotted lines define the upper and lower limits of the standard deviation. Let us start observing that, for w = 0, the dynamic structure factor of Eq. (2) is the elastic contribution to the coherent scattering. The elastic coherent structure factor is expressed for w = 0 in Eq. (2) as the limiting value of the time average of the density operator autocorrelation function:
s(i¢,0) = sco ( *( ) 1
shows noisy oscillations of amplitude comparable to those observed for the imaginary part. This is the reason why in the time domain it is not possible to compute a reliable value for the coherent elastic structure factor for a generic K point. It is to be noted that this is not true if the calculation is performed at Bragg points. At these points the limiting value is, in percentage, similar to the initial one and the amplitude of the fluctuations can be neglected with respect to the average value of the density operator. The long-time behaviour of the autocorrelation functions is therefore a consequence of statistical and numerical errors, that are known to be larger at longer correlation times [31]. We propose, therefore, to compute the elastic coherent scattering in the frequency domain, as done experimentally, by computing the full S(K,w) (Eq. (2)) and selecting the zero frequency value corresponding to the elastic coherent contribution. Clearly, the values obtained in this way are still affected by a large error but they can be improved by the use of appropriate window functions in the frequency domain. The results obtained with this procedure are valid for any K point. This approach finds its justification in the use in the experiments of energy analyzers with finite resolution. We have investigated two different kinds of windows and both give results in good agreement with the experimental observation for NaCN crystal [8]: a rectangular one (Eq. (8)) and a Gaussian averaging function (Eq. (9)):
T
1
0tlimr_~- ~ ~ (p(K,O);p(-K,t)).
t~0
=
Nc __E0S ( X , o,b ,
(8)
(7) Since at thermodynamic equilibrium this autocorrelation function is real, we can neglect the imaginary part. This is completely justified since any numerical discrepancy is only due to numerical errors. Two different behaviours are observed for the real part of the autocorrelation function depending on the K value. At the Bragg points the autocorrelation functions are noisily oscillating around the zero correlation time value. On the contrary, at a generic K point they are characterized by a rapid decay, followed by small amplitude oscillations compared to the zero time value. In both cases, the long time behaviour does not converge to a constant value (corresponding to the elastic contribution) but rather
1
s o(K)
Nm
Nm
ES(K,"j)
E e j=0
;=0
e
,
,
(9)
where Nc is the number of channels over which the average is performed, N O is the maximum number of channels resulting from the Fourier transform of the density operator autocorrelation function (Eq. (2)) and vT~ln2 is the half-width at half-maximum of the Gaussian. An estimate of the error can be obtained by comparing the value at equivalent symmetry points. The results are reported graphically in Fig. 3b for a Gaussian window of 5 cm -~. With a rectangular window similar results are obtained.
R. Chelli et al. / Chemical Physics Letters 274 (1997) 335-340
"~n
,
'
--
,
0
1
,
,
t
I
t
,'
,,,;
2
3
,
4
5
Q Fig. 4. Comparison between the experimental [8] (dashed line) and the calculated elastic coherent static structure factor fi'om Eq. (9) (full line). The Gaussian window is 5 cm - i. The data are reported along the 21r/a(Q 10) line, a being the cell parameter.
Fig. 4 shows the comparison between the experimental and the calculated elastic coherent structure factors along the (410) direction. The figure shows clearly the large improvement obtained by performing the calculation in the frequency domain. The improvement in the results is due to the fact that also the contribution from the low frequency inelastic scattering present in the experimental spectra is now
339
taken into account in the calculation. To further confirm the validity of this procedure we have computed the elastic scattering for the K points in the (001) plane compatible with our simulation box. A plot is reported in Fig. 5 where points with close intensity value were connected by smooth continuous lines. The agreement with Fig. 1 of Ref. [8] is impressive; the shape and the position of the maxima are practically the same except for the (440) Bragg region where we obtain a broader feature. In conclusion, a procedure is proposed in this work to bypass numerical difficulties entrenched in the calculation of the coherent elastic structure factor through Eq. (6). By use of the dynamic structure factor at to = 0 and appropriate window functions, a good agreement with experimental data is found.
Acknowledgements
This work was supported by the Italian Ministero dell'Universitfi e della Ricerca Scientifica e Tecnologica (MURST), by the Consiglio Nazionale delle Ricerche (CNR) and by the European Union (contract No. ERBFMGECT950017).
References
[1] [2] [3] [4] [5] [6] [7]
//
Qy
[8] [9] [10] 0
'
'
1
'
Q,
3
4
5
Fig. 5, Neutron scattering intensity contour lines in the (0 0 l) plane. Reciprocal lattice vectors are expressed as 2~r/a(Qx Q v O) where a is the cell parameter. The data were computed with a Gaussian window function with a half-width at half-maximum of 5 cm - I . The intensity on the isolines are reported in arbitrary units.
[11] [12] [13] [14] [15] [16] [17]
W. Dultz, H. Krause, Phys. Rev. B 18 (1978) 394. H.J. Verweel, J.M. Bijvoet, Z. Krist. 100 (1938) 201. J.M. Bijvoet, J.A. Lely, Rec. Trav. Chim. 59 (1940) 908. J.M. Rowe, D.G. Hinks, D.L. Price, S. Susman, J.J. Rush, J. Chem. Phys. 58 (1973) 2039. C.H. Wang, S.K. Satija, J. Chem. Phys. 67 (1977) 851. S.K. Satija, C.H. Wang, J. Chem. Phys. 70 (1979) 4437. J.M. Rowe, J.J. Rush, N. Vagelatos, D.L. Price, D.G. Hinks, S. Susman, J. Chem. Phys. 62 (1975) 4551. J.M. Rowe, S. Susman, Phys. Rev. B 29 (1984) 4727. R.E. Wasylishen, B.A. Pettitt, K.R. Jeffrey, J. Chem. Phys. 74 (1981) 6022. D.E. O'Reilly, E.M. Peterson, C.E. Scheie, P.K. Kadaba, J. Chem. Phys. 58 (1973) 3018. D. Fontaine, R.M. Pick, J. Phys. (Paris)40 (1979) 1105. K.H. Michel, K. Parlinski, Phys. Rev. B 31 (1985) 1823. K.H. Michel, J.M. Rowe, Phys. Rev. B 32 (1985) 5827. K.H. Michel, J. Chem. Phys. 84 (1986) 3451. K.H. Michel, J. Naudts, J. Chem. Phys. 68 (1978) 216. J.P. Ryckaert, A. Bellemans, G. Ciccotti, Molec. Phys. 44 (1981) 979. M.L. Klein, I.R. McDonald, Chem. Phys. Lett. 78 (1981) 383.
340
R. Chelli et al. / Chemical Physics Letters 274 (1997) 335-340
[18] R.M. Lynden-Bell, I.R. McDonald, M.L. Klein, Molcc. Phys. 48 (1983) 1093. [19] M.L. Klein, I.R. McDonald, M. Ferrario, J. Chem. Phys. 84 (1986) 3975. [20] J.P. Ryckaert, M.L. Klein, I.R. McDonald, Molec. Phys. 83 (1994) 439. [21] R. LeSar, R.G. Gordon, J. Chem. Phys. 77 (1982) 3682. [22] R. Bonaccorsi, C. Petrongolo, E. Scrocco, J. Tomasi, Chem. Phys. Lett. 3 (1969) 473. [23] J.E. Gready, G.B. Bacskay, N.S. Hush, Chem. Phys. 31 (1978) 467.
[24] M.L. Klein, I.R. McDonald, J. Chem. Phys. 79 (1983) 2333. [25] S.W. De Leeuw, J.W. Perram, E.R. Smith, Proc. R. Soc. Lond. 373 (1980) 27. [26] G. Ciccotti, M. Ferrario, J.P. Ryckaert, Molec. Phys. 47 (1982) 1253. [27] H.C. Andersen, J. Chem. Phys. 72 (1980) 2384. [28] L.V. Woodcock, Chem. Phys. Lett. 10 (1971) 257. [29] F.F. Abraham, S.W. Kock, R.C. Desai, Phys. Rev. Lett. 49 (1982) 923. [30] L. Van Hove, Phys. Rev. 95 (1954) 249. [31] R. Zwanzig, N.K. Ailawadi, Phys. Rev. 182 (1969) 280.