journal of MOLECULAR
, ELSEVIER
UQUIBS
Journal of Molecular Liquids 90 (2001) 205-214
I
i
i
www.elsevier.nl/locate/molliq
Density functional theory of liquids, some extensions and applications Toyonori Munakata Dpartment of Applied Mathematics and Physics, Kyoto University, Kyoto 606, Japan
The density functional theory(DFT) of nonuniform liquids plays an important role in classical many-body theory because of its (mathematical)simplicity and (physical)clarity. Some extensions and applications of DFT, including a new proposal for an m-body hypernetted-chaln(HNC) equation for the static structure of liquids are presented. © 2001 Elsevier Science B.V.All fights reserved. 1. I N T R O D U C T I O N The density-functional theory(DFT) of nonuniform fluids has been successfully applied in quantitative studies on solid-liquid transformations, interfacial and nucleation phenomena and so on.[1] It was originally developed for simple liquids, but now it is applied to study liquid crystals, polymers, defects in solids such as dislocations.[1] The purpose of this paper is to present some extensions of DFT, with which the author himself has been concerned, together with applications thereof. In Sec.2 we briefly review DFT for later convenience. In Sec.3 we consider dynamic extension of DFT by first deriving a Langevin-diffusion equation for the density field
n(r, t).
As applications of the dynamic DFT we consider polymer conformation in solvent. In Sec.4 we turn to static problems. First in See.4.1 DFT for liquids composed of deformable molecules is considered and in Sec.4.2 we propose an m-body HNC theory, which coincides
0167-7322/01/$ - see front matter © 2001 Elsevier Science B.V. All fights reserved. PII S0167-7322(01) 00123-4
2O6 with the ordinary HNC theory for m=2. At the moment we are still not fully successful in numerically solving complicated integral equations for the m(=3)-body HNC theory and we contend ourselves by showing the general structure of the theory and a dosed set of quations explicitly for the case rn = 3.
2. D E N S I T Y
FUNCTIONAL
THEORY(DFT)
D F T is a microscopic equilibrium theory for nonuniform liquids and it is based on the fact that one can define a free-energy
F[n] of a system as a functional of a density field
n(r).[1,2] The stationary or equilibrium density profile
n,t(r) is given as the solution to
the variational equation
~F/Sn(r) + ¢ ~ t ( r ) = U,
(1)
where Cext(r) denotes an external field on the system. When Eq.(1) has many solutions, as is often the case for studies of phase transitions, we recourse to the free-energy argument to decide which solution to choose as the (thermodynamically)stable one. From the above we notice that all the many-body problem is embedded in the functional form of
F[n]. Standard practice is to first decompose F[n] as an ideal gas part F/d and
an excess part
Fex and then to (Taylor)expand Fez around a reference liquid state with
uniform density no. Thus we have for a d-dimensional system
F[~] = F,~[n] + Fe~[~],
(2)
with F/~[n] =
kBT f drn(r){In[n(r)h a] - 1}),
Fex[n] - Fe{~)In] =
-kBT/2 f dr f dr' c2(Ir - r'[)An(r)An(r),
(3) (4)
where A in Eq.(3) is the de Brogli thermal wavelength. In Eq.(4) we expanded Fex upto second order in
An(r) =_n(r) - no and the Taylor coefficient c2(r) is the two-body direct
207 correlation function of the reference liquid. It is remarked here that DFT has a close relation with the structure of equilibrium fluids. That is, ff one has a reliable expression for the free-energy density functional Fin], one can derive a good equation for the radial distribution function g2(r) which represents two-body correlations in a fluid. To illustrate this interconnection based on the Percns idea[3,4], we consider a simple liquid with interparticle interaction ¢(r). First let us put a particle at the origin of the coordinate system. Then the (equilibrium) density n~q(r), which obeys the variational equation (1), just represents nog2(r). By employing the second order( or two-body) approximation (4), we have for g2(r) = 1 ÷ h2(r) the HNC eq. In g2(r) = - ¢ ( r ) / ( k s T ) + no f drh2(IP - r'l)ca(r'). = - ¢ ( r ) / ( k . T ) ÷ noh2 * c2(r).
(5)
This equation is combined with the two-body Ornstein-Zernike relation[2]
h2(r) = c2(r) + noh2 * c2(r).
(6)
to give a closed set of equation to decide g2(r) and c2(r). In Sec.4.2 we will formulate a general m-body HNC theory.
3. D Y N A M I C D F T Slow motion in dense liquids is mainly diffusive which is driven by a generalized force
-V(IF/$n(r). Under the assumption we derived some years ago the Langevin-diffusion equation[5a]
cOn(r, t)/Ot = D V . In(r, t ) V S ( F / k s T ) / S n ( r , t)] - V- JR,
(7)
where D denotes the diffusion constant and the random flow JR =-- n(r; t)l/2FR(r; t) is assumed to satisfy the fluctuation-dissipation theorem
< -~R(P, t)FR(P', ~') >-----2DS(r - r')5(t - t')I,
(8)
208 with I a d × dunit matrix. As to the general properties of the Langevin-diffusion equation (7), we comment on asymptotic behaviors of the solution. If we neglect effects of fluctuations by putting
JR ---- 0 in Eq.(7), it holds that dF[n]/dt <_ 0
(9)
and at the stationarity where dF[n]/dt = 0, we have the variational condition (1). By inclusion of effects of noise we can show[hal that the stationary distribution Peq[n] of the density profile n(r) is given by Peq[n] = exp{-j3F[n]}/Z.
Thus the long time trajec-
tory(solution) n(r, t) of Eq.(7) turns out to sample the density profile according to the Boltzmann factor without being trapped at a local minimum of the functional Fin]. As to applications of the dynamic DFT, we have considered three problems, namely, the density-density time correlation function(d-d tcf)[hb], and shear viscosity[6] of (supercooled)liquids and solvents effects on (hereto)polymer conformation[7]. Starting from the Langevin-diffusion equation (7) one can derive a closed integro-differential(mode-coupling) equation for the d-d tcf. By numerically solving the equation for a hard-core system, we could discuss how diffusion constant( instead of viscosity, an object of the ordinary modecoupling theory) changes as density of the liquid increases.[hb] Shear viscosity can be calculated by studying stationary state as follows: Let us consider a shear flow characterized by a shear rate % In the flow the two-body correlation function g2,shear(r) is not isotropic and expressed as (when 7 is very small)
g2,shear(r) = g2(r) + 7d(r),
(10)
where d(r) denotes the distortion of g2(r) due to the flow. Starting from Eq.(7) with Eq.(4), we found that g2(r) is given by the HNC theory and d(r) was obtained as a solution to a little complicated integral equation with g2(r) as an input to the equation.J6] Fortunately there are some experimental data for the distortion d(r) via nonequilibrium
2O9 molecular dynamics[8] and our result turned out to agree well with experiments up to around the freezing point. Finally in this section we consider solvent effects on polymer conformation based on Eq.(7). Our heteropolymer, which consists of P sites, are put in a liquid. So, in addition to the solvent density profile no(r,t), we have to consider P density profiles
ni(r,t)(i = 1,..,P) and the total free-energy F is a functional of these density proP E iad + F ~ . Based on this flee-energy functional, we have from Eq.(7) files, F = ~j=0 a (coupled) Langevin-diffusion equation for each density profile nj(r, t)(j = 0, 1, .., P). At the moment we have no reliable algorithm to solve Eq.(7) for long time. So to make our numerical scheme tractable we first employ a localization approximation ni(r, t) =
~(ri(t) - r)(i = 1, .., P) and change a Langevin-diffusion equation for ni(r, t) to a simple Langevin equation for ri(t), which takes the form
dri(t)/dt
= DiV,[- / drOo,,(Ir - ri(t) l)~n0(r, t)/(ksr) P
-
~
¢j,i(Iri(t) - rj(t)l)/(ksT)] + ]i(t),
(11)
1#i,1
where ¢id denotes the interaction between the species i and j. The first term on the right hand side(r.h.s) expresses the effect of interaction with the solvent and the remaining terms on the r.h.s denotes the Langevin dynamics under the Smoluchowski limit. Instead of solving the Langevin-diffusion equation for the solvent density n0(r, t), we solve the following HNC equation with the polymer sites held fixed at {r~(t)},
exp[f dr'c0.0(Ir -
neq(?']{~i(t)}) P
-
~ - ' ~ ¢ ~ / ( I r - r,(t)l)/(kBT)],
(12)
1
with c0,0(r) the direct correlation function of solvent. In summary we follow time evolution( from t to t + ~-) of the polymer-solvent system by first solving the Langevin eq.(ll) for small time interval ~- with no(r,t) held fixed. Then we solve the HNC eq.(12) for the
210 new density profile no(r, t + r) with the polymer sites held fixed at the position ri(t + r). We considered a two-dlmensional system with the soft-disc solvent. In Fig.1 we plot the initial(a) and large time(b) polymer conformation together with the density profile no(r, t), which forms a (coordination)shell structure, around polymer sites. The effective potential ¢0,i(r) is chosen to be attractive(repulsive) for the sites i = 1, .., 5(i = 6, .., 10). Reflecting the hydrophitic nature of the sites i = 1, .., 5, we observe that they expand and enclose the hydrophobic sites( except for sites 9 and 10). Although our calculation is still at a preliminary level, the hybrid procedure, which combines the Langevin dynamics with the HNC eq., seems .to work well for studying solvent effects on polymer conformation and we hope that this method may play some role in studies such as protein folding.
4. T W O - B O D Y C O R R E L A T I O N I N LIQUIDS 4.1 Molecular liquids As was remarked in Sec.2 DFT in its standard form (4) gives rise to the HNC theory for a two-body correlation function in liquids. Let us turn to (one component) molecular
Fig.1 Polymer conformation for initial(a) and after long time(b).
211 liquids, D F T for which was developed by Chandler et a1.[9] Within the framework of the site-site interaction model(SSIM)[2,9], each molecule is characterized by P sites(c~ = 1, .., P) and the important two-body correlation functions are the intermolecular radial distribution function ga,s(r) and the intramolecular one wa,s(r), which is normalized as
f drwa,s(r) = 1 and w~,~(r) = 5(r). From the above we note that wa,s(r) denotes the probability density for the site ~ of a molecule to be observed at r relative to the site a of the same molecule. So long as we know, the rigid molecular model with the restriction on the site-site length(/a,~)
=
-
(13)
has been widely used to study structure of molecular liquids. If the rigid 88IM is employed one can not deal with thermal expansion and isomerization of constituent molecules. In Ref.[10] we removed the rigidity constraint Eq.(13) and developed an integral equation theory, in which inter- and intra-molecular correlations are determined self-consistently. The numerical results for diatomic models were quite encouraging and we hope to extend the calculation to more realistic models for molecular liquids.
4.2 m - b o d y H N C t h e o r y Let us now generalize the argument to derive the HNC eq. (5) in Sec.2 and establish a theory to deal with (up to) m-body(m > 3) correlation functions in liquids. For the purpose we expand Fex up to the m-th order, instead of terminating the (Taylor)expansion at the second order as in Eq.(4), as
F [n]
Z:
(14)
2
where Fe(~)[n] contains j-th order direct correlation function cj (1, 2, .., j) with 1 denoting r l . Thus explicitly we have e.g.
212 The first step of our m-body HNC theory is to notice that for the m-body correlation function gin(l, .., m), which is normalized to unity at infinite separation, we have g,,(1, ..,m)
=
g,,-l(1,..,m -
1)
=
gl(2]l)gl(3]l, 2)g1(4]1, 2, 3)...g,(m]l, 2, .., m - 1).
(16)
If points (1,2,..,m) representing positions in space are regarded as point on time axis, Eq.(16) reminds us o f a non-Markovian stochastic process. Furthermore this non-Markovian property is similar in its origin to that in the random-walk interpretation of polymer conformation(the excluded volume effect). This is evident from the meaning of the conditional distribution, for example g1(411, 2, 3) denotes the one-body distribution function at 4 when three particles are located at 1,2, and 3. Tentatively assuming that cj(1, . . , j ) ( j = 2, .., m) are known, we follow the idea of Percus[3] (this time however maximally m - 1 particles are held fixed in space) to derive an equation for lngl(jll, 2, ..,j - 1)(j = 2, .., m) based on Eqs.(1),(2), and (14) and the fact that the external field appearing on the r.h.s, of the variational equation (1) is the sum of the field produced by particles located at 1,2,..j-1. From this we readily obtain j-1
lngl(j[1, 2, ..,j - 1) = - ~ ¢ ( i , j ) / ( k B T )
+ Fm[gl(jll, 2, ..,j - 1) - 1],
(17)
i=l ?n
F,,,[f] - Y ~ [ n ~ - l / ( i -
1)!1 ] 1'.. [ ( i -
1)'ci(l', 2', .., ( i - 1)', j ) f ( l ' ) . . f ( ( i -
1)'),
(18)
i=2
where f ( l ' ) _= gl(l']l, 2, .., j - 1)-1. Thus we have m - 1 equations,(17)(j=2,..,m) for m - 1 unknowns gl(j[1, 2, .., j - 1)(j = 2, .., m). As to the direct correlation functions which were assumed to be known above, we know that they are in a sense inverse functions of the correlation functions.J2,9] Typically this is expressed for the case m = 2 as the Ornstein-Zernike relation (6). Although the general(mbody) Ornstein-Zernike relation are too complicated to write down here, we can at least formally obtain the relation as an integral equation expressing c~ in terms of correlation
213 functions gj(1, .., j)(j = 2, .., m). The lowest Ornstein-Zernike relations(m=2,3) are most concisely expressed in term of Fourier transformation as
d2(q) = h2(q)/(1 + noh2(q)).
(19)
/'?/(ql, q2) = cs(ql, q2)G(ql, q2),
(20)
where G(ql, q2) = (i + noh2(ql))(1+ noh2(q~))(1+ n0h~(Jql + q2J)) and c3(q,,q2) denotes the Fourier transform of c3(1,2, 3) with 3 taken to be the originof the coordinate system. Similarly//(ql,q2) is defined as the Fourier transform of H(1,2,3)
=- h3(1,2,3) - h2(1, 2) - h2(2,3) - h2(1,3) -
no/d'~h~(1,~)h~(2,-~)h~(3,-~)
[h2(1, 2)h2(2, 3) + h2(1, 3)h2(3, 2) + h2(2,1)h2(1, 3)],
(21)
with h3(1,2,3) = gs(1,2,3) - 1. For m=3 we explicitly write the generalized HNC equations for gt(2Jl) = g2(1, 2) and In g1(2]1) = - ¢ ( 1 , 2)/(kBT) + noc2 * h2 + (no2/2) f dl'd2' c3(l', 2', 2) h2(1', 1)h2(2', 1).
(22)
In g(311, 2) =
(23)
-(ksT)-1(¢(ll
- 31) + ~(12 - 31)) + no f d1'c2(3, 1')h(l'Jl, 2)
+ @/2) f al'f d2%(1',2',3)h(l'll,
2)h(2'Jl, 2).
Comparing Eq.(22) with the 2-body HNC Eq.(5), we see that we have an extra threebody potential contribution produced by the density field at 1' and 2'. This corresponds to the bridge function, which is neglected in the 2-body HNC approximation.[2,4] Equation (22) seems to be first derived by Percus.[3] Now we have four integral equations (19),(20),(22) and (23) for four unknown functions c2,c3,g2 and gs or g(3]1,2) supplemented by g3(1, 2, 3) = gt(2tl)gl(311, 2)(see Eq.(16)). At the moment we are trying to solve the coupled set of equations for a one-dimensional soft-core system. We hope to report on numerical results in near future.
214 Recently DFT is gathering considerable interest in connection with dynamic extension[ll] and the glass transion.[12] In view of the fact that DFT is rather general and flexible and that we have no standard model for liquid, which allows for analytic approaches, DFT seems to be able to play important roles for better understanding of physics and chemistry of liquid-state. The author expresses sincere gratitude to Prof. F.Hirata(Institute for Molecular Science,Okazaki) and Prof.M.Kinoshita(Kyoto University) for useful discussions and comments.
REFERENCES
1. Y.Singh, Phys.Rep. 207(1991)351. 2. J.P.Hansen and I.R.McDonald, Theory of Simple LiquidsAcademic Press, New York,198( 3. J.Percus, in The Equilibrium Theory ofClassicalFluidseds. H.L.Frisch and J.L.Lebowitz Benjyamin,INC. New York 1964. 4. L.Verlet, Physica 30(1964)95: ibid.
31(1965)959. 5a) T.Munakata, Phys.Rev.E
50(1994)2347: b) T.Munakata, J.Phys.Soc.Jpn.58(1989)2434. 6. J.Araki and T.Munakata, Phys.Rev.E 52(1995)2577. See also S.Yoshida, F.Hirata and T.Munakata, Phys.Rev.E 54(1996)1763. 7. T.Takahashi and T.Munakata, Phys.Rev.E 56(1997)4344. 8. W.T.Ashurst and W.G.Hoover, Phys.Rev.A 11(1975)658. 9. D.Chandler,J.D.McCoy and S.J.Singer, J.Chem.Phys. 85(1986)5971. 10. T.Munakata, S.Yoshida and F.Hirata, Phys.Rev.E 54(1996)3687. 11. A.Yoshimori, Phys.Rev.E 59(1999)6535. 12. M.Cardenas,S.Franz and G.Parisi, J.Phys. A 31(1998)163.