diff --git a/bin/to_ezfio.exe b/bin/to_ezfio.exe index be798ba..4c4215a 100755 Binary files a/bin/to_ezfio.exe and b/bin/to_ezfio.exe differ diff --git a/doc/eplf.pdf b/doc/eplf.pdf index 10d62cb..11dd2eb 100644 Binary files a/doc/eplf.pdf and b/doc/eplf.pdf differ diff --git a/doc/intro.texi b/doc/intro.texi index 9971621..5235170 100644 --- a/doc/intro.texi +++ b/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. diff --git a/scripts/to_ezfio.py b/scripts/to_ezfio.py index 2bf9d9c..22876a7 100755 --- a/scripts/to_ezfio.py +++ b/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: