Computational implementation of integral equations

From SklogWiki
Jump to: navigation, search

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[edit]

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)[edit]

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

Perform the summation[edit]

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[edit]

\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}


\left. \gamma(12) \right. =\gamma (r,x_1x_2,y,z_1z_2).


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

\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})


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

Integrate over angles c_2(12)[edit]

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 
\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


Perform FFT from Real to Fourier space c_{mns}^{\mu \nu} (r) \rightarrow   \tilde{c}_{mns}^{\mu \nu} (k)[edit]

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

Conversion from axial reference frame to spatial reference frame[edit]

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( 
\right)g_{mns}^{\mu \nu} (r)

Fourier-Bessel Transforms[edit]

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

Conversion from the spatial reference frame back to the axial reference frame[edit]

\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( 
g_{\mu \nu}^{mnl}(r)

Ornstein-Zernike relation \tilde{c}_{mns}^{\mu \nu} (k)  \rightarrow   \tilde{\gamma}_{mns}^{\mu \nu} (k)[edit]

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{{\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 \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)

Conversion back from Fourier space to Real space[edit]

\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[edit]

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

Inverse Fourier-Bessel transform[edit]

\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

Change from spatial reference frame back to axial reference frame[edit]

\gamma^{mnl}_{\mu \nu} (r) \rightarrow   \gamma_{mns}^{\mu \nu} (r).

Ng acceleration[edit]

Angular momentum coupling coefficients[edit]


  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)