Computational implementation of integral equations: Difference between revisions
Carl McBride (talk | contribs) No edit summary |
Carl McBride (talk | contribs) No edit summary |
||
Line 15: | Line 15: | ||
is to perform a [[Fast Fourier transform |(fast) Fourier transform]] (FFT). | is to perform a [[Fast Fourier transform |(fast) Fourier transform]] (FFT). | ||
Note: convergence is poor for liquid densities. (See Ref.s 1 to 6). | Note: convergence is poor for liquid densities. (See Ref.s 1 to 6). | ||
==Picard iteration== | ==Picard iteration== | ||
===Closure relation=== | |||
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) {\bf Closure relation}: $\gamma_{mns}^{\mu \nu} (r) \rightarrow c_{mns}^{\mu \nu} (r)$\\ | |||
(Note: for linear fluids $\mu = \nu =0$)\\ | |||
~\\ | |||
i) Perform the summation | |||
\begin{equation} | |||
g(12)=g({\bf r}_{12},\omega_1,\omega_2)=\sum_{mns\mu \nu} g_{mns}^{\mu \nu}({\bf r}_{12}) \Psi_{\mu \nu s}^{mn}(\omega_1,\omega_2) | |||
\end{equation} | |||
where ${\bf 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 | |||
\begin{equation} | |||
\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) | |||
\end{equation} | |||
with $\overline{s} = -s$.\\ | |||
~\\ | |||
ii) Define the variables | |||
\begin{equation} | |||
x_1= \cos \theta_1 | |||
\end{equation} | |||
\begin{equation} | |||
x_2= \cos \theta_2 | |||
\end{equation} | |||
\begin{equation} | |||
z_1 = \cos \chi_1 | |||
\end{equation} | |||
\begin{equation} | |||
z_2 = \cos \chi_2 | |||
\end{equation} | |||
\begin{equation} | |||
y= \cos \phi_{12} | |||
\end{equation} | |||
Thus $\gamma(12)=\gamma (r,x_1x_2,y,z_1z_2)$.\\ | |||
~\\ | |||
iii) 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 | |||
\begin{equation} | |||
\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}) | |||
\end{equation} | |||
where | |||
\begin{equation} | |||
\hat{d}_{s \mu}^m (x) = (2m+1)^{1/2} d_{s \mu}^m(\theta) | |||
\end{equation} | |||
where $d_{s \mu}^m(\theta)$ is the angular, $\theta$, part of the | |||
rotation matrix $\mathcal{D}_{s \mu}^m (\omega)$,\\ | |||
and | |||
\begin{equation} | |||
e_s(y)=\exp(is\phi) | |||
\end{equation} | |||
\begin{equation} | |||
e_{\mu}(z)= \exp(i\mu \chi) | |||
\end{equation} | |||
For the limits in the summations | |||
\begin{equation} | |||
\begin{equation} | |||
L_1= \max (s,\nu_1) | |||
\end{equation} | |||
\begin{equation} | |||
L_2= \max (s,\nu_2) | |||
\end{equation} | |||
The above equation constitutes a separable five-dimensional transform. To rapidly evaluate | |||
this expression it is broken down into five one-dimensional transforms: | |||
\begin{equation} | |||
\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}) | |||
\end{equation} | |||
\begin{equation} | |||
\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}) | |||
\end{equation} | |||
\begin{equation} | |||
\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) | |||
\end{equation} | |||
\begin{equation} | |||
\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}) | |||
\end{equation} | |||
\begin{equation} | |||
\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}) | |||
\end{equation} | |||
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.\\ | |||
~\\ | |||
iv) Integrate over angles $c_2(12)$:\\ | |||
~\\ | |||
Use Gauss-Legendre quadrature for $x_1$ and $x_2$\\ | |||
Use Gauss-Chebyshev quadrature for $y$, $z_1$ and $z_2$\\ | |||
thus | |||
\begin{equation} | |||
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}) | |||
\end{equation} | |||
where the Gauss-Legendre quadrature weights are given by | |||
\begin{equation} | |||
w_i= \frac{1}{(1-x_i^2)}[P_{NG}^{'} (x_i)]^2 | |||
\end{equation} | |||
while the Gauss-Chebyshev quadrature has the constant weight | |||
\begin{equation} | |||
w=\frac{1}{NG} | |||
\end{equation} | |||
===Perform FFT from Real to Fourier space=== | |||
==Ng acceleration== | ==Ng acceleration== | ||
Revision as of 11:50, 30 May 2007
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;
- trial solution for
- calculate
- 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
Closure relation
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) {\bf Closure relation}: $\gamma_{mns}^{\mu \nu} (r) \rightarrow c_{mns}^{\mu \nu} (r)$\\ (Note: for linear fluids $\mu = \nu =0$)\\ ~\\ i) Perform the summation \begin{equation} g(12)=g({\bf r}_{12},\omega_1,\omega_2)=\sum_{mns\mu \nu} g_{mns}^{\mu \nu}({\bf r}_{12}) \Psi_{\mu \nu s}^{mn}(\omega_1,\omega_2) \end{equation} where ${\bf 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 \begin{equation} \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) \end{equation} with $\overline{s} = -s$.\\ ~\\ ii) Define the variables \begin{equation} x_1= \cos \theta_1 \end{equation} \begin{equation} x_2= \cos \theta_2 \end{equation} \begin{equation} z_1 = \cos \chi_1 \end{equation} \begin{equation} z_2 = \cos \chi_2 \end{equation} \begin{equation} y= \cos \phi_{12} \end{equation} Thus $\gamma(12)=\gamma (r,x_1x_2,y,z_1z_2)$.\\ ~\\ iii) 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 \begin{equation} \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}) \end{equation} where \begin{equation} \hat{d}_{s \mu}^m (x) = (2m+1)^{1/2} d_{s \mu}^m(\theta) \end{equation} where $d_{s \mu}^m(\theta)$ is the angular, $\theta$, part of the rotation matrix $\mathcal{D}_{s \mu}^m (\omega)$,\\ and \begin{equation} e_s(y)=\exp(is\phi) \end{equation} \begin{equation} e_{\mu}(z)= \exp(i\mu \chi) \end{equation} For the limits in the summations \begin{equation} \begin{equation} L_1= \max (s,\nu_1) \end{equation} \begin{equation} L_2= \max (s,\nu_2) \end{equation} The above equation constitutes a separable five-dimensional transform. To rapidly evaluate this expression it is broken down into five one-dimensional transforms: \begin{equation} \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}) \end{equation} \begin{equation} \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}) \end{equation} \begin{equation} \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) \end{equation} \begin{equation} \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}) \end{equation} \begin{equation} \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}) \end{equation} 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.\\ ~\\ iv) Integrate over angles $c_2(12)$:\\ ~\\ Use Gauss-Legendre quadrature for $x_1$ and $x_2$\\ Use Gauss-Chebyshev quadrature for $y$, $z_1$ and $z_2$\\ thus \begin{equation} 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}) \end{equation} where the Gauss-Legendre quadrature weights are given by \begin{equation} w_i= \frac{1}{(1-x_i^2)}[P_{NG}^{'} (x_i)]^2 \end{equation} while the Gauss-Chebyshev quadrature has the constant weight \begin{equation} w=\frac{1}{NG} \end{equation}
Perform FFT from Real to Fourier space
Ng acceleration
References
- M. J. Gillan "A new method of solving the liquid structure integral equations" Molecular Physics 38 pp. 1781-1794 (1979)
- Stanislav Labík, Anatol Malijevský and Petr Voncaronka "A rapidly convergent method of solving the OZ equation", Molecular Physics 56 pp. 709-715 (1985)
- F. Lado "Integral equations for fluids of linear molecules I. General formulation", Molecular Physics 47 pp. 283-298 (1982)
- F. Lado "Integral equations for fluids of linear molecules II. Hard dumbell solutions", Molecular Physics 47 pp. 299-311 (1982)
- F. Lado "Integral equations for fluids of linear molecules III. Orientational ordering", Molecular Physics 47 pp. 313-317 (1982)
- Enrique Lomba "An efficient procedure for solving the reference hypernetted chain equation (RHNC) for simple fluids" Molecular Physics 68 pp. 87-95 (1989)