Improved documentation

doc/intro.texi

@@ -1,221 +1,149 @@
 @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:
 
+The Electron Pair Localization Function (EPLF) is a three-dimensional function 
+measuring locally the electron pairing in a molecular
+system.@footnote{A. Scemama, P. Chaquin and M. Caffarel, @i{J. Chem. Phys.}
+@b{121}, 1725-1735 (2004).}
+An electron @math{i} located at @math{\vec{r}_i} is said to be paired to an electron @math{j}
+located at @math{\vec{r}_j} if electron @math{j} is the closest electron to @math{i}. The
+amount of electron pairing at point @math{{\vec r}} is therefore proportional to the
+inverse of
 @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})}
+ d({\vec r}) = \sum_{i=1,N} \langle \Psi | \delta(\vec{r}-\vec{r}_i) \min_{j\ne i} r_{ij} | \Psi \rangle
 $$
 @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
+where
+@math{d({\vec r})} is the shortest electron-electron distance at @math{{\vec r}},
+@math{\Psi({\vec r}_1,\dots,{\vec r}_N)} is the @math{N}-electron wave function
+and
+@math{r_{ij} = |{\vec r}_j - {\vec r}_i|}.
+
+There are two different types of electron pairs: pairs of electrons with
+the same spin @math{\sigma}, and pairs of electrons with opposite spins @math{\sigma}
+and @math{\bar \sigma}.
+Hence, two quantities are introduced:
 @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 \sigma}({\vec r})  =  \sum_{i=1,N} \langle \Psi | \delta(\vec{r}-\vec{r}_i) \min_{j\ne i ; \sigma_i = \sigma_j} r_{ij} | \Psi \rangle 
 $$
-
 $$
-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 
+ d_{\sigma {\bar \sigma}}({\vec r})  =  \sum_{i=1,N} \langle \Psi | \delta(\vec{r}-\vec{r}_i) \min_{j ; \sigma_i \ne \sigma_j} r_{ij} | \Psi \rangle 
 $$
 @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 electron pair localization function is bounded in the @math{[-1,1]} interval,
+and is defined as
+@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
+When the pairing of spin-unlike electrons is predominant, @math{d_{\sigma \sigma}
+> d_{\sigma {\bar \sigma}}} and @math{{\rm EPLF}({\vec r}) > 0}.
+When the pairing of spin-like electrons is predominant, @math{d_{\sigma \sigma}
+< d_{\sigma {\bar \sigma}}} and @math{{\rm EPLF}({\vec r}) < 0}.
+When the electron pairing of spin-like and spin-unlike electrons is equivalent,
+@math{{\rm EPLF}({\vec r}) \sim 0}.
 
-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}} .
+This localization function does not depend on the type of wave function, as
+opposed to the Electron Localization Function@footnote{A. D. Becke and K. E. Edgecombe
+@i{J. Chem. Phys.} @b{92}, 5397 (1990).} (ELF) where the
+wave function has to be expressed on a single-electron basis set.
+The EPLF can therefore measure electron pairing using any kind of wave
+function: Hartree-Fock (HF), configuration interaction (CI),
+multi-configurational self consistent field (MCSCF) as well as Slater-Jastrow,
+Diffusion Monte Carlo (DMC), Hylleraas wave functions, @i{etc}.
+Due to the presence of the @math{\min} function in the definition of @math{d_{\sigma
+\sigma}} and @math{d_{\sigma {\bar \sigma}}} these quantities cannot be evaluated
+analytically, so quantum Monte Carlo (QMC) methods have been
+used to compute three-dimensional EPLF grids via a statistical sampling of
+@math{|\Psi({\vec r}_1,\dots,{\vec r}_N)|^2}.
 
 @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
+Recently, we proposed an alternative definition of the EPLF which is
+analytically computable.@footnote{F. Alary, J.-L. Heully, A. Scemama, B. Garreau-de-Bonneval, K. I. Chane-Ching and M. Caffarel
+@i{Theor. Chem. Acc.} @b{126}, 243 (2010).}
+The @math{\min} function is approximated
+using:
 @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 } ) }
+\min_{j\neq i} r_{ij} = \lim_{\gamma \rightarrow \infty}
+\sqrt{ -{1\over\gamma} \ln f(\gamma;r_{ij})} 
 $$
-@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}}
+ f(\gamma;r_{ij}) =
+\sum_{j \neq i} e^{ -\gamma r_{ij}^2 }
 $$
 @end tex
 
-@subsection Single determinant wave functions
-
+As @math{f(\gamma ; r_{ij})} has small
+fluctuations in the regions of interest (bonding regions, lone pairs,@dots),
+the integrals
 @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'}
+\left\langle \Psi \left| \delta({\vec r}-{\vec r}_i) \left( \lim_{\gamma \rightarrow \infty} \sqrt{ -{1 \over \gamma} \ln f(\gamma;r_{ij})} \right) \right| \Psi \right\rangle
 $$
 @end tex
-
-@section Multi-determinant wave functions
-
+can be approximated as
 @tex
 $$
-\Psi = \sum_k c_k D_k
+\lim_{\gamma \rightarrow \infty} \sqrt{ -{1 \over \gamma} \ln 
+\langle \Psi | \delta({\vec r}-{\vec r}_i) f(\gamma;r_{ij}) | \Psi \rangle }
 $$
-$$
-\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:
+The expectation values of the minimum distances are now given by:
 @tex
 $$
-\gamma(\vec{r}) = (-\ln \epsilon) \left( {1\over N}{4\over 3}\pi \rho(\vec{r}) \right)^{2/3}
+d_{\sigma \sigma}(\vec{r})  \sim 
+\lim_{\gamma \rightarrow \infty} \sqrt{-{1 \over \gamma} \ln
+f_{\sigma \sigma}(\gamma;\vec{r})} 
+$$
+$$
+d_{\sigma \bar{\sigma}} (\vec{r})  \sim 
+\lim_{\gamma \rightarrow \infty} \sqrt{-{1 \over \gamma} \ln
+f_{\sigma \bar{\sigma}} (\gamma;\vec{r})}
 $$
 @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
+with the two-electron integrals:
 @tex
 $$
-d_{\rm max} = \sqrt{-{1 \over \gamma} \ln \epsilon}
+f_{\sigma \sigma}(\gamma;\vec{r}) = 
+\left \langle \Psi \left| \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 } \right| \Psi \right\rangle  
+$$
+$$
+f_{\sigma \bar{\sigma}} (\gamma;\vec{r})  = 
+\left \langle \Psi \left| \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 } \right| \Psi \right \rangle
 $$
 @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
+If @math{\gamma} is not large enough in regions where
+electrons are close to each other, the value of @math{f_{\sigma \sigma}} or
+@math{f_{\sigma \bar{\sigma}}} can be greater than 1. In that case, the
+approximation of the min function is irrelevant since its value becomes negative.
+On the other hand, if @math{\gamma} is too large @math{f_{\sigma \sigma}} and
+@math{f_{\sigma \bar{\sigma}}} have almost always extremely small values. In this
+case, the values of @math{f_{\sigma \sigma}} and @math{f_{\sigma \bar{\sigma}}} are in practice not
+accurately representable with 64 bit floating point numbers, as they may be
+rounded to zero.  We propose @math{\gamma} to depend on the density, such that the
+largest representable value of @math{d_{\sigma \sigma}} is equal to the radius of a
+sphere that contains almost 0.1 electron:
 @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}
+\gamma(\vec{r}) = \left({4\pi \over 3 n} \rho(\vec{r}) \right)^{2/3} (-\ln (\epsilon)) 
 $$
 @end tex
+where @math{n=0.1} and @math{\epsilon} is the smallest floating point number representable with
+64 bits (@math{\sim 2.10^{-308}}).
+
 
-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.
 

scripts/to_ezfio.py

@@ -95,15 +95,13 @@ def write_ezfioFile(res,filename):
   ezfio.ao_basis_ao_expo = expo
 
 # MOs
-  NumOrbSym = [ s[1] for s in res.symmetries ]
-  mo_tot_num = sum(NumOrbSym)
-  ezfio.mo_basis_mo_tot_num = mo_tot_num
-
   MoTag = res.mo_types[-1]
   if MoTag == 'Natural':
     res.mo_types.pop()
     MoTag = res.mo_types[-1]
-  print MoTag
+  mo_tot_num = len(res.mo_sets[MoTag])
+  ezfio.mo_basis_mo_tot_num = mo_tot_num
+
   try:
     newocc = [0. for i in range(mo_tot_num)]
     while len(res.occ_num[MoTag]) < mo_tot_num: