Computational implementation of integral equations
Integral equations are solved numerically. One has the Ornstein-Zernike relation, \(\gamma (12)\) and a closure relation, \(c_2 (12)\) (which incorporates the bridge function \(B(12)\)). The numerical solution is iterative;
- trial solution for \(\gamma (12)\)
- calculate \(c_2 (12)\)
- use the Ornstein-Zernike relation to generate a new \(\gamma (12)\) etc.
Note that the value of \(c_2 (12)\) is local, i.e. the value of \(c_2 (12)\) at a given point is given by the value of \(\gamma (12)\) at this point. However, the Ornstein-Zernike relation is non-local. The way to convert the Ornstein-Zernike relation into a local equation is to perform a (fast) Fourier transform (FFT). Note: convergence is poor for liquid densities. (See Ref.s 1 to 6).
[edit] Picard iteration
Picard iteration generates a solution of an initial value problem for an ordinary differential equation (ODE) using fixed-point iteration. Here are the four steps used to solve integral equations:
[edit] Closure relation \(\gamma_{mns}^{\mu \nu} (r) \rightarrow c_{mns}^{\mu \nu} (r)\)
(Note: for linear fluids \(\mu = \nu =0\))
[edit] Perform the summation
\[g(12)=g(r_{12},\omega_1,\omega_2)=\sum_{mns\mu \nu} g_{mns}^{\mu \nu}(r_{12}) \Psi_{\mu \nu s}^{mn}(\omega_1,\omega_2)\]
where \(r_{12}\) is the separation between molecular centers and \(\omega_1,\omega_2\) the sets of Euler angles needed to specify the orientations of the two molecules, with
\[\Psi_{\mu \nu s}^{mn}(\omega_1,\omega_2) = \sqrt{(2m+1)(2n+1)} \mathcal{D}_{s \mu}^m (\omega_1) \mathcal{D}_{\overline{s} \nu}^n (\omega_2)\]
with \(\overline{s} = -s\).
[edit] Define the variables
\[\left. x_1 \right.= \cos \theta_1\] \[\left. x_2\right.= \cos \theta_2\] \[\left. z_1 \right.= \cos \chi_1\] \[\left. z_2 \right.= \cos \chi_2\] \[\left. y\right.= \cos \phi_{12}\]
Thus \[\left. \gamma(12) \right. =\gamma (r,x_1x_2,y,z_1z_2)\].
[edit] Evaluate
Evaluations of \(\gamma (12)\) are performed at the discrete points \(x_{i_1}x_{i_2},y_j,z_{k_1}z_{k_2}\) where the \(x_i\) are the \(\nu\) roots of the Legendre polynomial \(P_\nu(cos \theta)\) where \(y_j\) are the \(\nu\) roots of the Chebyshev polynomial \(T_{\nu}(\ cos \phi)\) and where \(z_{1_k},z_{2_k}\) are the \(\nu\) roots of the Chebyshev polynomial \(T_{\nu}(\ cos \chi)\) thus
\[\gamma(r,x_{1_i},x_{2_i},j,z_{1_k},z_{2_k})= \sum_{\nu , \mu , s = -M }^M \sum_{m=L_2}^M \sum_{n=L_1}^M \gamma_{mns}^{\mu \nu} (r) \hat{d}_{s \mu}^m (x_{1_i}) \hat{d}_{\overline{s} \nu}^n (x_{2_i}) e_s(j) e_{\mu} (z_{1_k}) e_{\nu} (z_{2_k})\]
where
\[\hat{d}_{s \mu}^m (x) = (2m+1)^{1/2} d_{s \mu}^m(\theta)\]
where \(d_{s \mu}^m(\theta)\) is the angular, \(\theta\), part of the
rotation matrix \(\mathcal{D}_{s \mu}^m (\omega)\),
and
\[\left. e_s(y) \right.=\exp(is\phi)\]
\[\left. e_{\mu}(z) \right.= \exp(i\mu \chi)\]
For the limits in the summations
\[\left. L_1 \right.= \max (s,\nu_1)\]
\[\left. L_2 \right.= \max (s,\nu_2)\]
The above equation constitutes a separable five-dimensional transform. To rapidly evaluate this expression it is broken down into five one-dimensional transforms:
\[\gamma_{l_2m}^{n_1n_2}(r,x_{1_i})=\sum_{l_1=L_1}^M \gamma_{l_1 l_2 m}^{n_1 n_2}(r) \hat{d}_{m n_1}^{l_1} (x_{1_i})\]
\[\gamma_{m}^{n_1n_2}(r,x_{1_i},x_{2_i})=\sum_{l_2=L_2}^M \gamma_{l_2 m}^{n_1 n_2}(r,x_{1_i}) \hat{d}_{\overline{m} n_2}^{l_2} (x_{2_i})\]
\[\gamma^{n_1n_2}(r,x_{1_i},x_{2_i},j)=\sum_{m=-M}^M \gamma_{m}^{n_1 n_2}(r,x_{1_i},x_{2_i}) e_m(j)\]
\[\gamma^{n_2}(r,x_{1_i},x_{2_i},z_{1_k})=\sum_{n_1=-M}^M \gamma^{n_1 n_2}(r,x_{1_i},x_{2_i},j) e_{n_1}(z_{1_k})\]
\[\gamma(r,x_{1_i},x_{2_i},z_{1_k},z_{2_k})=\sum_{n_2=-M}^M \gamma^{n_2}(r,x_{1_i},x_{2_i},j,z_{1_k}) e_{n_2}(z_{2_k})\]
Operations involving the \(e_m(y)\) and \(e_n(z)\) basis functions are performed in complex arithmetic. The sum of these operations is asymptotically smaller than the previous expression and thus constitutes a ``fast separable transform". \(NG\) and \(M\) are parameters; \(NG\) is the number of nodes in the Gauss integration, and \(M\) the the max index in the truncated rotational invariants expansion.
[edit] Integrate over angles \(c_2(12)\)
Use Gauss-Legendre quadrature for \(x_1\) and \(x_2\) Use Gauss-Chebyshev quadrature for \(y\), \(z_1\) and \(z_2\). Thus
\[c_{mns}^{\mu \nu} (r) = w^3 \sum_{x_{1_i},x_{2_i},j,z_{1_k},z_{2_k}=1}^{NG} w_{i_1}w_{i_2}c_2(r,x_{1_i},x_{2_i},j,z_{1_k},z_{2_k}) \hat{d}_{s \mu}^m (x_{1_i}) \hat{d}_{\overline{s} \nu}^n (x_{2_i}) e_{\overline{s}}(j) e_{\overline{\mu}} (z_{1_k}) e_{\overline{\nu}} (z_{2_k})\]
where the Gauss-Legendre quadrature weights are given by
\[w_i= \frac{1}{(1-x_i^2)}[P_{NG}^{'} (x_i)]^2\]
while the Gauss-Chebyshev quadrature has the constant weight
\[w=\frac{1}{NG}\]
[edit] Perform FFT from Real to Fourier space \(c_{mns}^{\mu \nu} (r) \rightarrow \tilde{c}_{mns}^{\mu \nu} (k)\)
This is non-trivial and is undertaken in three steps:
[edit] Conversion from axial reference frame to spatial reference frame
\[c_{mns}^{\mu \nu} (r) \rightarrow c_{\mu \nu}^{mnl} (r)\]
this is done using the Blum transformation (Refs 7, 8 and 9):
\[g_{\mu \nu}^{mnl}(r) = \sum_{s=-\min (m,n)}^{\min (m,n)} \left( \begin{array}{ccc} m&n&l\\ s&\overline{s}&0 \end{array} \right)g_{mns}^{\mu \nu} (r)\]
[edit] Fourier-Bessel Transforms
\[c_{\mu \nu}^{mnl} (r) \rightarrow \tilde{c}_{\mu \nu}^{mnl} (k)\]
\[\tilde{c}_{\mu \nu}^{mnl} (k; l_1 l_2 l n_1 n_2) = 4\pi i^l \int_0^{\infty} c_{\mu \nu}^{mnl} (r; l_1 l_2 l n_1 n_2) J_l (kr) ~r^2 {\rm d}r\]
(see Blum and Torruella Eq. 5.6 in Ref. 7 or Lado Eq. 39 in Ref. 3), where \(J_l(x)\) is a Bessel function of order \(l\). `step-down' operations can be performed by way of sin and cos operations of Fourier transforms, see Eqs. 49a, 49b, 50 of Lado Ref. 3. The Fourier-Bessel transform is also known as a Hankel transform. It is equivalent to a two-dimensional Fourier transform with a radially symmetric integral kernel.
\[g(q)=2\pi \int_0^\infty f(r) J_0(2 \pi qr)r ~{\rm d}r\]
\[f(r)=2\pi \int_0^\infty g(q) J_0(2 \pi qr)q ~{\rm d}q\]
[edit] Conversion from the spatial reference frame back to the axial reference frame
\[\tilde{c}_{\mu \nu}^{mnl} (k) \rightarrow \tilde{c}_{mns}^{\mu \nu} (k) \] this is done using the Blum transformation
\[g_{mns}^{\mu \nu} (r) = \sum_{l=|m-n|}^{m+n} \left( \begin{array}{ccc} m&n&l\\ s&\overline{s}&0 \end{array} \right) g_{\mu \nu}^{mnl}(r)\]
[edit] Ornstein-Zernike relation \(\tilde{c}_{mns}^{\mu \nu} (k) \rightarrow \tilde{\gamma}_{mns}^{\mu \nu} (k)\)
For simple fluids:
\[\tilde{\gamma}(k)= \frac{\rho \tilde{c}_2 (k)^2}{1- \rho \tilde{c}_2 (k)}\]
For molecular fluids (see Eq. 19 of Lado Ref. 3)
\[\tilde[[:Template:\mathbf S]]_{m}(k) = (-1)^{m}\rho \left[{\mathbf I} - (-1)^{m} \rho \tilde{\mathbf C}_{m}(k) \right]^{-1} \tilde{\mathbf C}_{m}(k)\tilde{\mathbf C}_{m}(k)\]
where \(\tilde[[:Template:\mathbf S]]_{m}(k)\) and \(\tilde{\mathbf C}_{m}(k)\) are matrices with elements \(\tilde S_{l_1 l_2 m}(k), \tilde{C}_{l_1 l_2 m}(k), l_1,l_2 \geq m\).
For mixtures of simple fluids (see Ref. 10 Juan Antonio Anta PhD thesis pp. 107--109):
\[\tilde{\Gamma}(k) = {\mathbf D} \left[{\mathbf I} - {\mathbf D} \tilde{\mathbf C}(k)\right]^{-1} \tilde{\mathbf C}(k)\tilde{\mathbf C}(k)\]
[edit] Conversion back from Fourier space to Real space
\[\tilde{\gamma}_{mns}^{\mu \nu} (k) \rightarrow \gamma_{mns}^{\mu \nu} (r)\] (basically the inverse of step 2).
[edit] Axial reference frame to spatial reference frame
\[\tilde{\gamma}_{mns}^{\mu \nu} (k) \rightarrow \tilde{\gamma}^{mnl}_{\mu \nu} (k)\]
[edit] Inverse Fourier-Bessel transform
\[\tilde{\gamma}^{mnl}_{\mu \nu} (k) \rightarrow \gamma^{mnl}_{\mu \nu} (r)\] 'Step-up' operations are given by Eq. 53 of Ref. 3. The inverse Hankel transform is \[\gamma(r;l_1 l_2 l n_1 n_2)= \frac{1}{2 \pi^2 i^l} \int_0^\infty \tilde{\gamma}(k;l_1 l_2 l n_1 n_2) J_l (kr) ~k^2 {\rm d}k\]
[edit] Change from spatial reference frame back to axial reference frame
\[\gamma^{mnl}_{\mu \nu} (r) \rightarrow \gamma_{mns}^{\mu \nu} (r)\].
[edit] Ng acceleration
[edit] Angular momentum coupling coefficients
- Taro Tamura "Angular momentum coupling coefficients", Computer Physics Communications 1 pp. 337-342 (1970)
- J. G. Wills "On the evaluation of angular momentum coupling coefficients", omputer Physics Communications 2 pp. 381-382 (1971)
[edit] References
- M. J. Gillan "A new method of solving the liquid structure integral equations" Molecular Physics 38 pp. 1781-1794 (1979)
- Stanislav Labík, Anatol Malijevský and Petr Voncaronka "A rapidly convergent method of solving the OZ equation", Molecular Physics 56 pp. 709-715 (1985)
- F. Lado "Integral equations for fluids of linear molecules I. General formulation", Molecular Physics 47 pp. 283-298 (1982)
- F. Lado "Integral equations for fluids of linear molecules II. Hard dumbell solutions", Molecular Physics 47 pp. 299-311 (1982)
- F. Lado "Integral equations for fluids of linear molecules III. Orientational ordering", Molecular Physics 47 pp. 313-317 (1982)
- Enrique Lomba "An efficient procedure for solving the reference hypernetted chain equation (RHNC) for simple fluids" Molecular Physics 68 pp. 87-95 (1989)
- L. Blum and A. J. Torruella "Invariant Expansion for Two-Body Correlations: Thermodynamic Functions, Scattering, and the Ornstein—Zernike Equation", Journal of Chemical Physics 56 pp. pp. 303-310 (1972)
- L. Blum "Invariant Expansion. II. The Ornstein-Zernike Equation for Nonspherical Molecules and an Extended Solution to the Mean Spherical Model", Journal of Chemical Physics 57 pp. 1862-1869 (1972)
- L. Blum "Invariant expansion III: The general solution of the mean spherical model for neutral spheres with electostatic interactions", Journal of Chemical Physics 58 pp. 3295-3303 (1973)
- P. G. Kusalik and G. N. Patey " On the molecular theory of aqueous electrolyte solutions. I. The solution of the RHNC approximation for models at finite concentration", Journal of Chemical Physics 88 pp. 7715-7738 (1988)
