@section QMC Formulation
The EPLF@footnote{``Electron pair localization function, a practical tool to visualize electron
localization in molecules from quantum Monte Carlo data'', A. Scemama, P. Chaquin and M. Caffarel,
@i{J. Chem. Phys.} @b{121}, 1725-1735 (2004).
}
has been designed to describe local electron pairing in molecular
systems. It is defined as a scalar function defined in the three-dimensional
space and taking its values in the [-1,1] range. It is defined as follows:

@tex
$$
{\rm EPLF}(\vec{r}) = { d_{\sigma \sigma} (\vec{r}) - d_{\sigma {\bar \sigma}} (\vec{r}) \over d_{\sigma \sigma} (\vec{r}) + d_{\sigma {\bar \sigma}}
(\vec{r})}
$$
@end tex
where the quantity @math{ d_{\sigma \sigma} (\vec{r}) } [@math{resp. d_{\sigma
\bar \sigma} (\vec{r})} ] denotes the quantum-mechanical average of the
distance between an electron of spin @math{\sigma} located at @math{\vec r} and the closest
electron of same spin (resp., of opposite spin @math{\bar\sigma}).

The mathematical definition of these quantities can be written as
@tex
$$
d_{\sigma \sigma}(\vec{r}) = \int \Psi^2(\vec{r}_1,\dots,\vec{r_N}) \left[
\sum_{i=1}^N \delta(\vec{r}-\vec{r}_i) \min_{j\neq
i;\sigma_j=\sigma_i}|\vec{r}_i - \vec{r}_j| \right] d\vec{r}_1 \dots d\vec{r}_N
$$

$$
d_{\sigma {\bar \sigma}}(\vec{r}) = \int \Psi^2(\vec{r}_1,\dots,\vec{r}_N)
\left[ \sum_{i=1}^N \delta(\vec{r}-\vec{r}_i) \min_{j\ne
i;\sigma_j\neq\sigma_i}|\vec{r}_i - \vec{r}_j| \right] d\vec{r}_1 \dots
d\vec{r}_N 
$$
@end tex
where @math{\sigma} is the spin (\alpha or \beta), @math{\bar \sigma} is the
spin opposite to @math{\sigma} @math{\Psi(\vec{r}_1,\dots,\vec{r}_N)} is the
wave function, and @math{N} is the number of electrons.

In a region of space, if the shortest distance separating anti-parallel
electrons is smaller than the shortest distance separating electrons of same
spin, the EPLF takes positive values and indicates pairing of anti-parallel
electrons. In contrast, if the shortest distance separating anti-parallel
electrons is larger than the shortest distance separating electrons of same
spin, the EPLF takes negative values and indicates pairing of parallel
electrons (which in practice never happens). If the shortest distance
separating anti-parallel electrons is equivalent to the shortest distance
separating electrons of same spin, the EPLF takes values close to zero and
indicates no electron pairing.

The original formulation of EPLF is extremely easy to compute in the quantum
Monte Carlo framework. However, it is not possible to compute it analytically
due to the presence of the min function in the definitions of @math{d_{\sigma
\sigma}} and @math{d_{\sigma \bar \sigma}} .

@section Analytical Formulation

@subsection The EPLF operator

To evaluate exactly the average of the min function is not possible. We
proposed@footnote{``Structural and optical properties of a neutral Nickel bisdithiolene
complex. Density Functional versus ab-initio methods.'', F. Alary, J.-L.
Heully, A. Scemama, B. Garreau-de-Bonneval, K. I. Chane-Ching and M. Caffarel,
@i{Theoretical Chemistry Accounts}, DOI: 10.1007/s00214-009-0679-9}
to introduce a representation of this function in terms of gaussians
and to construct our @math{d}-functions in a slightly different way so that the
gaussian contributions can be exactly integrated out. The representation of the
min function we are considering here is written as
@tex
$$
\min_{j\neq i} |\vec{r}_i - \vec{r}_j| = \lim_{\gamma \rightarrow \infty}
\sqrt{ -{1 \over \gamma} \ln ( \sum_{j \neq i} e^{ -\gamma |\vec{r}_i -
\vec{r}_j|^2 } ) }
$$
@end tex
and we propose to define the new modified average-distances as
@tex
$$
d_{\sigma \sigma}(\vec{r}) \sim \lim_{\gamma \rightarrow \infty}
\sqrt{-{1 \over \gamma} \ln \int \Psi^2(\vec{r}_1,\dots,\vec{r}_N)
\sum_{i=1}^{N} \delta(\vec{r}-\vec{r}_i) \sum_{j \neq i ; \sigma_i =
\sigma_j}^{N} e^{ -\gamma |\vec{r}_i - \vec{r}_j|^2 } d\vec{r}_1 \dots
d\vec{r}_N}
$$

$$
d_{\sigma \bar{\sigma}} (\vec{r}) \sim \lim_{\gamma \rightarrow \infty}
\sqrt{-{1\over \gamma} \ln \int \Psi^2(\vec{r}_1,\dots,\vec{r}_N)
\sum_{i=1}^{N} \delta(\vec{r}-\vec{r}_i) \sum_{j\neq i ; \sigma_i \neq
\sigma_j}^{N} e^{ -\gamma |\vec{r}_i - \vec{r}_j|^2 } d\vec{r}_1 \dots
d\vec{r}_N} $$
@end tex

Note that these new definitions of the average-distances are different from the
previous ones essentially because @math{\ln \langle X \rangle \neq \langle \ln
X\rangle} where @math{X} denotes the sum @math{\sum_{j} e^{ -\gamma |\vec{r} -
\vec{r}_j|^2 }} and the symbol @math{\langle\dots\rangle} denotes the quantum
average. Remark that the equality is almost reached when the
fluctuations of @math{X} are small. In our applications it seems that we are
not too far from this regime and we have thus systematically observed that both
definitions of the @math{d}-functions lead to similar EPLF landscapes.

Now, introducing
@tex
$$
f_{\sigma\sigma}^\gamma(\vec{r}) = \sum_{i=1}^{N} \delta(\vec{r}-\vec{r}_i)
\sum_{j\ne i ; \sigma_i = \sigma_j}^{N} e^{ -\gamma |\vec{r}_i - \vec{r}_j|^2 }
$$

$$
f_{\sigma\bar{\sigma}}^\gamma(\vec{r}) = \sum_{i=1}^{N}
\delta(\vec{r}-\vec{r}_i) \sum_{j ; \sigma_i \neq \sigma_j}^{N} e^{ -\gamma
|\vec{r}_i - \vec{r}_j|^2 } $$
@end tex
we obtain
@tex
$$
{\rm EPLF}(\vec{r}) = {\sqrt{-{1\over \gamma} \ln \langle \Psi |
f_{\sigma\sigma}^\gamma(\vec{r}) | \Psi \rangle} - \sqrt{-{1\over \gamma} \ln
\langle \Psi | f_{\sigma\bar{\sigma}}^\gamma(\vec{r}) | \Psi \rangle} \over
\sqrt{-{1\over \gamma} \ln \langle \Psi | f_{\sigma\sigma}^\gamma(\vec{r}) |
\Psi \rangle} + \sqrt{-{1\over \gamma} \ln \langle \Psi |
f_{\sigma\bar{\sigma}}^\gamma(\vec{r}) | \Psi \rangle}}
$$
@end tex

@subsection Single determinant wave functions

@tex
$$
\Psi(r_1,\dots,r_N) = D_0 = |\phi_1,\dots,\phi_N|
$$

$$
\langle \Psi | f_{\sigma\sigma}^\gamma(\vec{r}) | \Psi \rangle = \sum_{\sigma =
\alpha,\beta} \sum_{i=1}^{N_\sigma} \sum_{j\ne i}^{N_\sigma} \int
\phi_i(\vec{r}) \phi_j(\vec{r'}) e^{-\gamma |\vec{r'}-\vec{r}|^2}
\phi_i(\vec{r}) \phi_j(\vec{r'}) d\vec{r'} 
$$
$$
- \int \phi_i(\vec{r})
\phi_j(\vec{r'}) e^{-\gamma |\vec{r'}-\vec{r}|^2} \phi_j(\vec{r})
\phi_i(\vec{r'}) d\vec{r'}
$$

$$
\langle \Psi | f_{\sigma\bar{\sigma}}^\gamma(\vec{r}) | \Psi \rangle =
\sum_{\sigma=\alpha,\beta} \sum_{i=1}^{N_\sigma} \sum_{j=1}^{N_{\bar \sigma}}
\int \phi_i(\vec{r}) \phi_j(\vec{r'}) e^{-\gamma |\vec{r'}-\vec{r}|^2}
\phi_i(\vec{r}) \phi_j(\vec{r'}) d\vec{r'}
$$
@end tex

@section Multi-determinant wave functions

@tex
$$
\Psi = \sum_k c_k D_k
$$
$$
\langle \Psi | f(\vec{r}) | \Psi \rangle = \sum_k \sum_l c_k c_l \langle D_k |
f(\vec{r}) | D_l \rangle $$
$$
\langle D_k | f_{\sigma\sigma}^\gamma(\vec{r}) | D_l \rangle = \sum_{\sigma =
\alpha,\beta} \sum_{i=1}^{N_\sigma} \sum_{j\ne i}^{N_\sigma} \int
\phi^k_i(\vec{r}) \phi^k_j(\vec{r'}) e^{-\gamma |\vec{r'}-\vec{r}|^2}
\phi^l_i(\vec{r}) \phi^l_j(\vec{r'}) d\vec{r'} 
$$
$$
- \int \phi^k_i(\vec{r})
\phi^k_j(\vec{r'}) e^{-\gamma |\vec{r'}-\vec{r}|^2} \phi^l_j(\vec{r})
\phi^l_i(\vec{r'}) d\vec{r'} $$
$$
\langle D_k | f_{\sigma\bar \sigma}^\gamma(\vec{r}) | D_l \rangle =
\sum_{\sigma=\alpha,\beta} \sum_{i=1}^{N_\sigma} \sum_{j=1}^{N_{\bar \sigma}}
\int \phi^k_i(\vec{r}) \phi^k_j(\vec{r'}) e^{-\gamma |\vec{r'}-\vec{r}|^2}
\phi^l_i(\vec{r}) \phi^l_j(\vec{r'}) d\vec{r'} $$
@end tex

@subsection Expression of the gamma parameter

In our definition of the modified EPLF the parameter @math{\gamma} is supposed to be (very)
large. In practice, if it is chosen too small, the integrals @math{\langle \Psi |
f_{\sigma{\bar \sigma}}^\gamma(\vec{r}) | \Psi \rangle} can take values larger
than 1 in regions of high density, giving a negative value for the shortest
average distance. If it is chosen too large numerical instabilities appear
since the values of the integrals become extremely small. We propose to use a
value of @math{\gamma} which depends on the electron density @math{\rho(\vec{r})} as:
@tex
$$
\gamma(\vec{r}) = (-\ln \epsilon) \left( {1\over N}{4\over 3}\pi \rho(\vec{r}) \right)^{2/3}
$$
@end tex
where @math{\epsilon} is the smallest possible floating point number that can
be represented with 64 bits (@math{\sim 2.10^{-308}}). The reason for this
choice is the following.
The largest possible value of the approximate shortest distance can be computed as
@tex
$$
d_{\rm max} = \sqrt{-{1 \over \gamma} \ln \epsilon}
$$
@end tex
If the density is conidered constant in a sphere of radius @math{d_{\rm max}}, the number
@math{N} of electrons contained in the sphere is
@tex
$$
N = {4\over 3} \pi d_{\rm max}^3 \rho(r)
$$
@end tex
The maximum possible distance can be parametrized by a target constant number @math{N}
of electrons in the sphere
@tex
$$
d_{\rm max}(r) = \left( {1\over N} {4\over 3} \pi \rho(r) \right)^{-1/3}
$$
@end tex

If @math{N} is chosen small enough (0.01 electron), the value of
@math{\langle \Psi | f(r) | \Psi \rangle} is very unlikely to exceed 1, and the
value of @math{\gamma} will also be sufficiently large in the valence regions.