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;

1. trial solution for $\gamma (12)$
2. calculate $c_{2}(12)$
3. 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).

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:

Closure relation $\gamma _{mns}^{\mu \nu }(r)\rightarrow c_{mns}^{\mu \nu }(r)$

(Note: for linear fluids $\mu =\nu =0$)

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$.

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_{1}x_{2},y,z_{1}z_{2})$.

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_{2}m}^{n_{1}n_{2}}(r,x_{1_{i}})=\sum _{l_{1}=L_{1}}^{M}\gamma _{l_{1}l_{2}m}^{n_{1}n_{2}}(r){\hat {d}}_{mn_{1}}^{l_{1}}(x_{1_{i}})$

$\gamma _{m}^{n_{1}n_{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_{1}n_{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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle e_m(y)} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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". Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle NG} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle M} are parameters; $NG$ is the number of nodes in the Gauss integration, and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle M} the the max index in the truncated rotational invariants expansion.

Integrate over angles Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle c_2(12)}

Use Gauss-Legendre quadrature for Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x_1} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x_2} Use Gauss-Chebyshev quadrature for $y$, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle z_1} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle z_2} . Thus

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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}}$

Perform FFT from Real to Fourier space Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle c_{mns}^{\mu \nu} (r) \rightarrow \tilde{c}_{mns}^{\mu \nu} (k)}

This is non-trivial and is undertaken in three steps:

Conversion from axial reference frame to spatial reference frame

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle c_{mns}^{\mu \nu} (r) \rightarrow c_{\mu \nu}^{mnl} (r)}

this is done using the Blum transformation (Refs 7, 8 and 9):

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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)}

Fourier-Bessel Transforms

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle c_{\mu \nu}^{mnl} (r) \rightarrow \tilde{c}_{\mu \nu}^{mnl} (k)}

${\tilde {c}}_{\mu \nu }^{mnl}(k;l_{1}l_{2}ln_{1}n_{2})=4\pi i^{l}\int _{0}^{\infty }c_{\mu \nu }^{mnl}(r;l_{1}l_{2}ln_{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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle J_l(x)} is a Bessel function of order Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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.

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle g(q)=2\pi \int_0^\infty f(r) J_0(2 \pi qr)r ~{\rm d}r}

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle f(r)=2\pi \int_0^\infty g(q) J_0(2 \pi qr)q ~{\rm d}q}

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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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)}

Ornstein-Zernike relation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \tilde{c}_{mns}^{\mu \nu} (k) \rightarrow \tilde{\gamma}_{mns}^{\mu \nu} (k)}

For simple fluids:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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)

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \tilde{{\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 {\mathbf {S} }}_{m}(k)$ and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \tilde{\mathbf C}_{m}(k)} are matrices with elements Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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):

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \tilde{\Gamma}(k) = {\mathbf D} \left[{\mathbf I} - {\mathbf D} \tilde{\mathbf C}(k)\right]^{-1} \tilde{\mathbf C}(k)\tilde{\mathbf C}(k)}

Conversion back from Fourier space to Real space

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \tilde{\gamma}_{mns}^{\mu \nu} (k) \rightarrow \gamma_{mns}^{\mu \nu} (r)}

(basically the inverse of step 2).

Axial reference frame to spatial reference frame

${\tilde {\gamma }}_{mns}^{\mu \nu }(k)\rightarrow {\tilde {\gamma }}_{\mu \nu }^{mnl}(k)$

Inverse Fourier-Bessel transform

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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}

Change from spatial reference frame back to axial reference frame

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \gamma^{mnl}_{\mu \nu} (r) \rightarrow \gamma_{mns}^{\mu \nu} (r)} .

Ng acceleration

• Kin-Chue Ng "Hypernetted chain solutions for the classical one-component plasma up to Gamma=7000", Journal of Chemical Physics 61 pp. 2680-2689 (1974)

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)

References

1. M. J. Gillan "A new method of solving the liquid structure integral equations" Molecular Physics 38 pp. 1781-1794 (1979)
2. Stanislav Labík, Anatol Malijevský and Petr Voncaronka "A rapidly convergent method of solving the OZ equation", Molecular Physics 56 pp. 709-715 (1985)
3. F. Lado "Integral equations for fluids of linear molecules I. General formulation", Molecular Physics 47 pp. 283-298 (1982)
4. F. Lado "Integral equations for fluids of linear molecules II. Hard dumbell solutions", Molecular Physics 47 pp. 299-311 (1982)
5. F. Lado "Integral equations for fluids of linear molecules III. Orientational ordering", Molecular Physics 47 pp. 313-317 (1982)
6. Enrique Lomba "An efficient procedure for solving the reference hypernetted chain equation (RHNC) for simple fluids" Molecular Physics 68 pp. 87-95 (1989)
7. 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)
8. 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)
9. 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)
10. 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)