Density functional theory of solution dynamics

joumaJ of"

MOLECULAR

LIQUIDS ELSEVIER

Journal of Molecular Liquids, 65/66 (1995) 131-138

Density Functional Theory of Solution Dynamics Toyonori M u n a k a t a

Department of Applied Mathematics and Physics, K y o t o University, K y o t o 606 J a p a n

Abstract The density functional theory for classical(equilibrium) statistical mechanics is generalized to deal with various dynamical processes associated with density fluctuations in liquids and solutions. This is effected by deriving a Langevin-diffusion equation for the density field. As applications of our theory we consider density fluctuations in both supercooled liquids and molecular liquids and transport coefficients.

I. Introduction The density functional theory(DFT) 1-4 has become a useful method to study freezing transition from a quantitative standpoint.

The DFT is now generalized to deal with complex systems such as

molecular liquids s and liquid crystals 4. It is to be noted that the DFT is basically concerned with the (quasi)equilibrium density profile neq(r) and the corresponding free-energy Feq --- F[neq(r)], where

F[n(r)] denotes the free-energy functional of the system and plays the central role in the DFT. If we could introduce dynamics to the DFT and follow time evolution of the density field n(r,t), this might enable us to study some dynamic processes such as density fluctuations and transport coefficients. In this paper we present a time-dependent(TD) DFT together with some (preliminary) applications of the TD-DFT. In Section 2 the DFT is briefly reviewed and in Section 3 a Langevin-diffusion(L-D) equation for the density field n(r,t) is presented. In Section 4 we consider, as applications of the TD-DFT, (A) density fluctuations in liquids and solutions and (B)mass flow around a fixed particle to calculate transport coefficients of liquids. Section 5 contains some remarks and summary of this paper.

II. D e n s i t y F u n c t i o n a l T h e o r y ( D F T )

First let us consider a one-component liquid, whose density profile is denoted as n(r).

0167-7322J95/$09.50 9 1995 Elsevier Science B.V. All rights reserved. SSD/0167-7322 (95) 00847-0

From a

132 variational principle 2,3, we obtain an equation to determine the equilibrium density profile

$f[n]/6n(r)- p + r

neq(r) :

= 0,

(1)

where Celt(r) is the external field and p is the chemical potential. We may say that all the difficulties of a many-body problem are embedded in the free-energy functional methods heretofore 3. All the approaches divide (excess) part - ~ [ n ] , with

Fid exactly

F[n] into

F[n] and

with h the thermal wavelength [h 2/2rmkBT]

1/2. Standard

Taylor expansion around the uniform liquid state n(r) = =

Fid[n], and

an interaction

given by

Fid[n] = kBT f drn(r)[ln(n(n)h a) -

.[n]

there have been proposed many

an ideal gas part

(2)

1],

perturbational approaches employ a functional

nL, a

constant,

E(1/rn[)fdrl.-.fdrm6m~/6n(rl)...6n(rrn)lnL(n(rl)--nL)...(n(rm)--nL)

=_ (kBT)E(1/m!) f dr,.../drm

era(r1,'" ,rm)(n(rl)

which is usually truncated after m = 2. e In this case

F[n]'~Fid[n] -

Fin] takes

-

I'lL)...(n(bmirm)

-

I'lL) ,

(3)

the form

(kBT f dvlc,[n(v,)-nL]

+ (kBT/2)fdrlfdr2e2([r2--r~l)[n(rl)--nLl[n(r2)--nL]},

(4)

where use is made of the fact that for a uniform system cl is a constant and c2 is a function of Iv2 - r l I. It is not difficult to derive the relation between c2 and the radial distribution function g(r) =

h(r)+ 1, 2,3 (5)

h(r) = c2(r) + nL f dr'h(r')c~(Ir - r'l).

Equation (5) is the familiar Ornstein-Zernike relation and hereafter we denote the direct correlation function c2(r) as

c(r). Based on (1) and

(4) with (2), one can discuss freezing transition of one-component

liquids. 1-3 For an S-component mixture, which is specified by S density fields,

nj(r)(j = 1,... ,S),

we

have, instead of (4),

Fin] ~_ F~jFid[nj]- {kBtEif arlcl,j

[ n i ( r ) l - nj,L]

+ (kBT/2,Eij fdrl/dr2ci.j(lr2--rl,[ni(rl)--ni.L][nj(r2)--nj,L]},

(6,

where c i j ( r ) denotes the direct correlation function between the species i and j.7 For molecular liquids, the DFT is more complicated due to the constraints imposed by the fixed bond length. 5

III. Time-Dependent

Density Functional Theory (TD-DFT)

We now present a Langevin-diffusion(L-D) equation for the density field n ( r , t ) and the corresponding

133

Fokker-Planck(F-P) equation for the distribution functional f[n(r),t] and discuss general properties of the TD-DFT. s,9 By combining a number-conservation law, the D F T and theory of Brownian motion, we derived the following equations for the density n(r,t) and the m o m e n t u m density g(r,t) : cg(r,t)/Ot = - V . g ( r , t ) / m , cgg(r,t)/~ = - n ( r , t ) V 6 F / 6 n ( r , t ) -

(7)

r 0 g ( r , t) + y ( r , t ) .

(8)

where - V 6 F / 6 n ( r ; t ) represents a generalized force on a particle at r x~ and the random force f ( v , t ) satisfies the fluctuation-disspation (FD) theorem of the form

< f i ( r , t ) f j ( r ' , t ' ) > = 2 m k B T r o n ( r , t ) ~ ( r - r')6(t - t')6i~,

(9)

with i and j denoting the Cartesian components. Since we are interested in long time behavior, we employ an adiabatic approximation for (8), yielding g(r,t) = {-n(r,t)V6F/6n(r,t) + y(r,t)}/r0.

(10)

From (7) and (10) we finally obtain the L-D equation,

cgn(r,t)/Ot = - V . { - n ( r , t ) V 6 F / 6 n ( r , t ) + f ( r , t ) } / ( m r 0 ) = - 1 7 . { i s + J R } ,

(11)

with the FD theorem (9). Js and JR in (11) denote the systematic and the random current ,respectively. From the FD theorem (9) it is seen that the random current j a ( r ; t ) = f ( r ; t ) / r n r o is a multiplicative noise it and one must specify how one interprets the noise. Here for our purpose it is to be treated as an Ito type. Following a routine procedure n to derive a F-P equation from a Langevin equation, we see that the distribution functional f[n(r),t] evolves in time, with D - k B T / ( m r o ) , according to

Of~at= -/dr6n~J(f

),

J ( f ) = D { f l f V . n ( r ) V 6 F / 6 n ( r ) + V . n(r)V6 f /6n(r)} . When f is proportional to exp(-flF[n]), V . n ( r ) V 6 f / 6 n ( r ) =

V. n ( r ) { - f l f V 6 F / 6 n ( r ) }

(12) (13) =-/3fVn(r)V

6F/6n(r) and we confirm that the stationary solution is given by e x p ( - f l F ) . In other words, the L-D equation (11) actually samples,in a steady state, the density field n(r) according to the weight exp(-/~F). The multiplicativeness of the noise f ( r , t) can be interpreted based on a simple hopping diffusion model and can be related to the internal noise proposed by Mikhailov. t2 General properties of the T D - D F T are most concisely represented by the following two H-theorems the proofs of which are given elsewhere. 9 First we neglect the random current JR and consider the diffusion equation

On(r,t)/Ot = ~ D V . n(r,t)X76F/6n(r,t) = - V . ] s ( r , t ) .

(14)

First H-theorem : When the density field n ( r , t ) evolves in time according to (14), F[n] decreases in time according to dF/dt = - 0 3 D ) - t f dr { j s ( r , t ) } ~ / n ( r , t ) < 0 until ] s ( r , t ) vanishes and it holds that

6F/6n(r) = p, representing the variational condition (1) in the D F T to determine the equilibrium density

134

field for the case ~,=t(r) = 0. N o w we turn to the full L-D equation (I I). Second H-theorem : W h e n the distribution functional fin,t]evolves in time according to the F-P equation (12) the generalized free-energy functional Fa[f] defined by

Fat/] = / DnFtnlf[n;t] + kBT / Dnftn]ln(ftn]),

(15)

decreases in time monotonically according to

dFa/dt =

-

/ Dn / dr {(kBT//Dn(v))x/2:is(r) - (DkBTn(r)//)l/2v(6//6n(r))} 2 < 0

(16)

until f[n, t] takes the form

/,,In] = const, exp[-F[nl/kBT].

(17)

We note that f Dn denotes the integration over the function space of n(v). When the integrand of (16) is zero, we have (17). Comparing the two theorems it is seen that the noise JR prevents the density field

n(r,t) from b~eing trapped in one of the local minima of the functional Fin]. la The results obtained above is readily generalized to mixtures.

V. A p p l i c a t i o n s of T D - D F T

In this section we apply T D - D F T developed in the previous section to study dynamic density fluctuations (A) and transport coefficients(B). In these studies the approximation (4) or (6) for the free-energy functional is employed since we have no reliable information on higher-order direct correlation functions cn(n > 3).

(A) Density fluctuations in liquids and solutions

Inserting (6) with (2) into (11) we obtain the following L-D equation :

cOn(r,t)/Ot = nV2n - DV. n(r,t)V f dr~c([r - r~l)[n(r',t) - nL]- V . JR

(18)

We note that the second term on the rhs of (18) describes a diffusion process, which is induced by the Vlasov field VF -

-kBTfdr'c(Ir-

r'[)[n(r',t)- nL]. Equation(18) with or without the random

current T M together with its generalization to two-component system 15 and to polar liquids 16'17 have been playing important roles and now generally called a Smoluchowski-Vlasov equation. As an application of the L-D equation (18) with the F-D theorem (9), we calculated the dynamic structure factor r of a simple liquid

G(q,t) = < n(q,t)n(-q,O) > / < n(q)n(-q) > .

(19)

135

with n(q,t)=_ f dr { n ( r , t ) - nL } e x p ( - - i q - r ) / N 1/2. Fourier-transformation of (18) and the FD theorem (9) yield

On(q, t)/Ot = -7(q)n(q,t) + Ek V(q, k)n(k,t)n(q - k,t) + r
(20)

> = - 2 D q . q'6(t - t') {6q+q,,0 + n(q + q',t)/NX/2},

(21)

where 7(q) = Dq2/s(q), with s(q) denoting the static structure factor, c(q) = 1 - 1/s(q), and V(q,k) =

Dc(k)q. k / N ~/2. If we neglect effects of both nonlinearity in(20) and the multiplicativeness of the noise, n(q,t) becomes a simple Ornstein -Uhlenbeck process ~1 and G(q,t) and the dynamic structure factor G"(q,w) = (1/2) f dtG(q,t)exp(iwt) are given by Go(q,t) = exp(-v(q)t) and Gg(q, w) = 7(q)/[v(q) 2 + w2],

(22)

respectively, where the subscript 0 on G means that we regard it as the zeroth aproximation to G.

In

calculating G(q,t) based on (20) and (21) we follow a nonlinear theory of fluctuations by Mori and Fujisaka. xr

Numerical calculation is performed for a hard sphere system characterized by a packing

fraction p = r~3nL/6, is Percus-Yevick approximation is used to supply the structural information. As p becomes large we observe narrowing of the central peak, reflecting slowing down of the density fluctuations due to nonlinear coupling in (20). At p = 0.53 a dynamic 'instability' occured where an effective diffusion constant crosses zero. As is well known the hard sphere system freezes at p ~_ 0.5 and the D F T in its various version has been applied to study the transition. 19 It is interesting to note that the interaction part of the free-energy(4) gives rise to the equilibrium transition on the one hand and it also gives rise to the instability through the nonlinear coupling in the L-D equation (20). Before leaving this section it is noted that if we take the free-energy functional (6) for an S-component mixture we have a set of L-D equations for S density fields ni(r,t)(i = 1,... ,S), which lead to

On,(r,t)/Ot=D, [ V 2 n , ( r , t ) - V

.n,(v,t)V { E j / d , ' c , ~ ( l r -

r'l)6nj(r', t)}]-

V . in,,

(23)

Furthermore if we take the free-energy functional as given by Chandler et hi. 5 for molecular(polyatomic) liquids, we obtain from (11) a L-D equation for each of the atomic species constituting the molecules, which is very similar to (20). After liniarization as done in (22), we obtain a linear diffusion equation, which was used by Hirata lr to investigate the dynamic structure factors of water.

(B) Some Transport Coefficients In this subsection we put a particle at an origin in a velocity field u ( r ) and study effects of flow on a stationary density profile n,t(r). When there is no flow u(r) = o, n,t(r) is obviously given by nLg(r), with g(r) a radial distribution function. Due to the flow u(r), this equilibrium distribution is distorted and from the distortion we can calculate some transport coefficients, like viscosity 7] and friction constant ~, as we show below. We consider a one-component system and neglect random current JR in (18) to obtain

8n(r,t)lbt = - V - j , ( r , t ) ,

(24)

136

j,(v,t) = DVn - a D ~ V t f dr'v~]] (Jr - r ' l ) , ( r ' , t ) + r where v,/l(r ) = -kBTc(r) and r

+ nu,

(25)

two-body interaction, represents effects of the particle fixed at

r = o. The last term on the rhs of (25) represents particle flow due to the velocity field, u ( r ) . We are interested in a stationary density profile n , t ( r ) around the fixed particle. First we consider the equilibrium solution n , q ( r ) when u = 0. In this case the particle flow j0 = o and from (25) we readily obtain ln[gCr)] __

ln[-,,C.)/-L] = -~ f dr'v,ss(l~ -

v'l)nL[g(r')

-

1] -- flr

(26)

where the boundary condition n,q(r) --* nL as r ---,oo is taken into account. From the relation (5) between h and c, we observe immediately that (26) is equivalent to the H N C equation to determine g(r). Thus our theory, when applied to anequilibrium situation, g!ves the H N C

result for g(r).z

Now let us turn to effects of the flow field u(r) on nat(r). We consider, to be concrete, shear flow

u(r) = 7yeffi where ez denotes a unit vector in the x-direction. In this case we assume a solution of the form n , i ( r ) - nLg(r)[1 4- w(r)7/D -t- o(7)]

(27)

where o(7)/~ --+ 0 as 7 ~ 0. Inserting (27) into (24) with i)n/Ot - O, we obtain

gV2w--nLgV2 i dr'c(Ir-r'l)g(r')w(r')+Vg.Vw-nLVg.V

i dr'c(Ir-r'l)g(r')w(,')

zyg'(r)/7. (28)

with use of Fourier transformation it is seen that (28) has a solution of the form w(r) = z y a ( r )

(29)

and a(r) satisfies a complicated integro-differential equation, which we do not write down here. The final step to obtain shear viscosity 71 is to calculate the xy component of the stress tensor axy, which is expressed, on the one hand, in terms of ,7 and 7 as =! =

'7")'

(30)

and also microscopically as T

ffxy = where r

(nLI2) ] drno,(r)xyr

(31)

- de~dr and we have neglected the kinetic contribution, which is very small at a liquid

density. Use of the solution (29) together with (30) and (31) gives

r] = (n2L/2D) / dr(xy)~ g(r)r

(32)

If we fix a particle in a uniform flow u ( r ) = uoe~, we can calculate the friction constant ( from a stationary force on a fixed particle by following similar lines of reasoning as above. Our preliminary results for shear viscosity r / o f a soft-core system with r

= e ( a / r ) 12 shows that

main contributions in the integrand in (32) come from the region r ~_ (V/N) ~/a (of course in the highdensity state) and that r/increases sharply as a function of p. =_ (elkBT)I/4(Nvra/V). Although detailed

137

analysis of 77 and ~ based on our TD-DFT needs some time to be completed, it seems that our approach is promising in view of the fact that we have at the moment nearly no (microscopic) theory for r/and (, especially for dense complex systems like water.

V. S o m e R e m a r k s a n d S u m m a r y

In this paper we gave a dynamic extension of the DFT, by deriving a L-D equation (11) with the /

fluctuation-dissipation theorem (9). We showed that the stochastic equation correctly samples the density field according to the probability exp {-/~F[n]},(17), based on the second H-theorem (16). At this point we note however that our TD-DFT is phenomenological and it is desirable to have a first-principle dynamics generalization of DFT. As applications of TD-DFT, we considered density fluctuations in liquids and transport coefficients. We are currently trying to solve the L-D equation in a real space-time to study slow dynamics in supercooled liquids. This study may be considered to be a dynamic counterpart of the work by Dasgupta and Ramaswamy is and as a similar attempt we mention the work by Lust and Valls 2~ We expect that our L-D equation has many fields of application. One example is dynamics in molecular liquids, as noted at the end of subsection IV(A), which is a very important field in connection with chemical reactions in solutions but, at the same time, is very complex to deal with from first principles.

References 1 R.Evans, Adv.Phys. 28 143(1979). 2 A.D.J.Haymet, Ann.Rev.Phys.Chem. 38 89(1987). 3 D.W.Oxtoby, In 'Liquid, Freezing, and the Glass Transition', (Eds. J.P.Hansen,D.Levesque and J.Zinn-Justin) (Elsevier:New York,1990) 4 Y.Singh, Phys.Rep. 207 351(1991). 5 D.Chandler, J.D.McCoy, and S.J.Singer, J.Chem.Phys. 85 5971;5977(1986). 6 T.V.Ramakrishnan, and M.Yussouff, Phys.Rev. B19 2775(1979). 7 J.P.Hansen, and I.R.McDonald.' Theory of Simple Liquid s'(Academic, New York, 1986). 8 T.Munakata, J.Phys.Soc.Jpn. 58 2434(1989). 9 T.Munakata, Phys.Rev.E(1994) to appear 10 T.R.Kirkpatrick and P.G.Wolynes, Phys.Rev. A35 3072(1987). 11 C.W.Gardiner ' Handbook od Stochastic Methods' (Springer, Berlin,1982) 12 A.S.Mikhailov, Phys.Rep. 84 307(1989).

13 C.Dasgupta and S.C.Ramaswamy, Physica A186 314(1992). 14 T.Munakata, J.Phys.Soc.Jpn. 43 1723(1977). 15 B.Bagchi, Physica 145 A 273(1987). 16 D.F.Calefand P.G.Wolynes, J.Chem.Phys. 78 4145(1983).

17 F.Hirata, J.Chem.Phys. 96 4619(1992). 18 H.Mori and H.Fujisaka, Prog.Theor.Phys. 49 764(1973). 19 T.Munakata, J.Phys.Soc.Jpn. 59 1299(1990). 20 A.D.J.ttaymet and D.W.Oxtoby, J.Chem.Phys. 84 1769(1986). 21 L.M.Lust and V.Valls, Phys.Rev.E 48 1787(1993).