2017: SklogWiki celebrates 10 years on-line

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

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


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

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({\begin{array}{ccc}m&n&l\\s&\overline {s}&0\end{array}}\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}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 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({\begin{array}{ccc}m&n&l\\s&\overline {s}&0\end{array}}\right)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 }}_{{\mu \nu }}^{{mnl}}(k)

Inverse Fourier-Bessel transform[edit]

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

'Step-up' operations are given by Eq. 53 of Ref. 3. The inverse Hankel transform is

\gamma (r;l_{1}l_{2}ln_{1}n_{2})={\frac  {1}{2\pi ^{2}i^{l}}}\int _{0}^{\infty }{\tilde  {\gamma }}(k;l_{1}l_{2}ln_{1}n_{2})J_{l}(kr)~k^{2}{{\rm {d}}}k

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

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