NEMD “q-kick response” method for the evaluation of the dynamic liquid structure factor

Journal of Molecular Liquids, 54 (1992) 321-332

321

Elsevier Science Publisher~ B.V., Am~terdRm

N E M D "q-Kick Resl nse" method for the evaluation of the Dynamic Liquid Structure Factor W . A . ]3. E v a n s a n d C . J. O j e d n , Physics Laboratory, University of Kent, Canterbury,

Kent CT2

7NR,

UK.

Abstract A new non-equilibrium molecular dynamics (NEMD) technique to obtain the Van-Hove d y n a m i c l i q u i d s t r u c t u r e f a c t o r , S ( q , ~ ) , is p r e s e n t e d . The technique depends on evaluating the system's density response to a de17sity-couplln~ time-impulsive and space-sinusoidA1 field of wave v e c t o r , q, b y t h e p e r t u r b a t i o n - a n d - d i f f e r e n c e method. Results pertaining to a L e n n a r d - J o n e s f l u i d m o d e l o f l i q u i d A r g o n a r e p r e s e n t e d as a n i l l u s t r a t i o n . 1.

INTRODUCTION We here present a new MD technique for the evaluation of the dynamic and static structure factors, S(q,t~) and S(q), by what we call the "q-Kick Response" method The idea for this occurred to one of us (WABE) some years ago I and some preliminary results were then obtained. However, the work then begun was not finalised and written up till now. As the method is p r o m i s i n g , quite distinct from o t h e r a,~4 c l a s s i c a l techniques for evaluating t h e s a m e q u a n t i t i e s a n d w i t h its p o t e n t i a l v i r t u a l l y u n e x p l o r e d , w e s h a l l p r e s e n t it's p r i n c i p l e s f u l l y h e r e w i t h a n i l l u s t r a t i v e e x a m p l e o f its a p p l i c a t i o n to the Lennard-Jones Fluid model of liquid Argon. 2.

THEORETICAL DEVELOPMENT The linear response of a system to a perturbation that couples to t h e d e n s i t y o f t h e s y s t e m is c l o s e l y r e l a t e d t o t h e Y a n - H o v e ~ s D y n a m i c L i q u i d S t r u c t u r e F a c t o r , S ( q , e ) , a n d h e n c e a l s o t o S ( q ) [ = f d e S ( q , t a ) 1. T h e r e l a t i o n , w h i c h is v a l i d e v e n q u a n t u m m e c h a n i c a l l y , is b e s t u n d e r s t o o d in teauls of the exact relation between the spectral function, D(qja), of the density linear response kernel that governs the ultrasonic response o f t h e s y s t e m ( f o r e x a m p l e ) - a n d S ( q , ~ ) , w h i c h is t h e s p e c t r a l f u n c t i o n of the space time density-density correlation function viz. S(q,ta)----

D(q,¢0)/(1

- e -~!~)

This has been shown explicitly, for example, in an earlier paper o f E v e n s a n d P o w l e s 6. T h e a b i n t t i o t h e o r e t i c a l d e v e l o p m e n t we give (below) will also make the above relation transparent. AJthough our

0167-7322/92/~05.00

01992 - - Elsevier Science Publishers B.V. All righta reserved

:]22

m e t h o d is technique ourselves principles

easily generalised to heterogenous systcans (when an equivalent can yield the partta/ structure factors), we currently confine t o a s i m p l e 1 - c o m p o n e n t f l u i d , a s o u r a i m h e r e is t o p r e s e n t t h e underlying the new technique. For a 1-component fluid we have

S(q,to) =

1 ~.

J'i

2dxt

e "t

< ~'q(t) ~'_q(O)>

(2.1)

where the circumflex denotes operators representation. Specifically, for any operator, 6(0

=

elJ(t/hOe

-D(t/h

pq being

the Fourier

p(r)

] . ~ ~(1" - t ] ) :

:

so that

component

in 0

the

~q(t) = eDtt/hpq of the density

Heisenberg

time

e -Utt/l~

operator

viz.

(2.2)

I__ X P q e t q ' r tl q 1,ff

Le.

=

or, in second

p(r)

quantisation,

=

(2.2a)

Oq ~

X

akak~l

(2.3)

k

The above definition of S(q,t~) makes clear the conventional method, used by Rahman 2 and others for the classical evaluation of S(q,t~). B y m o l e c u l a r d y n a m i c s simulation one evaluates

and then S(q,t~) simply results by the Fourier transformation of the above time(or ensemble-) averaged correlation function. From =~_q and time translational invariance, it is e a s i l y d e d u c e d that ,t)" - F ( q , - t ) k e . t h e r e a l p a r t o f F ( q , t ) is m a n i f e s t l y e v e n i n t a n d its i m a g i n a r y p a r t is o d d . I n a n i s o t r o p i c f l u i d , o n l y t h e m o d u l u s ( n o t d i r e c t i o n ) of q can affect the answer and, classically, the xj(t)'s are taken as the real p o s i t i o n s o f t h e p a r t i c l e s a t t i m e , t. S o , if w e p u t q = - q o n t h e r.h.s, o f (2.4), w e o b t a i n a n e q u i v a l e n t f o r m o f F ( q , t ) w h i c h is i t s o w n c o m p l e x c o n j u g a t e . T h u s , c l a s s i c a l l y , F ( q , t ) is a p u r e l y r o e / a n d , h e n c e , e v e n f u n c t i o n o f t. Another w a y o f s e e i n g t h i s is t o n o t e t h a t p q c o m m u t e s with the interaction potential energy part of the Hamiltonian but not with the kinetic energy. But. classically, when one integrates over phase space, the k i n e t i c e n e r g y t e r m " d r o p s o u t " i.e. g i v e s t h e i d e a l g a s c o n t r i b u t i o n as the interaction depends only on particle positions and not on their velocities. Hence, classically evaluated, F(q,t) will be real and even, which implies t h a t its t i m e F o u r i e r T r a n s f o r m , S(q,t~), w i l l a l s o b e a r e a l a n d e v e n f u n c t i o n o f t~. T h i s is w e l l - k n o w n to be in error for the resulting function

~q

323 cannot obey the detailed balance relation, S(q,~)= e--Oh~S(q,-~), nor the "well-known" f i r s t m o m e n t (f-) s u m r u l e , J d ~ 6) S (q,ca)= h q 2 / 2 m . S c h o f i e l d v has proposed a correction scheme based on accounting for the "recoil momentum" that rectifies these deficiencies. The alternative approach proposed here needs no analogous c o r r e c t i o n s i n c e , a s w i l l b e d e m o n s t r a t e d , t h e S(q,¢~) o b t a i n e d m a n i f e s t l y obeys both detailed balance and the f-sum rule. We shall classically e v a l u a t e t h e l i n e a r r e s p o n s e o f t h e s y s t e m to a d e n s i t y - c o u p l i n g external field by the "perturbation and difference" molecular dynamic simulation technique first advocated by Ciccotti and Jaccucci s - and then utUise the e~act quantum mechanical relation between this response and S(q,¢~). T h e i d e a h e r e is t h a t , a t l e a s t i n a c o m p u t e r simulation, we can i m a g i n e o u r s y s t e m p e r t u r b e d b y a " f i e l d " t h a t c o u p l e s w i t h t h e m a s s (or particle) density in our system - in an exactly analogous way to that in which normal electric fields couple with charge density. In heterogeneous systems the field may be imagined to couple to only one particle type a n d , i n t h i s w a y , t h e r e s p o n s e a l l o w s u s t o d e d u c e t h e partial s t r u c t u r e factors in an exactly analogous way - which we plan to investigate more f u l l y i n f u t u r e w o r k . In t h i s p a p e r w e r e p o r t r e s u l t s o n l y f o r o n e s t a t e o f a "single species" Lennard-Jones system (conventionally considered as an excellent model for liquid Argon and other heavier Group 0 fluids). By conventional quantum response theory we describe the external time-dependent " f i e l d " , E ( r , t ) = -Vcp(r,t~ t h a t c o u p l e ~ w i t h t h e n u r n / ~ a r d ~ s i t y v i a a n a d d i t i o n a l p e r t u r b a t i o n t e r m , K(t), i n t h ~ H a m i l t o n i a n v i z . K(t)

~

fdr d

p (r) qJ (r, t)

~

~.

~(rj(t),t)

(2.5)

J-X

w h e r e , c l a s s i c a l l y , r j ( t ) is t h e p o s i t i o n o f t h e j~dl p a r t i c l e a t t i m e t. T h e a b o v e is q u i t e g e n e r a l . F o r o u r p u r p o s e s o f e v a l u a t i n g S (q,¢~) w e c o n s i d e r the field to be a S-function in time and sinusoidally varying in space Le w i t h a w e l l - d e f i n e d q - v e c t o r t h a t is a r b i t r a r y e x c e p t t h a t , t o p r e s e r v e t h e simulation's periodic boundary c o n d i t i o n s , it h a s t o b e c o m m e n s u r a t e w i t h t h e d i m e n s i o n s o f o u r s i m u l a t i o n b o x i.e. q

=

27r [ n,~i + n y i + n = k 1, t Lx Ly Lz ]

42.6)

w h e r e o u r b o x h a s d i m e n s i o n s L x , Ly, L z ( n o r m a l l y e q u a l ) a n d n x , n y , n = a r e a r b i t r a r y i n t e g e r s - a t l e a s t o n e o f w h i c h is n o n - z e r o . F o r o b v i o u s reasons, what we termed "q-kick response" in the title to this paper derives from this kind of field or potential which does not violate the periodicity of our simulation box. Thus, with cpOr,t)

=

Eo cos(qar) ~;(t)/q

which corresponds to a "kick"-field, o f p e r t u r b a t i o n , K(t) b e c o m e s

-Vq) = EoS(t)sin(q.r)~.

For this kind

324

K(t)

=

Eo S(t)fdrcos(q.r)p{r) q d

=

E o S(t) ~ . c o s ( q . r j ( t ) ) t:1 q

(2.6)

Now the general response to this kind of perturbation may be deduced from Von-Neumann's density matrix evolution equation ihall(t) = [ J L + K(t), ll(t) ] (2.7) ll_ ..I ~t w h i c h is e q u i v a l e n t t o S c h r ~ d i n g e r ' s e q u a t i o n b u t m o r e c o n v e n i e n t here a s w e s u p p o s e o u r s y s t e m t o s t a r t o f f f r o m a n e q u i l i b r i u m s t a t e w h e r e its i n i t i a l (t----0) d e n s i t y m a t r i x is 11(0) = o - B N / Z , / t b e i n g t h e H a m i l t o n i a n of the system in the absence of any external perturbations. The part of the r e s p o n s e t h a t is l i n e a r i n K(t) (or E o ) d o m i n a t e s p r o v i d e d E o is s u f f i c i e n t l y small - as in our simulations. Formally, (see, for example, Evans 9 ) this linear response of the density matrix may be written in the manner p i o n e e r e d b y M a t s u b a r a I°, K u b o ll e t c . ll(t) - e - i J i t / h i i ( 0 ) e iJ{t/l~ ~

_ i e-iat/hf

w h e r e it s h o u l d b e n o t e d t h a t

t .to

Y,.(tt, t2)

= e L l i t ' / l ~ K ( t i ) e -i~t~#l~

(2.8) i m p l i e s t h a t t h e l i n e a r r e s p o n s e

of any operator,

42.9)

O , is

• ,c,<,>o/ w h e r e w e h a v e ulilised the c y c l i c t n v a r t a n o e o f T r a c e . N o t e the terms o n the left o f ('2.10) represent the perturbed and u n p e r t u r b e d e x p e c t a t i o n values o f the operator, O, at t i m e t i.e. the response o f O. In our case li(0) = e - ~ / Z i.e. c o m m u t e s w i t h JL. C o n s e q u e n t l y the r.h.s, o f (2.10) m a y be written

F o r t h e c a s e c o n s i d e r e d h e r e K(ti), a s g i v e n b y (2.6), is p r o p o r t i o n a l t o 8(ti), s o t h a t i f w e l o o k at t h e l i n e a r d e n s i t y r e s p o n s e (i_e. 0 = p ( r ) ) t o t h e q - k i c k f i e l d w e f i n d , u s i n g 42.6) a n d 42.11)

~< p(r,t,>

=

I)< p(r,t)>

=

-@0(t)fdrtcos(q.rt)Tr(

ze-aX [ ' ~ ( r , t ) , p ( r D ] )

-iEo0(t, Tr(e-a"[~(r.t),Pq+P-q])

42.12, 42.13,

2

w h i c h , in v i e w o f m o m e n t u m c o n s e r v a t i o n in h o m o g e n o u s fluids, only t h e q ' t h a n d - q ' t h F o u r i e r c o m p o n e n t s o f p(r,t) w i l l c o n t r i b u t e . U s i n g 42-9-) and appreciating that we have isotropic symmetry so that expectation v a l u e s d o n o t d e p e n d o n t h e d i r e c t i o n o f q, w e f i n d

825

hqt)

~-

w h i c h is o u r b a s i c e q u a t i o n f o r t h e l i n e a r r e s p o n s e e v a l u a t e d d u r i n g t h e simulations by the "perturbation and differonoa" method (Ciccotti and J a c c u c c i a ) . N o t e t h a t t h i s r e s p o n s e is c o s i n u s o i d a l i n s p a c e , a n d h e n c e , if integrated over the simulation v o l u m e , is o r t h o g o n a l to "other" similar responses there may be with different q-values (commensurate, of course, with the periodicity of the simulation box)- This moans that in one simulation we may simultaneously apply several =q-kick" perturbations at the start and still be able to differentiate the various responses simply by multiplying the total local response, 8 (ta(r,t)>, by the appropriate cos(q.r) and integrating o v e r t h e e n t i r e b o x v o l u m e . T h i s sug~Bests w e l o o k a t t h e t t t f f e r e m c e i n g(q,t) =-

2q < [dr EoN J

p(r,t)cosCq.r)>

= - 2___2__q__~ . ( c o s ( q . r l ( t ) ) ) E o N i-x

(2.15)

a s e v a l u a t e d i n a " p e r t u r b e d " r u n ( w h e r e a q - k i c k i m p u l s e f t e l d is a p p l i e d at t h e s t a r t ) a n d a l s o a s e v a l u a t e d f o r t h e u n p e r t u r b e d system for time "segments" starting from the s~me initial state. The mean of g(q,t) in the l a t t e r c ~ o s h o u l d b e z e r o a n d its F l u c t u a t i o n s w o u l d b e n e g l i g i b l e i n a g e n u i n e l y m a c r o s c o p i c s y s t m , N -> ~,. H o w e v e r , i n M D s i m u l a t i o n s , N is typically 100-~1000 particles, and here the fluctuations are generally c o m p a r a b l e o r l a r g e r t h a n t h e l i n e a r r e s p o n s e i t s e l f - w h i c h is w h y C i c c o t t i and Jaccucci's "perturbation and difference" m e t h o d a is e s s e n t i a l to distinguish the response from the ongoing fluctuations. Hence the response, t h e l i n e a r p a r t o f w h i c h is d e f i n e d i n (2`14), b e c o m e s g d t f r ( q , t ) = -2_2_q_ EoN

(oas(q.ri(t))

-- c o . ~ ( q . r l ( t ) ) )

(2.16)

whore r~(t) signifies the t'th particle's position at time t since the start of t h e " p e r t u r b e d " r u n , a n d r~(t} is t h e o o r r e s p o n m n 8 q u a n t i t y daring the u n p e c t u r b e d r u n . H e n c e , f r o m 42-14) a n d n o t i n g t h a t J d r c o s a ( q . r ) = t ) / 2 , w e s e e t h e l i n e r r e s p o n s e f o r e a c h "1rick" q - v a l u e r e d u c e s t o gdlf:q.t)

=

NhiO(t)Tr(~__e-aX[~q(t). "~_q(0)])

(2.17)

which is intimately related to S(q,m) as M i n e d in (2`1).T h e conne~;tion is most clearly seen by comparing their spectral forms, obtained by expanding the above expectation values in terms of the complete basi~ set of eigenvectors o f Jt i.e. I n > s u c h t h a t J t l n > = E n i n > . T h u s e x p a n d e d 42.1) b e c o m e s S(q,t~)

:'-~S-~ dt e'~t <~q(t)

~t_q(0)>

2x

==

o-ar"mlia N m.,, Z

Em)/h)

(2.18)

326

is c l e a r l y s e e n to b e a p o s i t i v e s e m i d e f i n i t e q u a n t i t y at all t o - v a l u e s . L i k e w i s e f r o m (2.17),gdifr(q,t), w r i t t e n i n t e r m s o f t h e e i g e n s t a t e s o f M., is g d l r r ( q , t ) = iO(t) ~ [ o -I~Em e -13E" (Em-E,,)t/l~

Nl~

.

Z-

-I I12 e'

i.e. h a s t h e s p e c t r a l f o r m i 0(t) ~

g dtff(q,t)=

dto e - i t ° t D ( q , ( o )

42.19)

where, clearly D(q,(0) =

~ ~

[ e - i ~ E ' ~ e-I~E" ] I ( m l p q l n > l 2 S ( o - t e n -

Em)/i~)

N o t e r e v e r s i n g t h e m , n l a b e l s i n t h e a b o v e s h o w s D(qda) t o b e m a n i f e s t l y real a n d o d d in t0. B e c a u s e o f t h e 8 - f u n c t i o n t h e a b o v e is, e q u i v a l e n t l y , D(q,o) =

1 X . . (I3Emz [1 - e-~'~'] II 2 S ( o - t E n - Em)/h) N

a n d s o is s i m p l y r e l a t e d to S(q, toX as g i v e n in (2.18), viz. D(q, to) = [ 1 - e-~h~]S(q,to) (2.~0) W e n o t e t h e s i g n o f D(q,¢~) is a l w a y s t h e s a m e as t h a t o f to a n d that, b e c a u s e D ( q , ~ ) is a n o d d f u n c t i o n o f ¢~, S ( q , ~ t h u s o b t a i n e d , m a n i f e s t l y o b e y s t h e d e t a i l e d b a l a n c e r e l a t i o n , S ( q , ~ ) = e - ~ n ~ S ( q , - ~ ) . H e n c e ( 2 ~ 0 ) c a n b e r e a d as D(q,co) = S(q,(o)+ S ( q , - o ) i.e.D(q,(0) is ( t w i c e ) t h e o d d part o f S(q,ca). T h e e x a c t r e l a t i o n s h i p , (2_20), h a s p r e v i o u s l y b e e n e x p l o i t e d b y E v a n s a n d P o w l e s 6 i n r e l a t i n g u l t r a s o n i c r e s p o n s e p r o p e r t i e s to n e u t r o n scattering and X-ray measurements. W e n o w s h o w w h y t h e f - s u m r u l e is o b e y e d . T h i s f o l l o w s b e c a u s e it is e a s y to d e d u c e t h e f o r m o f g d i f f ( q , t ) f o r v e r y s h o r t t i m e s Le. just a f t e r t h e q - k i c k i m p u l s e field, E o S ( t ) s i n ( q a 0 ~ , w h i c h n a t u r a l l y p r o d u c e s a n i m p u l s i v e c h a n g e o f v e l o c i t y , E o s i n ( q . r l ( 0 ) ) ~ / m , to t h e i'th p a r t i c l e o f m a s s m s i t e d at rlC0 ) w h e n t h e i m p u l s e w a s a p p l i e d . F o r a s h o r t d u r a t i o n after ( b e f o r e c o l l i s i o n s c l o u d t h e issue), t h e d i f f e r e n c e , ri'(t ) - ri(t), b e t w e e n t h e i'th p a r t i c l e ' s p o s i t i o n in t h e " p e r t u r b e d " a n d " n o n - p e r t u r b e d " e v o l u t i o n s is e n t i r e l y a c c o u n t e d f o r b y t h i s c h a n g e o f v e l o c i t y i.e. r i ' ( t ) - ri(t) = tEo~ sin(q.rl(0))/m Le_ t h e d i f f e r e n c e is l i n e a r in L H e n c e (2.16) g i v e s N

gmff(q,t) =-

2 q ,=~ ( - 2 s i n ( ~ q . ( r j ' ( t ) ÷ rj(t))) s i n ( ~ q . ( r t ' ( t ) - r j ( t ) ) ) ) EoN

w h i c h , in v i e w o f t h e above., is, a s y m p t o t i c a l l y for s m a l l t i m e s , g i v e n b y N

gdir~q,t)=

tq a 2 ~-'(sin2(q.rj(0))) m N J~i

=

tq 2 m

327

s i n c e t h e m e a n o f o v e r all p a r t i c l e s w i l l b e 1 / 2 . H o w e v e r , f o r f i n i t e N 4 = 2 5 6 in t h e s i m u l a t i o n s described i n S e c t i o n 3 b e l o w ) , s m a l l random departures from this mean value depending on the initial positions would be expected in the initial gradients of the various gdiff(q,t) response segments and, indeed, these were observed. The scatter in the initial slopes was typically (for N=2S6) about 3~ and, when avera~;ed over many segments, t h e a b o v e p r e d i c t e d m e a n i n i t i a l s l o p e o f q a / m w a s , a s e x .p~:ted, o b t a i n e d v e r y a c c u r a t e l y . F r o m 42.19), w e a l s o s e e t h a t , f o r t s m a l l , e -1~°t --~ 1 - ic0t , a n d , a s D(q,co) is o d d i n co, w e o b t a i n t o l e a d i n g o r d e r gatrr(q,t)

=

t

O~

1~fJ_.~ cka (o D(q,~a)

C o m p a r i s o n s h o w s t h e f i r s t m o m e n t o f D ( q ~ a ) is h q a / m . S i n c e , a s a r g u e d a b o v e , t h e o d d p a r t o f S(q,¢o) is ~D(q,¢o), w e o b t a i n

Le. S(q,¢o), t h u s o b t a i n e d f r o m g d i f f ( q , t ) , m a n i f e s t l y o b e y s t h e f - s u m r u l e . T h o u g h o u r e v a l u a t i o n o f g d f f r ( q , t ) b y s i m u l a t i o n is c l a s s i c a l , w e e m p h a s i z e t h a t , if w e c o u l d e v a l u a t e t h e l a t t e r q u a n t u m mechanically, the above procedure to deduce the dynamic structure factor would be e~ract/y the same. 3.

SIMULATION

DETAILS

AND

RESULTS

g d t f i ~ q , t ) w a s e v a l u a t e d b y M D s i m u l a t i o n a n d t h e i n v e r s e o f 42.19) w a s u s e d t o o b t a i n D ( q ~ ) a n d h ~ n c ~ S ( q ~ ) u s i n g (2,20). E x t e n d i n g the d e f i n i t i o n o f g d i r f ( q , t ) t o n e g a t i v e t i m e s a s a n o d d f u n c t i o n o f t, t h e i n v e r s e o f (2.19) y i e l d s D(q,~a) =

Le.

- ih foe(it e l°t gdtff(q,t) 27[ a_~

=

h fo~(it gdifr(q,t)sin(c0t) Jo

roe S(q,(0) ~--- (1 -1 e -°h'°) ~ J o dt g d t f f ( q , t ) s i n ( ( o t )

(3.1)

(3.2)

T h o u g h o u r e v a l u a t i o n o f g d i f f ( q , t ) b y M D s i m u l a t i o n is classical, w e e m p h a s i z e that (3.2) is a x a c t - i.e. t h e s a m e p r e s c r i p t i o n w o u l d g i v e us t h e eoract q u a n t u m S(q,t0) f r o m t h e e x a c t q u a n t u m g d i f f ( q , t ) - i f o n l y w e knew how to evaluate the latter (without approximation). I n o u r s i m u l a t i o n s w e u s e d , a s is c o n v e n t i o n a l , V e r l e t ' s u n i t s f o r t h e Lennard-Jones fluid i.e. l e n g t h is i n u n i t s o f o, e n e r g y i n u n i t s o f e a n d mass in units of ml48. The potential used was

which

is a " c u t - o f f " L e n n a r d - J o n e s

potential at rout = 30 with the constants

328

a a n d b c h o s e n s u c h t h a t b o t h V ( r ) a n d its d e r i v a t i v e ( t h e f o r c e ) v a n i s h e s a t the cut-off as this improves the numericai accuracy. We used the GearN o r d s i e c k 5"th o r d e r a l g o r i f l ' l m I w i t h a t i m e s t e p o f .05 ~ V e r l e t " t i m e u n i t s 0.e. 1/(moa/48e)') for most of our simulations. A cubic box containing 256 particles and of side length 6.718a correspond!nB to a number density of 0 ~ 1 4 2 o - a w a s u s e d in t h e s i m u l a t i o n s . T h e m e a n t e ~ - n p e r a t u r e w a s 0.72r.. W e c h o s e t h i s s t a t e s i n c e it w a s t h e o n e u s e d b y L e v e s q u e " V e r l e t a n d K u r k i j a r v i 4 i n t h e i r c o r r e l a t i o n f u n c t i o n a p p r o a c h t o a c l a s s i c a l S(q,t~}. P u t t i n g i n s o m e w i d e l y u s e d v a l u e s viz~ e~ 119~K, o = 0 _ 3 4 0 5 n m a n d m=40a.m.u., that describe liquid Argon rather well as a Lennard-Jones F l u i d ( s e e , f o r e x a m p l e , H i r s c h f e l d e r , C u r t i s s a n d B i r d ~ p. 1110), t h e s e v a l u e s a r e s e e n to c o r r e s p o n d to l i q u i d A r g o n at a b o u t ~ :and d e n s i t y 1.42 g m s / c c . Le. a s t a t e c l o s e t o t h e c r i t i c a l p o i n L P u t t i n g i n t h o s e A r g o n v a l u e s w e f i n d t h a t h = 0 . 2 0 5 ~ / m e o a / 4 8 ". T h i s v a l u e w i l l b e n e e d e d , a c c o r d i n g t o (39.), i n o r d e r t o d e d u c e t h e S(q,t~) c u r v e s f r o m t h e g d t f f ( q , t ) ' s . T h u s w h e r e a s g d i f f ( q , t ) m a y b e c l a i m e d to b e p u r e l y a p r o p e r t y o f a c l a s s i c a l L e n n a x d J o n e s f l u i d , t h e v a r i o u s S(q~aJ) p l o t s w e d e d u c e f r o m t h e m c a n n o t - a s t h e y require in addition a value for h to be specified - and would, differ for e x a m p l e if w e w e r e t o a p p l y t h e L e n n a r d J o n e s m o d e l t o o b t a i n S ( q ~ ) f o r (say) Xenon or Krypton. We would conjecture that the same"Lennard-Jones" gdlff(q,t) can account for the S(q~)'s of correspond/rig liquid states of the noble gases with the only difference bein~ the detailed balance d!~tortion t h a t d e p e n d s o n t h e p a r a m e t e r hl~lm~'oai48" - a n d w e p l a n t o i n v e s t i g a t e t h i s i n t h e n e a r f u t u r e . O f c o u r s e , i f w e t a k e t h e l i m i t h ~ 0 , (39.) w o u l d still b e m e a n i n g f u l a n d y i e l d s w h a t w e m i g h t t e r m Scl(q,t~) = D ( q , ~ ) / [ ~ h t ~ . W e b e l i e v e t h i s is w h a t s h o u l d c o r r e s p o n d w i t h t h e r e s u l t s o f p r e v i o u s " c l a s s i c a l " simulations of S(q,~) as done, for example, by Levesque, Verlet and K u r k i j a r v i 4. T h i s d e f i n i t i o n o f Sc~(q,t~ ) c l e a r l y i m p l i e s t h a t

S(q,t0) = [$ht~ Sci(q,~)/(l - e -~h~) (3.4) recovers our previous S(q,t~) f o r m which, as w a s s h o w n , satisfies detailed balance a n d the f - s u m rule. This w o u l d appear to be a n alternative (and possibly simpler) a p p r o a c h to understanding Schofield's v prescription. C}iven our b o x dimensions, the possible c o m m e n s u r a t e q-values are the v m 2 o u s posiP.on vectors of points in q-space foJ'ming a cubic lattice of side len~th 0.9352o -t Le. the lowest f e w h a v i n g m o d u l l of 0.9352, L323, 1.620, 1.8"/0,2.091, 2.449,... etc. M o s t l y w e simulated with just a single q-kick perturbation (with a typical strength of E o = 0.01 - at w h i c h v a l u e the observed response w a s virtually entirely linear) over s e g m e n t s of about 7.SVerlet time units (of ~/(mo2/48e) ") w h i c h w a s sufficient for gmff(q,t) to rise a n d decay d o w n exept possibly at the v e r y lowest q-value. W e also applied u p to 8 slmuultaneous q-kicks obtaining virtually identical results, the slight penalty being some additional "noise" due to interference between the various responses since they are only perfectly orthogonal in an infinite N system. A s far as could be judged, the time average of this "noise" w a s zero - so that this m o d e of w o r k i n g is certainly e c o n o m i c a l o n c o m p u t e r time a n d is to be r e c o m m e n d e d . Typical s e g m e n t s of gdiff(q,t) are s h o w n in Figure I w h e r e it c a n b e seen that the actual responses rose linearly initially

329

with time with the gradient exF, o c t e d f r o m t h e f - s u m r u l e . A t .02 larger times they differed much more in our 255 particle system .01 so that we had to average over very many (~25 or more) to get a 0 r e a s o n a b l e m e a n o f g d i f f ( q , t ) - It might be thought these wide -.01 variations at the larger times might be due to the cumulation of numerical round-off errors but t h i s is n o t s o - a s w a s p r o v e d b y F i g u r e 1. 5 0 " s e g m e n t s " o f g d i f f ( q , t ) (with their mean shown in heavy line) simulating the same segment in double precision and obtaining a virtually identical shape_ Thus the curves we obtained are 8enuino linear responses of the small system we simulated - and the large variations in these responses at the larger times illustrate how sensitively these depend on the initial state. The averages of 8dif~q.t) are shown in Figure 2 for q - v a l u e s o f 0_935, L87, 4~65 a n d 8~42 o -~. T h e y w e r e l e a s t s q u a r e s f i t t e d t o a Gaussian × [Odd Polynomial} Fiffunction of the form gdiff(q,t)

=

NL

n.~ b~e -sat2/2 H~_,(at)

(3.5)

w h e r e H ~ _ l ( a t ) a r e o d d H e a - m i t e P o l y n o m i a l s - a, b e i n g t h e " i n v e r s e - w i d t h " o f t h e G a u s s i a n , is t h e o n l y n o n - l i n e a r p a r a m e t e r i n t h e fit. T h e a b o v e form has the advantage that with a comparatively few linear tet~s, NL (~5), a f i t t o g d U ~ q , t ) C~'m b e o b t a i n e d t h a t is a s g o o d a s t h e u n c e r t a i n t y i n t h e d a t a . F u r t h e r its F o u r i e r - L a p l a c e T r a n s f o r m , (3.2), t o o b t a i n D ( q , ~ ) a n d h e n c e S(q,(~) is a n a l y t i c a n d l e a d s t o S(q,t0)

I

(1 - h e - ~ h ~ ) Y r-~--~al

-.~, b . (--1)" e - ° u / 2 a 2 H 2 ~ , ( t ~ / a )

(3J5)

In the above, we have to introduce a value for h and, as mentioned e a r l i e r , w e u s e d t h e v a l u e o f 0 7 - 0 5 t h a t is a p p r o p r i a t e t o $ r g o n . T h e p l o t s on the r~. in Figure 2 show the resulting S(q,~) [full llne] and Sol(q,~ ) [da-~hed l i n e ] o b t a i n e d u s i n g t h e g d w f { q , t ) d a t a s h o w n i n t h e p l o t t o i t s immediate left. The results are in very reasonable agreement with measured dynarnic structure factors for Argon (Henshaw~ P a g e t4, s e e E g e l s t a f f L~ a l ~ ) . A t t h e l o w e s t q o f 0 . 9 3 5 o- l w e s e e c l e a r l y t h e r e m n a n t s o f t h e w e l l k n o w n Rayleigh and Brillouin peaks that are well-defined in the long wavelength (hydrodynamic) l i m i t , s e e , f o r e x a m p l e , M o u n t a i n ~6. I n d e e d t h e i r p o s i t i o n w o u l d a p p e a r t o i n d i c a t e a v e l o c i t y o f s o u n d o f a b o u t 0 . 7 " / u n i t s o f ~/(48eYm)" or about 842m~-L in excellent (possibly fortuitousl) agreement with the mea~'red experimental value, 863m.s-~ for the sound velocity in liquid Arg-. at 83~1K (see RowHnson~V). At higher q-,~aluesthe Brlllouin peaks get s m a l l e r u n t i l a t l a r g e q a s i n g l e l a r g e p e a k d o m i n a t e s w h i c h is w h a t h a p p e n s

830

gdtrr(q.t)

S(q,o) ( qo = 0.93,.q

) and Scl(q,~ ) (- - -)

.02 -

.01 ."

_

-..

0 0

0

• 0

qo = 1.87 .03 .02

.01 0 •

J

II

II

•

0

0

-2

t---~

-1

0

1

2 (0->

qo=4.6B .05

.1

o ~

"

=.8

•

O

. •

0

i

i

[

1

. . . . . . . . .

•

2

. . . . . . . . .

I

3

. . . . . . . . .

I

4

. . . . . . . . .

t->

-5

0

¢~--> 5

qo=8.42

.5

.2

0 0

1

2

3

4

t->

'o -

0

~

~t~ --->

Figure 2. g,.ff(q,t) averaged over typically 25 segments with corresponding S(q,o) and Scl(q,~0) plots (Argon Parameters, h= 0.205)

e x p e r i m e n t a l l y . O f c o u r s e , a l l t h e S(q,t~) c u r v e s o b e y t h e f - s u m r u l e . F o r e a c h S(q, to) c u r v e o b t a i n e d , a v a l u e o f S ( q ) = j ' d t o S(q,to), c a n b e deduced which can be compared with the S(q) obtained either experimentally or by the normal MD method of evaluating static structure factor, g(r), and Fourier Transformation. In MD, g(r) can be evaluated up to modestly large r-values depending on box size. We show in Figure 3 the S(q) curve o b t a i n e d i n t h i s m a n n e r - w h i c h , i n c i d e n t a l l y , is i n f a i r a g r e e m e n t w i t h t h e A r g o n S ( q ) p l o t ( a t t h e l a r g e r q - v a l u e s w h e r e it is r e a s o n a b l y a c c u r a t e ) . A t the smaller q-values the S(q) obtained v i a g ( r ) is i n h e r e n t l y 2 i n a c c u r a t e - a n d it s e e m s l i k e l y that the q-kick response values are more trustworthy here. The q - k i c k v a l u e s a r e generally t o o small at q> 3 as the points marked in Figure 3 show. We belie~,e that this may be due to 0 the (3aussian in the parametric 0 5 10 q _~ 1 5 form causing the fit to underestimate the contribution F i g u r e 3. S ( q ) c o m p a r i s o n s from the larger q-values. Since q - k i c k v a l u e s (o) S ( q ) = .J'dt~ S(q,ta) is e q u i v a l e n t t o S(q) derived from g(r) (~) S(q) =

kaT j~odt gditt(q,t) coth(xtkBT/h)

43.7)

it is c l e a r t h a t i f t h e g d t f f ( q , t ) f i f f u n c t i o n c u t s o f f t o o q u i c k l y a s i g n i f i c a n t contribution to S(q) would be lost ( as the coth factor will be unity for large t). It c a n b e a r g u e d , a t l e a s t i n t h e h y d r o d y n a m i c limit, that the asymptotic c u t - o f f is a d a m p e d e x p o n e n t i a l - a n d p l a u s i b l y t h i s is t r u e i n g e n e r a l . T h i s a s p e c t is a l s o c u r r e n t l y b e i n g f u r t h e r i n v e s t i g a t e d i n N > 2 5 6 s i m u l a t i o n s , and through the use of alternative fiffunctions whose asymptotic behaviour at large times more closely resembl~ the decay of gdiff(q,t). 4. ACKNOWLEDGEMENTS It is b o t h a p l e a s u r e a n d a p r i v i l e d g e t o a c k n o w l e d g e the constant help, encouragement and stimulation, both on this paper's contents and al~o on diverse other topics, that one of us ObrABE) has received from Professor Jack G. Powles over the past many years whilst we worked together. His inexhaustible enthusiasm and capacity for research remains undiminished and has been an inspiration to all who collaborated, were supervised or were taught by him. We wish him very many more years of fruitful research during his "retirement". W A B E is a l s o g r a t e f u l f o r p a s t d i s c u s s i o n s w i t h D r . M . P . A l l e n a n d P r o f . (3. R i c k a y z e n o n t h e t o p i c s o f t h i s p a p e r . C J O g r a t e f u l l y a c k n o w l e d g e s the support of a European Social Fund Studentship.

332

REFERENCES 1. The author discussed

2. 3. 4, 5, 6.. 7, 8. 9, 10. ll_ 12. 13. 14_ lS_ 16_ 17.

this method with Dr_ M. P_ Allen during the mid 80’s who, subsequently, alluded to it in the NEMD section of his book viz, M, P_ Allen and D_ 3, Tildesley, **Computer Simulation of Liquids”, Oxford Sci, Pub (1987) see page 253A, Rahman in “Neutron Inelastic Scattering” (International Atomic Ener y Agency, ViennaA Vol_ I (1968), 561. D-G. %o unds, M. L. Klein and G. N_ Patey. J_ Chem_ Phys, 72 (1980). 5348. D, Levesque. L_ Verlet and J, Kurkijarvi, Phys. Rev, A7 (1973). 1690_ L_ Van Hove, Phys, Rev_ 95 (1954). 249_ W_ A. B. Evans and J. G. Powles, J_ Phys- A, 7 (1974). 1944. P_ Schofield, Phys. Rev_ Lett, 4. (1960), 239, G_ Ciccotti and G_ Jaccucci. Phys, Rev. Let&, 35 (1975), 789, W_ A_ B. Evans, Proc_ Phys. Sod, 88 (1967). 723. T_ Matsubara, Prog. Theor, Phys. 14 (1955). 351, R_ Kubo, J. Phys. Sot. (Japan) 12 (1957), 570, J_ 0, Hirschfelder, C. F_ Curtiss and R B_ Bird, “Molecular Theory of Gases and Liquids” (1954), Wiley and Sons, G_ D_ Henshaw, Phys_ Rev, 105 (1957). 976, D_ I. Page, ‘The Structure Factor of Liquid Argon at 845K”, (1972). Atomic Energy Research Establishment Report No, AERE - R 6828 P_ A_ Egelstaff, “An Introduction to the Liquid State”, A_P_ (l967)_ R_ D. Mountain, Rev_ Mod Phys. 38 (1966), 205. J_ S. Rowlinson. “Liquids and Liquid Mixtures”, Butter-worth, (1971) p 46.