From a0eeaf84d32446173689354c98672d4668978ad1 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 29 Apr 2010 17:15:07 +0200 Subject: [PATCH] Improved gamma parameter --- doc/eplf.texi | 57 +++++++++++ doc/ezfio.texi | 23 +++++ doc/intro.texi | 221 ++++++++++++++++++++++++++++++++++++++++ doc/wf.texi | 64 ++++++++++++ src/eplf_function.irp.f | 64 +++++++----- 5 files changed, 402 insertions(+), 27 deletions(-) create mode 100644 doc/eplf.texi create mode 100644 doc/ezfio.texi create mode 100644 doc/intro.texi create mode 100644 doc/wf.texi diff --git a/doc/eplf.texi b/doc/eplf.texi new file mode 100644 index 0000000..8cd06c1 --- /dev/null +++ b/doc/eplf.texi @@ -0,0 +1,57 @@ +\input texinfo @c -*-texinfo-*- +@c %**start of header (This is for running texinfo on a region.) +@setfilename eplf.info +@settitle Electron Pair Localization Function +@afourpaper +@c %**end of header (This is for running texinfo on a region.) +@macro exe +{@acronym{EPLF} } +@end macro + +@ignore +@ifinfo +@format +START-INFO-DIR-ENTRY +* EPLF: (eplf.info). Electron Pair Localization Function +END-INFO-DIR-ENTRY +@end format +@end ifinfo +@end ignore + +@ifinfo +This file documents @exe, a code to compute the Electron Pair +Localization Function.@refill +@end ifinfo + +@titlepage +@title @exe +@subtitle Electron Pair Localization Function: User's Guide +@subtitle May, 2010 +@author @email{scemama@@irsamc.ups-tlse.fr,Anthony Scemama} +@end titlepage + +@ifinfo +@node Top, (dir), (dir) +@top General Introduction +@c Preface or Licensing nodes should come right +@c after the Top node, in `unnumbered' sections, +@c then the chapter, `What is foogol'. +This file documents @exe, a code to compute the Electron Pair +Localization Function.@refill +@end ifinfo + +@menu +* Introduction:: Definition the Electron Pair Localization Function +* Using EPLF :: A basic introduction to using @exe. +@end menu + +@node Introduction +@chapter Introduction +@include intro.texi + +@node Using EPLF +@chapter Using the @exe program +@include wf.texi +@include ezfio.texi + +@bye diff --git a/doc/ezfio.texi b/doc/ezfio.texi new file mode 100644 index 0000000..10bc9c4 --- /dev/null +++ b/doc/ezfio.texi @@ -0,0 +1,23 @@ +@macro ezfio +{@acronym{EZFIO} } +@end macro + +@section Preparing the @ezfio file + +The @exe program uses the @ezfio (Easy Fortran Input/Output) library +generator@footnote{@url{http://ezfio.sourceforge.net}} to handle data persistence. +An @ezfio ``file'' is a directory which contains all the needed input/output data +needed by the code. These files are accessible via a simple @acronym{API} both +in Fortran and Python. + +The output file containing the computed wave function has to be transformed +into an @ezfio file in order to be used by @exe. This step is realized by +calling the @command{to_ezfio.py} Python script followed by the name of the +file containing the wave function: +@example +to_ezfio.py test.out +@end example +An @ezfio file is produced, named @file{test.out.ezfio}. + + + diff --git a/doc/intro.texi b/doc/intro.texi new file mode 100644 index 0000000..9971621 --- /dev/null +++ b/doc/intro.texi @@ -0,0 +1,221 @@ +@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. + diff --git a/doc/wf.texi b/doc/wf.texi new file mode 100644 index 0000000..b02242d --- /dev/null +++ b/doc/wf.texi @@ -0,0 +1,64 @@ +@section Wave function preparation + +@macro gamess +{@acronym{GAMESS} } +@end macro + +@macro mcscf +{@acronym{MCSCF} } +@end macro + +@macro scf +{@acronym{SCF} } +@end macro + +@macro cas +{@acronym{CAS-SCF} } +@end macro + +@macro rhf +{@acronym{RHF} } +@end macro + +@macro ci +{@acronym{CI} } +@end macro + +Output files of Gaussian, Molpro and @gamess can be read to build the wave function files. +A major constraint is to realize @emph{single point} a calculation. + +@subsection Using Gaussian + +In the Gaussian input file, use the keywords @code{GFPRINT} and @code{pop=Full}. +In the case of @cas wave functions, use the @code{#p} keyword and the @code{SlaterDet} +attribute of the @code{CAS} keyword. + +@subsection Using Molpro + +Use the following options in the Molpro input file: +@itemize @bullet +@item @code{print,basis;} +@item @code{gprint,civector;} +@item @code{gprint,orbital;} +@item @code{gthresh,printci=0.;} for @mcscf calculations +@end itemize + +An @rhf calculation is mandatory before any @mcscf calculation. Be sure to +print @emph{all} molecular orbitals using the @code{orbprint} keyword. + +@subsection Using @gamess + +For @mcscf calculations, first compute the @mcscf single-point wave function +with the @acronym{GUGA} algorithm. Then, put the the @mcscf orbitals in the +@gamess input file, and run a single-point @acronym{GUGA} @ci calculation with +the following keywords: +@itemize @bullet +@item +@code{PRTTOL=0.0} in the @code{$GUGDIA} group +@item +@code{NPRT=2} in the @code{$CIDRT} group +@item +@code{PRTMO=.T.} in the @code{$GUESS} group +@end itemize + + diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f index 6ae4787..ba007f8 100644 --- a/src/eplf_function.irp.f +++ b/src/eplf_function.irp.f @@ -4,9 +4,10 @@ BEGIN_PROVIDER [ real, eplf_gamma ] ! Value of the gaussian for the EPLF END_DOC include 'constants.F' - real :: eps + real :: eps, N + N = 1. eps = -real(dlog(tiny(1.d0))) - eplf_gamma = (4./3.*pi*density_value_p)**(2./3.) * eps + eplf_gamma = (4./(3.*N)*pi*density_value_p)**(2./3.) * eps END_PROVIDER BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ] @@ -40,13 +41,12 @@ BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ] do k=1,ao_num - if (abs(mo_coef(k,i)) /= 0.) then + if (mo_coef(k,i) /= 0.) then do l=1,ao_num t = mo_coef(k,i)*ao_eplf_integral_matrix(l,k) if (abs(ao_eplf_integral_matrix(l,k))>1.d-16) then do j=i,mo_num - mo_eplf_integral_matrix(j,i) = mo_eplf_integral_matrix(j,i) + & - t*mo_coef_transp(j,l) + mo_eplf_integral_matrix(j,i) += t*mo_coef_transp(j,l) enddo endif enddo @@ -70,43 +70,53 @@ END_PROVIDER END_DOC integer :: i, j - double precision :: thr - thr = 1.d-12 / eplf_gamma eplf_up_up = 0.d0 eplf_up_dn = 0.d0 PROVIDE mo_coef_transp - do j=1,elec_beta_num - do i=1,elec_beta_num + do j=1,mo_closed_num + do i=1,mo_closed_num eplf_up_up += 2.d0*mo_value_p(i)* ( & mo_value_p(i)*mo_eplf_integral_matrix(j,j) - & mo_value_p(j)*mo_eplf_integral_matrix(i,j) ) - enddo - - do i=elec_beta_num+1,elec_alpha_num - eplf_up_up += mo_value_p(i)* ( & - mo_value_p(i)*mo_eplf_integral_matrix(j,j) - & - mo_value_p(j)*mo_eplf_integral_matrix(i,j) ) + eplf_up_dn += 2.d0*mo_value_p(j)**2 * & + mo_eplf_integral_matrix(i,i) enddo enddo - do j=elec_beta_num+1,elec_alpha_num - do i=1,elec_alpha_num - eplf_up_up += mo_value_p(i)* ( & - mo_value_p(i)*mo_eplf_integral_matrix(j,j) - & - mo_value_p(j)*mo_eplf_integral_matrix(i,j) ) - enddo - enddo + integer :: k,l,m,n,p + double precision :: ckl + do k=1,det_num + do l=1,det_num + + ckl = det_coef(k)*det_coef(l) + + do p=1,2 + do m=1,elec_num_2(p)-mo_closed_num + j = det(m,k,p) + do n=1,elec_num_2(p)-mo_closed_num + i = det(n,l,p) + eplf_up_up += ckl*mo_value_p(i)* ( & + mo_value_p(i)*mo_eplf_integral_matrix(j,j) - & + mo_value_p(j)*mo_eplf_integral_matrix(i,j) ) + enddo + enddo + enddo + + do m=1,elec_beta_num-mo_closed_num + j = det(m,k,2) + do n=1,elec_alpha_num-mo_closed_num + i = det(n,k,1) + eplf_up_dn += ckl * ( mo_value_p(i)**2 * mo_eplf_integral_matrix(j,j) & + + mo_value_p(j)**2 * mo_eplf_integral_matrix(i,i) ) + enddo + enddo - do j=1,elec_beta_num - do i=1,elec_alpha_num - eplf_up_dn += mo_value_p(i)**2 * & - mo_eplf_integral_matrix(j,j) enddo enddo - eplf_up_dn = 2.d0*eplf_up_dn + END_PROVIDER