Computational implementation of integral equations: Difference between revisions

From SklogWiki
Jump to navigation Jump to search
No edit summary
mNo edit summary
 
(24 intermediate revisions by 2 users not shown)
Line 1: Line 1:
Integral equations are solved numerically.
One has the [[Ornstein-Zernike relation]], <math>\gamma (12)</math>
and a [[closure relations | 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:
===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>\left. \gamma(12) \right. =\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 polynomials |Legendre polynomial]] <math>P_\nu(cos \theta)</math>
where <math>y_j</math> are the  <math>\nu</math> roots of the [[Chebyshev polynomials |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>\left. e_s(y) \right.=\exp(is\phi)</math>
:<math>\left. e_{\mu}(z) \right.= \exp(i\mu \chi)</math>
For the limits in the summations
:<math>\left. L_1 \right.= \max (s,\nu_1)</math>
:<math>\left. L_2 \right.= \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====
:<math>c_{mns}^{\mu \nu} (r)  \rightarrow  c_{\mu \nu}^{mnl} (r)</math>
this is done using the Blum transformation (Refs 7, 8 and 9):
:<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 in Ref. 7 or Lado Eq. 39 in Ref. 3),
where <math>J_l(x)</math> is a [[Bessel functions |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  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.
:<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====
:<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>
===Ornstein-Zernike relation <math>\tilde{c}_{mns}^{\mu \nu} (k)  \rightarrow  \tilde{\gamma}_{mns}^{\mu \nu} (k)</math>===
For simple fluids:
:<math>\tilde{\gamma}(k)= \frac{\rho \tilde{c}_2 (k)^2}{1- \rho  \tilde{c}_2 (k)}</math>
For molecular fluids (see Eq. 19 of Lado Ref. 3)
:<math>\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)</math>
where <math>\tilde{{\mathbf S}}_{m}(k)</math> and <math>\tilde{\mathbf C}_{m}(k)</math> are matrices
with elements <math>\tilde S_{l_1 l_2 m}(k), \tilde{C}_{l_1 l_2 m}(k), l_1,l_2 \geq m</math>.
For mixtures of simple fluids  (see Ref. 10 Juan Antonio Anta PhD thesis pp. 107--109):
:<math>\tilde{\Gamma}(k) =  {\mathbf D}  \left[{\mathbf I} -  {\mathbf D}  \tilde{\mathbf C}(k)\right]^{-1} \tilde{\mathbf C}(k)\tilde{\mathbf C}(k)</math>
===Conversion back from Fourier space to Real space===
:<math>\tilde{\gamma}_{mns}^{\mu \nu} (k)  \rightarrow \gamma_{mns}^{\mu \nu} (r)</math>
(basically the inverse of step 2).
====Axial reference frame to spatial reference frame====
:<math>\tilde{\gamma}_{mns}^{\mu \nu} (k) \rightarrow  \tilde{\gamma}^{mnl}_{\mu \nu} (k)</math>
====Inverse Fourier-Bessel transform====
:<math>\tilde{\gamma}^{mnl}_{\mu \nu} (k) \rightarrow  \gamma^{mnl}_{\mu \nu} (r)</math>
'Step-up' operations are given by Eq. 53 of Ref. 3. The inverse Hankel transform is
:<math>\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</math>
====Change from  spatial reference frame back to  axial reference frame====
:<math>\gamma^{mnl}_{\mu \nu} (r) \rightarrow  \gamma_{mns}^{\mu \nu} (r)</math>.
==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)]
==Angular momentum coupling coefficients==
*[http://dx.doi.org/10.1016/0010-4655(70)90034-2  Taro Tamura  "Angular momentum coupling coefficients", Computer Physics Communications  '''1''' pp.  337-342 (1970)]
*[http://dx.doi.org/10.1016/0010-4655(71)90030-0 J. G. Wills  "On the evaluation of angular momentum coupling coefficients", omputer Physics Communications  '''2''' pp. 381-382 (1971)]
==References==
==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/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/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)]
#[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)]
#[http://dx.doi.org/10.1063/1.1676864 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)]
#[http://dx.doi.org/10.1063/1.1678503 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)]
#[http://dx.doi.org/10.1063/1.1679655 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)]
#[http://dx.doi.org/10.1063/1.454286    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)]
[[category: integral equations]]

Latest revision as of 15:57, 31 January 2008

Integral equations are solved numerically. One has the Ornstein-Zernike relation, and a closure relation, (which incorporates the bridge function ). The numerical solution is iterative;

  1. trial solution for
  2. calculate
  3. use the Ornstein-Zernike relation to generate a new etc.

Note that the value of is local, i.e. the value of at a given point is given by the value of 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 [edit]

(Note: for linear fluids )

Perform the summation[edit]

where is the separation between molecular centers and the sets of Euler angles needed to specify the orientations of the two molecules, with

with .

Define the variables[edit]

Thus

.

Evaluate[edit]

Evaluations of are performed at the discrete points where the are the roots of the Legendre polynomial where are the roots of the Chebyshev polynomial and where are the roots of the Chebyshev polynomial thus

where



where is the angular, , part of the rotation matrix , and

For the limits in the summations

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

Operations involving the and 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". and are parameters; is the number of nodes in the Gauss integration, and the the max index in the truncated rotational invariants expansion.

Integrate over angles [edit]

Use Gauss-Legendre quadrature for and Use Gauss-Chebyshev quadrature for , and . Thus

where the Gauss-Legendre quadrature weights are given by

while the Gauss-Chebyshev quadrature has the constant weight

Perform FFT from Real to Fourier space [edit]

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

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

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

Fourier-Bessel Transforms[edit]

(see Blum and Torruella Eq. 5.6 in Ref. 7 or Lado Eq. 39 in Ref. 3), where is a Bessel function of order . `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.


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

this is done using the Blum transformation

Ornstein-Zernike relation [edit]

For simple fluids:

For molecular fluids (see Eq. 19 of Lado Ref. 3)

where and are matrices with elements .

For mixtures of simple fluids (see Ref. 10 Juan Antonio Anta PhD thesis pp. 107--109):

Conversion back from Fourier space to Real space[edit]

(basically the inverse of step 2).

Axial reference frame to spatial reference frame[edit]

Inverse Fourier-Bessel transform[edit]

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

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

.

Ng acceleration[edit]

Angular momentum coupling coefficients[edit]

References[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)