Editing Computational implementation of integral equations

Integral equations are solved numerically.
One has the [[Ornstein-Zernike relation]], <math>\gamma (12)</math>
and a closure relation, <math>c_2 (12)</math> (which
incorporates the [[bridge function]] <math>B(12)</math>).
The numerical solution is iterative; 

# trial solution for  <math>\gamma (12)</math>
# calculate  <math>c_2 (12)</math>
# use the [[Ornstein-Zernike relation]] to generate a new  <math>\gamma (12)</math> ''etc.''

Note that the value of  <math>c_2 (12)</math> is '''local''', ''i.e.''
the  value of  <math>c_2 (12)</math> at a given point is given by
the value of  <math>\gamma (12)</math> 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 |(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:
===1. Closure relation <math>\gamma_{mns}^{\mu \nu} (r) \rightarrow c_{mns}^{\mu \nu} (r)</math>===
(Note: for linear fluids <math>\mu = \nu =0</math>)

====Perform the summation====

:<math>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)</math>

where <math>r_{12}</math> is the separation between molecular centers and 
<math>\omega_1,\omega_2</math> the sets of [[Euler angles]] needed to specify the orientations of the two molecules, with

:<math>\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)</math>

with <math>\overline{s} = -s</math>.

====Define the variables====

:<math>\left. x_1 \right.= \cos \theta_1</math>
:<math>\left. x_2\right.= \cos \theta_2</math>
:<math>\left. z_1 \right.= \cos \chi_1</math>
:<math>\left. z_2 \right.= \cos \chi_2</math>
:<math>\left. y\right.= \cos \phi_{12}</math>

Thus 
:<math>\gamma(12)=\gamma (r,x_1x_2,y,z_1z_2)</math>.

====Evaluate====
Evaluations of  <math>\gamma (12)</math> are performed at the discrete points <math>x_{i_1}x_{i_2},y_j,z_{k_1}z_{k_2}</math>
where the <math>x_i</math> are the <math>\nu</math> roots of the [[Legendre polynomial]] <math>P_\nu(cos \theta)</math>
where <math>y_j</math> are the  <math>\nu</math> roots of the [[Chebyshev polynomial]] <math>T_{\nu}(\ cos \phi)</math>
and where <math>z_{1_k},z_{2_k}</math>  are the  <math>\nu</math> roots of the [[Chebyshev polynomial]]
<math>T_{\nu}(\ cos \chi)</math>
thus

:<math>\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})</math>

where

:<math>\hat{d}_{s \mu}^m (x) = (2m+1)^{1/2} d_{s \mu}^m(\theta)</math>

where <math>d_{s \mu}^m(\theta)</math> is the angular, <math>\theta</math>, part of the
rotation matrix  <math>\mathcal{D}_{s \mu}^m (\omega)</math>,
and

:<math>e_s(y)=\exp(is\phi)</math>


:<math>e_{\mu}(z)= \exp(i\mu \chi)</math>

For the limits in the summations

:<math>L_1= \max (s,\nu_1)</math>

:<math>L_2= \max (s,\nu_2)</math>

The above equation constitutes a separable five-dimensional transform. To rapidly evaluate
this expression it is broken down into five one-dimensional transforms:

:<math>\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})</math>

:<math>\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})</math>

:<math>\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)</math>

:<math>\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})</math>

:<math>\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})</math>

Operations involving the <math>e_m(y)</math> and <math>e_n(z)</math> 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".
<math>NG</math> and <math>M</math> are parameters; <math>NG</math> is the number of nodes in the Gauss integration, and <math>M</math> the the max index in the truncated rotational invariants expansion.

====Integrate over angles <math>c_2(12)</math>====

Use [[Gauss-Legendre quadrature]] for <math>x_1</math> and <math>x_2</math>
Use [[Gauss-Chebyshev  quadrature]] for <math>y</math>, <math>z_1</math> and <math>z_2</math>
thus

:<math>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})</math>

where the Gauss-Legendre quadrature weights are given by

:<math>w_i= \frac{1}{(1-x_i^2)}[P_{NG}^{'} (x_i)]^2</math>

while the  Gauss-Chebyshev  quadrature has the constant weight

:<math>w=\frac{1}{NG}</math>

===Perform FFT from Real to Fourier space<math>c_{mns}^{\mu \nu} (r) \rightarrow   \tilde{c}_{mns}^{\mu \nu} (k)</math>====

This is non-trivial and is undertaken in three steps:
#Conversion from axial reference frame to spatial reference frame, ''i.e.''

:<math>c_{mns}^{\mu \nu} (r)  \rightarrow   c_{\mu \nu}^{mnl} (r)</math>

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

:<math>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)</math>

#'''Fourier-Bessel Transforms''': <math>c_{\mu \nu}^{mnl} (r) \rightarrow \tilde{c}_{\mu \nu}^{mnl} (k)</math>

:<math>\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</math>

(see Blum and Torruella Eq. 5.6 \cite{JCP_1972_56_00303} or Lado Eq. 39 \cite{MP_1982_47_0283}),
where <math>J_l(x)</math> is a [[Bessel function]] of order <math>l</math>.
`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.

<math>g(q)=2\pi \int_0^\infty f(r) J_0(2 \pi qr)r ~{\rm d}r</math>


<math>f(r)=2\pi \int_0^\infty g(q) J_0(2 \pi qr)q ~{\rm d}q</math>


#Conversion from the spatial reference frame back to the  axial reference frame
''i.e.''

<math>\tilde{c}_{\mu \nu}^{mnl} (k)  \rightarrow   \tilde{c}_{mns}^{\mu \nu} (k) 
</math>
this is done using the Blum transformation

<math>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)</math> 

OZ Equation} $  \tilde{c}_{mns}^{\mu \nu} (k)  \rightarrow   \tilde{\gamma}_{mns}^{\mu \nu} (k)$\\
~\\
For simple fluids: 
\begin{equation}
\tilde{\gamma}(k)= \frac{\rho \tilde{c}_2 (k)^2}{1- \rho  \tilde{c}_2 (k)}
\end{equation}
For molecular fluids (see Eq. 19 of Lado \cite{MP_1982_47_0283})
%(see derivation in the thesis of Juan Antonio Anta pp. 105--107):
%\begin{equation}
%\tilde{\Gamma}_{\chi}(k) = (-1)^{\chi}\rho \left[{\bf I} - (-1)^{\chi} \rho \tilde{\bf C}_{\chi}(k) \right]^{-1} \tilde{\bf C}_{\chi}(k)\tilde{\bf C}_{\chi}(k)
%\end{equation}
\begin{equation}
\tilde{{\bf S}}_{m}(k) = (-1)^{m}\rho \left[{\bf I} - (-1)^{m} \rho \tilde{\bf C}_{m}(k) \right]^{-1} \tilde{\bf C}_{m}(k)\tilde{\bf C}_{m}(k)
\end{equation}
where $\tilde{{\bf S}}_{m}(k)$ and $\tilde{\bf 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 \cite{JCP_1988_88_07715} and the thesis of Juan Antonio Anta pp. 107--109):
\begin{equation}
\tilde{\Gamma}(k) =  {\bf D}  \left[{\bf I} -  {\bf D}  \tilde{\bf C}(k)\right]^{-1} \tilde{\bf C}(k)\tilde{\bf C}(k)
\end{equation}
~\\
4) {\bf 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).\\
i) axial reference frame to spatial reference frame: $ \tilde{\gamma}_{mns}^{\mu \nu} (k) \rightarrow  \tilde{\gamma}^{mnl}_{\mu \nu} (k)$\\
ii) 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  \cite{MP_1982_47_0283}.\\
The inverse Hankel transform is
\begin{equation}
\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
\end{equation}
iii) Change from  spatial reference frame back to  axial reference frame:$   \gamma^{mnl}_{\mu \nu} (r) \rightarrow   \gamma_{mns}^{\mu \nu} (r)$.

==Ng acceleration==

*[http://dx.doi.org/10.1063/1.1682399  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)]
\section{Angular momentum coupling coefficients}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\cite{CPC_1970_1_0337,CPC_1971_2_0381}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Clebsch-Gordon coefficients and Racah's formula}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The Clebsch-Gordon coefficients are defined by
\begin{equation}
\Psi_{JM}= \sum_{M=M_1 + M_2} C_{M_1 M_2}^J \Psi_{M_1 M_2},
\end{equation}
where $J \equiv J_1 + J_2$ and satisfies
\begin{equation}
(j_1j_2m_1m_2|j_1j_2m)=0
\end{equation}
for $m_1+m_2\neq 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 (Giulio Racah (1909 - 1965)),
\begin{equation}
V(j_1j_2j;m_1m_2m)
\end{equation}
(See also the Racah W-coefficients, sometimes simply called the Racah coefficients).
\cite{CPC_1974_8_0095}

==References==
#[http://dx.doi.org/10.1080/00268977900102861 M. J. Gillan "A new method of solving the liquid structure integral equations" Molecular Physics '''38''' pp. 1781-1794 (1979)]
#[http://dx.doi.org/10.1080/00268978500102651 Stanislav Labík,  Anatol Malijevský and Petr Voncaronka "A rapidly convergent method of solving the OZ equation", Molecular Physics '''56''' pp. 709-715 (1985)]
#[http://dx.doi.org/10.1080/00268978200100202 F. Lado "Integral equations for fluids of linear molecules I. General formulation", Molecular Physics '''47''' pp. 283-298 (1982)]
#[http://dx.doi.org/10.1080/00268978200100212 F. Lado "Integral equations for fluids of linear molecules II. Hard dumbell solutions", Molecular Physics '''47''' pp. 299-311 (1982)]
#[http://dx.doi.org/10.1080/00268978200100222 F. Lado "Integral equations for fluids of linear molecules III. Orientational ordering", Molecular Physics '''47''' pp. 313-317 (1982)]
#[http://dx.doi.org/10.1080/00268978900101981 Enrique Lomba "An efficient procedure for solving the reference hypernetted chain equation (RHNC) for simple fluids" Molecular Physics '''68''' pp. 87-95 (1989)]