Editing Computational implementation of integral equations

Jump to navigation Jump to search
Warning: You are not logged in. Your IP address will be publicly visible if you make any edits. If you log in or create an account, your edits will be attributed to your username, along with other benefits.

The edit can be undone. Please check the comparison below to verify that this is what you want to do, and then publish the changes below to finish undoing the edit.

Latest revision Your text
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]]
Please note that all contributions to SklogWiki are considered to be released under the Creative Commons Attribution Non-Commercial Share Alike (see SklogWiki:Copyrights for details). If you do not want your writing to be edited mercilessly and redistributed at will, then do not submit it here.
You are also promising us that you wrote this yourself, or copied it from a public domain or similar free resource. Do not submit copyrighted work without permission!

To edit this page, please answer the question that appears below (more info):

Cancel Editing help (opens in new window)