Consider two lines of finite length $L$ containing a homogeneous density
of LJ centers. Let $\vec{e}_{1,2}$ be the direction vectors,
$\vec{r}_{12}$ the vector between the centers of the two line
segments, and $\lambda, \, \mu$ parameters giving the positions of
the interacting points along $1$ and $2$. The squared distance between
any two such points is given by
$r^{2}= \left|
r_{12}+\mu \cdot \vec{e}_{2}-\lambda \cdot \vec{e}_{1}
\right|^{2}
$.
The total interaction energy between the two lines is then
We define the triple Gaussian approximation function as
$
g_{3G}(r)=\sum_{i=1}^{3}A_{i} \exp\{-B_{i} r^{2}\}
$
and demand that it equals $u_{LJ}(r)$ at the points
$r=0.7; 0.8; 0.9; 1.0; r_{min}; 1.6$,
where $r_{min}=1.12246$ is the minimum point of the LJ
potential. The six parameters are then given by
Figure 1:
The LJ potential approximated by a sum of three Gaussians:
(a) logarithm of $u_{LJ}(r)$ and $g_{3G}(r)$ in the range $r=0.7 - 0.99$;
(b) $u(r)$ and $g_{3G}(r)$, range $r= 0.99 - 3$.
Solid line: $u_{LJ}$; dashed: $g_{3G}$; the dashed-dot lines denote the deviation
$g_{3G}-u_{LJ}$.
In place of $\lambda, \mu$ we introduce new parameters $\gamma, \delta $ along the
two carrier lines, such that $\gamma = \delta = 0$ at the points of least line-line
distance $r_{0}$. The squared distance between two points along the line segments
is then
with $d = \vec{e}_{1} \cdot \vec{e}_{2}$. Thus each of the three Gauss terms
contributes
$
G_{i}=C_{i} \int \limits_{\gamma_{a}}^{\gamma_{b}} d \gamma
\int \limits_{\delta_{a}}^{\delta_{b}} d \delta \,
\exp \left\{ - B_{i} \left(
\gamma^{2}+\delta^{2}-2 \gamma \delta d \right)
\right\}
$
with $C_{i} \equiv A_{i}\exp\{-B_{i} r_{0}^{2} \}$.
This is the integral of a bivariate Gaussian over a rectangular
region. It can be
decomposed into a linear combination of four quadrant integrals:
$G_{i}= H_{i}(b,b,d)-H_{i}(a,b,d)-H_{i}(b,a,d)+H_{i}(a,a,d) $, where
where $\mathrm{He}_{k}$ are the Hermite polynomials
defined by the recurrence relation
$\mathrm{He}_{n+1}(x)=x \mathrm{He}_{n}(x)-n \mathrm{He}_{n-1}(x)$
with $\mathrm{He}_{0}(x)=1$ and $\mathrm{He}_{1}(x)=x$. (This variant
is sometimes called the "probabilists' Hermite polynomials".)
$n(x)$ denotes the 1D normal distribution with cumulative probability
$N(x)$. (Note that $N(x)=(1/2) [1+erf(x/\sqrt{2}]$.)
Vasicek's series, convergent for $0.5 \leq \rho^{2} \leq 1.0$:
Since for $xy \neq 0$
$
N_{2}(x,y,\rho) =
N_{2} \left(
x,0,\mathrm{sign}(x) ( \rho x - y)/\sqrt{x^{2}-2\rho x y +y^{2}}
\right)
+N_{2}\left(
y,0,\mathrm{sign} (y)( \rho y - x)/\sqrt{x^{2}-2\rho x y +y^{2}}
\right)
-s_{0}
$
where $s_{0}=0$ or $0.5$ for $xy > 0$ or $xy < 0$, respectively, it
suffices to evaluate $N_{2}(x,0,\rho)$.
(where $N(x)=\left(1 / \sqrt{2 \pi}\right)
\int_{- \infty}^{x} du \exp\{ -u^{2}/2\} $,
the error integral in 1 D) and accumulate the sum
$
Q = \sum_{k=0}^{m} \, A_{k}
$
where $C_{i} \equiv A_{i}\exp\{-B_{i} r_{0}^{2} \}$, and
$r_{0}$ is the shortest distance between the carrier lines.
The terms $G_{i}$ are to be evaluated using either the
tetrachoric or Vasicek's series.
Number of terms? --
in most circumstances the series may be truncated at
$m=1$ to $8$, with an absolute energy error below $10^{-6}$.
Computing times? --
A comparison of computing times is given in Table 1.
$m=5$ was used in both configurations, tee and end-to-end.
LJ lines
6 LJ sites
Gay-Berne
end-to-end
0.858e-04
0.216e-04
0.172e-05
tee
0.412e-04
0.216e-04
0.119e-05
Computing times (seconds) per pair interaction, for
various models. (Pentium M CPU; 1.7 GHz; F77 compiler.)
The "tee" configuration corresponds to zero correlation in
the bivariate Gaussian representation, while an almost parallel
(e.g. end-to-end) configuration relates to $\rho \rightarrow 1$
and thus requires the use of a Vasicek series; see text.
For comparison, the potential between two 6-LJ-sites particles
as well as the Gay-Berne potential are computed.
Compare with other models:
In the figure the LJL interaction is compared
to multisite potentials with $9 \dots 72$ LJ sites;
"tee" configuration; axis length $L=4$; total
"LJ load" $Q=1$.
Figure 2:
The LJL model as the limiting case of multi-LJ-sites potentials.
"tee" position; length of axis: $L=4$.
Solid line: LJL; dashed: 9LJ; remaining graphs, approaching LJL:
$n=18, 36$ and $72$.
Next we compare the LJL model with the 4LJ potential in
the three basic pair configurations.
Figure 3:
Comparison of the LJ lines (LJL) model
with the four-LJ-sites (4LJ) potential with axis length
$L=2$. Curves at left: side-by-side configuration
with $d \equiv \vec{e}_{1} \cdot \vec{e}_{2}=0.99$;
middle: tee configuration; right: end-to-end with $d$ as in
side-by-side.
Solid lines: LJL potential, dashed: 4LJ potential.
Pair Force:
The pair force
$F_{12}$ acting from 2 on 1 is given by:
where
$ \nabla_{\textstyle e_{1}} \equiv
\left\{ \partial / \partial e_{1}^{\alpha} \right\}
$
is the derivative with respect to the direction vector
$\vec{e}_{1}$. To the best of my knowledge at this date (Sep-07)
this derivative has the following form:
Remarks regarding perturbation theory:
Hard-body correlate: hard spherocylinder model?
Danger: do NOT divide the elementary interaction in the
definition of the LJL potential according to WCA,
using just the repulsive part of the LJ potential to
define the reference interaction. The attractive part
is important as it compensates part of the repulsion.
Example:
$L=4$; end-to-end configuration with distance $r=0.7$
has a considerable probability in the full LJL system
but not in the purely repulsive case.
Thus the assumption of identical structure in the reference
and target systems is not met.
The hard body correlate must be found in another way.
Any ideas?