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 $\underset{mns}{\overset{\mu \nu }{\gamma }}(r)\to 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})\underset{\mu \nu s}{\overset{mn}{\mathrm{\Psi }}}({\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

$\underset{\mu \nu s}{\overset{mn}{\mathrm{\Psi }}}({\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

$x_1=\cos {\theta }_1$

$x_2=\cos {\theta }_2$

$z_1=\cos {\chi }_1$

$z_2=\cos {\chi }_2$

$y=\cos {\varphi }_{12}$

Thus

$\gamma (12)=\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\varphi )$ 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\underset{mns}{\overset{\mu \nu }{\gamma }}(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

$e_s(y)=\exp (is\varphi )$

$e_{\mu }(z)=\exp (i\mu \chi )$

For the limits in the summations

$L_1=max(s,{\nu }_1)$

$L_2=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:

$\underset{l_{2}m}{\overset{n_1{n}_2}{\gamma }}(r,x_{1_i})=\sum _{l_1=L_1}^M\underset{l_1{l}_{2}m}{\overset{n_1{n}_2}{\gamma }}(r){\hat{d}}_{m{n}_1}^{l_1}(x_{1_i})$

$\underset{m}{\overset{n_1{n}_2}{\gamma }}(r,x_{1_i},x_{2_i})=\sum _{l_2=L_2}^M\underset{l_{2}m}{\overset{n_1{n}_2}{\gamma }}(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\underset{m}{\overset{n_1{n}_2}{\gamma }}(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.

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}^{\text{'}}(x_i){]}^2$

while the Gauss-Chebyshev quadrature has the constant weight

$w=\frac{1}{NG}$

Perform FFT from Real to Fourier space $c_{mns}^{\mu \nu }(r)\to {\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

$c_{mns}^{\mu \nu }(r)\to c_{\mu \nu }^{mnl}(r)$

this is done using the Blum transformation \cite{JCP_1972_56_00303,JCP_1972_57_01862,JCP_1973_58_03295}:

$g_{\mu \nu }^{mnl}(r)=\sum _{s=-min(m,n)}^{min(m,n)}(\begin{array}{ccc}m& n& l\\ s& & 0\end{array})g_{mns}^{\mu \nu }(r)$

Fourier-Bessel Transforms

$c_{\mu \nu }^{mnl}(r)\to {\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\underset{0}{\overset{\mathrm{\infty }}{\int }}{c}_{\mu \nu }^{mnl}(r;l_1{l}_{2}l{n}_1{n}_2)J_l(kr)r^2\mathrm{d}r$

(see Blum and Torruella Eq. 5.6 \cite{JCP_1972_56_00303} or Lado Eq. 39 \cite{MP_1982_47_0283}), 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 \cite{MP_1982_47_0283}. 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 \underset{0}{\overset{\mathrm{\infty }}{\int }}f(r)J_0(2\pi qr)r\mathrm{d}r$

$f(r)=2\pi \underset{0}{\overset{\mathrm{\infty }}{\int }}g(q)J_0(2\pi qr)q\mathrm{d}q$

Conversion from the spatial reference frame back to the axial reference frame

${\tilde{c}}_{\mu \nu }^{mnl}(k)\to {\tilde{c}}_{mns}^{\mu \nu }(k)$

this is done using the Blum transformation

$g_{mns}^{\mu \nu }(r)=\sum _{l=|m-n|}^{m+n}(\begin{array}{ccc}m& n& l\\ s& & 0\end{array})g_{\mu \nu }^{mnl}(r)$

Ornstein-Zernike relation ${\tilde{c}}_{mns}^{\mu \nu }(k)\to {\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 \cite{MP_1982_47_0283})

${\tilde{S}}_m(k)=(-1{)}^m\rho {[I-(-1{)}^m\rho {\tilde{C}}_m(k)]}^{-1}{\tilde{C}}_m(k){\tilde{C}}_m(k)$

where ${\tilde{S}}_m(k)$ and ${\tilde{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\ge m$. For mixtures of simple fluids (see \cite{JCP_1988_88_07715} and the thesis of Juan Antonio Anta pp. 107--109):

$\tilde{\mathrm{\Gamma }}(k)=D{[I-D\tilde{C}(k)]}^{-1}\tilde{C}(k)\tilde{C}(k)$

Conversion back from Fourier space to Real space

${\tilde{\gamma }}_{mns}^{\mu \nu }(k)\to \underset{mns}{\overset{\mu \nu }{\gamma }}(r)$ (basically the inverse of step 2).

Axial reference frame to spatial reference frame

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

Inverse Fourier-Bessel transform

${\tilde{\gamma }}_{\mu \nu }^{mnl}(k)\to \underset{\mu \nu }{\overset{mnl}{\gamma }}(r)$ `Step-up' operations are given by Eq. 53 of \cite{MP_1982_47_0283}. The inverse Hankel transform is

$\gamma (r;l_1{l}_{2}l{n}_1{n}_2)=\frac{1}{2{\pi }^2{i}^l}\underset{0}{\overset{\mathrm{\infty }}{\int }}\tilde{\gamma }(k;l_1{l}_{2}l{n}_1{n}_2)J_l(kr)k^2\mathrm{d}k$

Change from spatial reference frame back to axial reference frame

$\underset{\mu \nu }{\overset{mnl}{\gamma }}(r)\to \underset{mns}{\overset{\mu \nu }{\gamma }}(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)

Clebsch-Gordon coefficients and Racah's formula

The Clebsch-Gordon coefficients are defined by

${\mathrm{\Psi }}_{JM}={\sum }_{M=M_1+M_2}{C}_{M_1{M}_2}^J{\mathrm{\Psi }}_{M_1{M}_2},$

where $J\equiv J_1+J_2$ and satisfies $(j_1{j}_2{m}_1{m}_2|j_1{j}_{2}m)=0$ for $m_1+m_2\ne m$. They are used to integrate products of three spherical harmonics (for example the addition of angular momenta). The Clebsch-Gordan coefficients are sometimes expressed using the related Racah V-coefficients, $V(j_1{j}_{2}j;m_1{m}_{2}m)$ (See also the Racah W-coefficients, sometimes simply called the Racah coefficients).

• Robert E. Beck and Bernard Kolman "Racah's outer multiplicity formula", Computer Physics Communications 8 pp. 95-100 (1974)

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)