Improved gamma parameter

This commit is contained in:
Anthony Scemama 2010-04-29 17:15:07 +02:00
parent cc7163ebd3
commit a0eeaf84d3
5 changed files with 402 additions and 27 deletions

57
doc/eplf.texi Normal file
View File

@ -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

23
doc/ezfio.texi Normal file
View File

@ -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}.

221
doc/intro.texi Normal file
View File

@ -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.

64
doc/wf.texi Normal file
View File

@ -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

View File

@ -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