Journal of Molecular Liquids 120 (2005) 159 – 162 www.elsevier.com/locate/molliq
Wavelet method for solving integral equations of simple liquids Maxim V. Fedorova, Gennady N. Chuevb,* a
Theory and Computation Group, Centre for Synthesis and Chemical Biology, Conway Institute of Biomolecular and Biomedical Research, Department of Chemistry, University College Dublin, Belfield, Dublin 4, Ireland b Institute of Theoretical and Experimental Biophysics, Russian Academy of Science, Pushchino 142290, Russia Available online 21 August 2004
Abstract A new algorithm is developed to solve integral equations for simple liquids. The algorithm is based on the discrete wavelet transform of radial distribution functions. The Coifman 2 basis set is employed for the wavelet treatment. Using the algorithm, we have calculated structural and thermodynamic properties of a Lennard–Jones fluid in a wide range of energy and size parameters of the fluid. D 2004 Published by Elsevier B.V. Keywords: Lennard–Jones fluid; Integral equation; Wavelets
1. Introduction
2. Theory
Wavelets are basic elements of multiresolution analysis [1,2] for treating multiscale processes, which provided simultaneously a homogeneous underlying description of processes, and the capacity to control and vary the resolution of this description at will. Recently we have used discrete wavelets to parameterize the radial distribution functions of hydrated ions and hydrophobic solutes [3]. Our study indicates that wavelets are a suitable and powerful instrument for approximating distribution functions in the classical domain. The goal of this work is to apply discrete wavelets for calculations of thermodynamic properties and structure of simple liquids. Simple liquids are very attractive for application of the integral equation (IE) theory due to analytical solutions obtained in some cases. Thus, we will use discrete wavelets to calculate radial distribution function (RDF) of simple liquids within the framework of the IE theory.
The IE theory has proved to be a successful tool for treating simple liquids. From a theoretical point of view, the IE method consists in obtaining the RDF by solving the set of equations formed by the Ornstein–Zernike (OZ) integral equation and a closure. The OZ equation couples the total correlation function. h(r)=g(r)1 and the direct correlation function c(r). For an isotropic liquid with density n˜, the OZ equation is given by:
* Corresponding author. E-mail addresses:
[email protected] (M.V. Fedorov)8
[email protected] (G.N. Chuev). 0167-7322/$ - see front matter D 2004 Published by Elsevier B.V. doi:10.1016/j.molliq.2004.07.060
hðrÞ ¼ cðrÞ þ qcðjR rjÞ4hðRÞ;
ð1Þ
where the symbol * means the convolution in the real space. The closure relation couples the same quantities and the interacting potential U(r) and can be written in terms of the indirect correlation function c(r)=h(r)c(r) as: hðrÞ ¼ exp½ bU ðrÞ þ cðrÞ þ B½cðrÞ 1;
ð2Þ
where b=(k BT)1 is the inverse temperature. The closure relation introduces a third function, namely, the bridge function B[c(r)]. Since there is no exact expression for the dependence B[c(r)], the approximation to the bridge function is the key to contemporary IE theories. The list of such approximating closures is still expanding and is well documented in the literature [4]. For the numerical
160
M.V. Fedorov, G.N. Chuev / Journal of Molecular Liquids 120 (2005) 159–162
calculations, the Fourier representation of the OZ equation is more useful: cðk Þ ¼ q2 cðk Þ=ð1 qcðk ÞÞ;
ð3Þ
where c(k) and c(k) are the three-dimensional Fourier transforms (FT) of the functions c(r) and c(r). Current versions of the numerical scheme for solving IE suggest the solution to be divided into the coarse and fine parts and the algorithm to be a hybrid of the Newton–Raphson and the Picard schemes for the coarse and fine parts of the solution, respectively [5,6]. We restrict our discussion to orthonormal compact supported wavelets [1,2,7,8]. Like the plane waves in Fourier analysis, the wavelets are the basis for expanding square-integrable function in L 2(R). But unlike the harmonic functions, wavelets have dual localization characteristics in both real and reciprocal spaces. We use the discrete wavelet transform (DWT) [1,9] to parameterize RDF. Any square integrable function f(r) can be expanded as a wavelet series [9]: f ðr Þ ¼
X s
aj0s uj0s ðrÞ þ
l X X j0
djs wjs ðrÞ;
ð4Þ
transform (DWT). In the expansion (4), the first term gives a coarse approximation of RDF at the resolution j 0 and the second term gives a sequence of successive details. In practice, we actually do not use the infinite resolution therefore the sequences are cut off at a resolution j max, {s i } has also a finite number of terms S. In numerical calculation of the coefficients (4), we can avoid the direct integration using the algorithm of fast wavelet transform (FWT) [1]. The elegant pyramidal procedure can be performed for this purpose. We represent c(r) in terms of choosing wavelet basis functions. Let us denote the complete set of the coefficients {a j0s } and {d js } as the matrices a and d, respectively. Thus, with the required accuracy we can substitute the c(r) in the integral equation by the obtained approximation c ap(a,d). For the given c ap(a,d), we can calculate direct correlation function c ap(r) via (2) and then find c(k) as a dependence on the matrices (a,d) using direct FT and Eq. (3), after that the next approximation c(a,d) is obtained via inverse FT. Then new a and d are obtained by DWT. Finally our numerical scheme proposes the following cycle: ain ; din Yð4ÞYcap ða; dÞYð2ÞYcap ðrÞYðFT ÞYcðk ÞYð3Þ
s
where u(r) is the scaling function and w(r) is the mother wavelet of a discrete wavelet basis set {u(r),w(r)}, j, s are algebraic integers and j 0 is a chosen level of resolution. The functions u j0s (r) and w js (r) are generated with the use of pair of functions {u(r),w(r)} by integer translation and stretching (or squeezing), s is the translation parameter, while j is the scaling factors indicating the wavelet stretching (squeezing). Eq. (1) provides the discrete wavelet
Ycðk ÞY N ðFTÞ1 YcðrÞYðDWTÞYaout ; dout
ð5Þ
The efficiency of the method depends on the number of approximating coefficients. Due to the special choice of the wavelet basis, this number can be rather small, providing a rapid convergence of iterations. To illustrate our scheme, we have investigated the uniform Lennard–Jones fluid with the interacting potential U(r) given by U(r)=4e((r/R)12(r/R)6), where j and q are
Fig. 1. Difference in the total correlation functions obtained by the wavelets and by the conventional IE method for qr 3=0.8 and be=0.5. The insets shows the RDFs near the first peak, which are obtained by wavelets (solid line) and by the conventional IE method (dashed line).
M.V. Fedorov, G.N. Chuev / Journal of Molecular Liquids 120 (2005) 159–162
161
Fig. 2. Amplitudes of the coefficients of the wavelet expansion for the solutions found by wavelets (solid line) and by the conventional IE method.
Table 1 Thermodynamic parameters calculated by the wavelets with HNC closure (in parentheses are the same quantities obtained by the conventional method) qj3
be
E
Z
N
v
0.2
0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0
0.41 (0.40) 0.94 (0.98) 1.661 (1.66) 2.42 0.77 (0.76) 1.91 (1.91) 3.10 (3.15) 4.27 (4.61) 0.95 (1.0) 2.82 (2.74) 4.56 (4.55) 6.43 (6.42) 8.52 (8.38) 1.40 (0.95) 3.35 (3.2) 5.60 (5.58) 8.07 (8.04) 10.64 (10.54) 0.54 (0.37) 2.91 (2.81) 5.34 (5.51) 8.15 (8.35) 11.27 (11.27)
1.17 (1.18) 1.08 (0.96) 0.93 (0.71) 0.84 1.58 (1.61) 1.22 (1.21) 0.86 (0.74) 1.06 (2.78) 2.65 (2.55) 1.93 (2.16) 1.54 (1.56) 0.93 (0.88) 0.40 (0.18) 3.21 (4.39) 4.22 (4.54) 4.19 (4.26) 3.65 (3.75) 3.47 (3.11) 9.71 (7.67) 9.54 (9.40) 10.93 (10.40) 11.31 (11.04) 11.43 (11.43)
4.39 (4.48) 4.71 (4.77) 4.57 (5.03) 5.35 7.51 (7.54) 7.81 (7.84) 7.66 (7.97) 8.24 (8.22) 9.57 (9.69) 10.20 (9.98) 9.90 (10.08) 9.95 (10.14) 10.41 (10.16) 13.45 (11.0) 12.02 (11.26) 11.09 (11.33) 11.08 (11.38) 11.95 (11.4) 11.86 (11.71) 12.71 (11.88) 11.74 (11.91) 11.44 (11.92) 11.73 (11.94)
1.46 (1.41) 1.07 (0.98) 0.63 (0.49) 0.15 2.56 (2.47) 1.88 (1.70) 1.05 (0.77) 0.22 (0.14) 4.82 (4.67) 4.27 (3.99) 3.26 (2.84) 2.05 (1.50) 0.78 (0.13) 9.02 (8.79) 9.74 (9.28) 8.79 (7.85) 8.79 (7.85) 7.76 (6.60) 16.21 (15.8) 20.41 (19.52) 23.09 (21.65) 24.98 (22.96) 26.36 (23.76)
0.4
0.6
0.8
1.0
162
M.V. Fedorov, G.N. Chuev / Journal of Molecular Liquids 120 (2005) 159–162
the size and energy parameters, respectively. We have used a grid with a number of points n=2048 and step size dr=0.025 2. The integral equations were solved iteratively as shown in the previous section until the required accuracy is achieved. The precision parameter for the numerical solution was equal to 104. Since the main criterion of the IE theory is concerned with calculations of various thermodynamic quantities of LJ liquids, we investigated the influence of the wavelet set not only on the quality of the RDFs, but also on calculations of the internal energy E, the compressibility factor Z, the coordination number N, and the isothermal compressibility v using common equations for these parameters [4]. We have applied the Coifman 2 wavelet set. We have used the improved Picard method to calculate the matrices a and d. The method converged rather fast even for high densities. For low and middle densities, we took the hard-sphere approximations as initial estimates. We took the middle density solutions of LJ fluid as initial data for ain,din in the case of high densities.
3. Results and discussion Using the IE method based on wavelets, we have calculated the RDF and thermodynamic properties of the LJ fluid. The properties of the LJ were evaluated in a wide range of the size and energy parameters of the fluid. For comparison, we have calculated all these quantities by the conventional IE method based on the hybrid scheme [6]. We should mention that the computational cost for our scheme is less than 100 times that of the conventional one. Fig. 1 depicts the differences dc(r)=c ap(r)c IE(r) and dh(r)=h ap(r)h IE(r) where c ap(r) and h ap(r) refer to the functions obtained by the wavelet method, while c IE(r) and h IE(r) correspond to the calculations by the conventional IE method. As it is seen, the wavelet approximation is rather close to that obtained by the conventional method, a relatively pronounced difference takes place only at the region of the first peak; however, due to high values of the RDF in this region, the above difference does not influence the accuracy of the calculations (see insets in Fig. 1). To reveal the properties of the wavelet approximation, we compare the amplitude of the wavelet coefficients obtained by our method and that extracted from the conventional IE solutions expanded in the same wavelet series (Fig. 2). As is seen, the amplitude of the coefficients rapidly decreases as the resolution rises. Because of the method of the wavelet approximation, the high-level coefficients of details equal zero, while the exact expansion of the IE solutions results in the infinite series. On the other hand, the difference between the wavelet approximation and the exact approximation is small even at the coarse resolution. This clearly indicates that our scheme yields correct results. The comparison with the thermodynamic quantities obtained by the conventional method indicates a good agreement between the data (Table 1).
In itself, the fact of the converging of the simple Picard method for reasonable number of iterations is remarkable. It could be considered as a validation of high stability of our scheme. Indeed, the current opinion about the Picard method is its unfitness [5,6] due to the slow and unreliable convergence. The method is usually presented as divergent in the case of high densities. The base of high stability is the lucky choice of the basis sets. The solving of the IE in a finite grid demands reliable representation of the correlation functions in both real and reciprocal spaces. At that, the number of significant coefficients should be as small as possible. The orthogonal discrete wavelets have finite support in both spaces that is the main reason of their amazing approximating properties. The hybrid algorithms [5,6] use the basis sets well localized only in one space. The roof functions in the Gillan scheme [5] have finite support in the real space but are spread in the reciprocal one that leads to the numerical artifacts in the Fourier representation of approximated functions. On the contrary, the plane waves [6] are well localized in the k-space but infinite-dimensional in the real one. In the case of finite grid, it leads to the numerical boundary artifacts, i.e. false pulsations of the approximated functions in the real space. The applied scheme of wavelet decomposition allows us to treat OZ equation with a small number of approximating coefficients. Due to this fact, we hope the wavelets can be applied not only to approximate RDF but also to calculate them directly by the density functional methods. The application of DFT to molecular liquids [10] indicates that the method is very promising to treat solvation phenomena. Discrete wavelets can be applied to calculate distribution functions of quantum-classical system treated by the density functional method. Recently, we have used the Coifman basis set for approximating the correlation functions and consequent minimization of the free energy functional of the two coupled electrons in polar liquid [11]. We suppose that the wavelet study can provide further progress in this direction. References [1] S.G. Mallat, A Wavelet Tour of Signal Processing, 2nd ed., Academic Press, San Diego, 1999. [2] Y. Meyer, Wavelets: Algorithms and Applications, SIAM, Philadelphia, 1993. [3] G.N. Chuev, M.V. Fedorov, Phys. Rev., E. 68 (2003) 1 (027702). [4] J.-P. Hansen, I.R. McDonald, Theory of Simple Liquids, Academic, London, 1986. [5] M.J. Gillan, Mol. Phys. 38 (1979) 1781. [6] S. Labik, A. Malijevsky, P. Vonka, Mol. Phys. 56 (1985) 709. [7] D.L. Donoho, Appl. Comput. Harmon. Anal. 1 (1993) 1. [8] R.T. Ogden, Essential Wavelets for Statistical Applications and Data Analysis, Birkhauser, Boston, 1997. [9] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, 1992. [10] D. Chandler, J.D. McCoy, S.J. Singer, J. Chem. Phys. 85 (1986) 5971. [11] G.N. Chuev, V. Fedorov, Int. J. Quantum Chem. 99 (2004) (in press).