EDF: Computing electron number probability distribution functions in real space from molecular wave functions

Computer Physics Communications 178 (2008) 621–634 www.elsevier.com/locate/cpc

EDF: Computing electron number probability distribution functions in real space from molecular wave functions ✩ E. Francisco ∗ , A. Martín Pendás, M.A. Blanco Departamento de Química Física y Analítica, Facultad de Química, Universidad de Oviedo, 33006 Oviedo, Spain Received 17 July 2007; accepted 29 November 2007 Available online 18 January 2008

Abstract

 3 Given an N-electron molecule and an exhaustive partition of the real space (R3 ) into m arbitrary regions Ω1 , Ω2 , . . . , Ωm ( m i=1 Ωi = R ), the edf program computes all the probabilities P (n1 , n2 , . . . , nm ) of having exactly n1 electrons in Ω1 , n2 electrons in Ω2 , . . . , and nm electrons (n1 + n2 + · · · + nm = N) in Ωm . Each Ωi may correspond to a single basin (atomic domain) or several such basins (functional group). In the later case, each atomic domain must belong to a single Ωi . The program can manage both single- and multi-determinant wave functions which are read in from an aimpac-like wave function description (.wfn) file (T.A. Keith et al., The AIMPAC95 programs, http://www.chemistry. mcmaster.ca/aimpac, 1995). For multi-determinantal wave functions a generalization of the original .wfn file has been introduced. The new format is completely backwards compatible, adding to the previous structure a description of the configuration interaction (CI) coefficients and the determinants of correlated wave functions. Besides the .wfn file, edf only needs the overlap integrals over all the atomic domains between the molecular orbitals (MO). After the P (n1 , n2 , . . . , nm ) probabilities are computed, edf obtains from them several magnitudes relevant to chemical bonding theory, such as average electronic populations and localization/delocalization indices. Regarding spin, edf may be used in two ways: with or without a splitting of the P (n1 , n2 , . . . , nm ) probabilities into α and β spin components. Program summary Program title: edf Catalogue identifier: AEAJ_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEAJ_v1_0.html Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 5387 No. of bytes in distributed program, including test data, etc.: 52 381 Distribution format: tar.gz Programming language: Fortran 77 Computer: 2.80 GHz Intel Pentium IV CPU Operating system: GNU/Linux RAM: 55 992 KB Word size: 32 bits Classification: 2.7 External routines: Netlib Nature of problem: Let us have an N-electron molecule and define an exhaustive partition of the physical space into m three-dimensional regions. The edf program computes the probabilities P (n1 , n2 , . . . , nm ) ≡ P ({np }) of all possible allocations of n1 electrons to Ω1 , n2 electrons to Ω2 , . . . , and nm electrons to Ωm , {np } being integers. ✩ This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect (http://www.sciencedirect. com/science/journal/00104655). * Corresponding author. Tel.: +(34) 985 10 30 39; fax: +(34) 985 10 3125. E-mail address: [email protected] (E. Francisco).

0010-4655/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2007.11.015

622

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

Solution method: Let us assume that the N-electron molecular wave function, Ψ (1, N), is a linear combination of M Slater determinants,  C ψr (1, N). Calling Srs Ψ (1, N) = M r r Ωk the overlap matrix over the 3D region Ωk between the (real) molecular spin-orbitals (MSO) in s ), edf finds all the P ({n })’s by solving the linear system r r ψr (χ1 , . . . χN ) and the MSOs in ψs , (χ1s , . . . , χN p m   M m   n   tk k P ({np }) = Cr Cs det tk Srs (1) Ωk , {np }

r,s

k

k

where tm = 1 and t1 , . . . , tm−1 are arbitrary real numbers. Restrictions: The number of {np } sets grows very fast with m and N, so that the dimension of the linear system (1) soon becomes very large. Moreover, the computer time required to obtain the determinants in the second member of Eq. (1) scales quadratically with M. These two facts limit the applicability of the method to relatively small molecules. Unusual features: Most of the real variables are of precision real*16. Running time: 0.030, 2.010, and 0.620 seconds for Test examples 1, 2, and 3, respectively. References: [1] A. Martín Pendás, E. Francisco, M.A. Blanco, Faraday Discuss. 135 (2007) 423–438. [2] A. Martín Pendás, E. Francisco, M.A. Blanco, J. Phys. Chem. A 111 (2007) 1084–1090. [3] A. Martín Pendás, E. Francisco, M.A. Blanco, Phys. Chem. Chem. Phys. 9 (2007) 1087–1092. [4] E. Francisco, A. Martín Pendás, M.A. Blanco, J. Chem. Phys. 126 (2007) 094102. [5] A. Martín Pendás, E. Francisco, M.A. Blanco, C. Gatti, Chemistry: A European Journal 113 (2007) 9362–9371. © 2007 Elsevier B.V. All rights reserved. PACS: 31.15.-p; 31.10.+z; 31.15.Ar; 02.70.-c Keywords: Quantum theory of atoms in molecules; Electron probability distribution; Molecular wave function; Chemical bonding theory

1. Introduction According to the principles of Quantum Mechanics, the multi-electron wave function, Ψ (1, . . . , N ), contains all the information relevant to a molecular system. Actually, the information contained in Ψ (1, . . . , N) is richer than is usually required in most quantum chemical discussions. For instance, since the standard non-relativistic molecular Hamiltonian involves only one- and twoparticle terms, the second order density matrix, γ 2 (x1 , x2 ; x1 , x2 ), obtained after integrating Ψ Ψ over all the spin and spatial coordinates of electrons 2 to N , is enough to extract from it all of the energetic properties that appear in Quantum Chemistry. In the above integration process, however, some information is lost. The answer to a question as simple as “which is the probability that exactly n out of the N electrons are located inside a region Ω of the physical space while the remaining N − n electrons are outside that region?” cannot be given in terms of γ 2 (x1 , x2 ; x1 , x2 ). In fact, to compute this probability, one needs the complete wave function Ψ (1, . . . , N ) [2], using



2 N P (n) = (2) dx1 · · · dxn Ψ (1, . . . , N) dxn+1 · · · dxN , n Ω

Ω

where

Ω

= R3

P (n) =

− Ω, or the diagonal density matrices of orders from n to N [3], using

N−n 1  (−1)i γ n+i (x1 , . . . , xn+i ) dx1 · · · dxn+i . n! i! i=0

(3)

Ω

In more general terms, let us assume that we partition exhaustively the physical space R3 into m arbitrary disjoint domains Ω1 , Ω2 , . . . , Ωm , and we want to know the probability that exactly n1 electrons are within Ω1 , n2 electrons within Ω2 , . . . , and nm electrons (n1 + n2 + · · · + nm = N ) within basin Ωm , where all ni ∈ N . To obtain these probabilities or multi-variate electron distribution functions (EDF) P (n1 , n2 , . . . , nm ), one needs again Ψ (1, . . . , N): P (n1 , n2 , . . . , nm ) = N !Λ Ψ Ψ dx1 · · · dxN , (4) D

where D is a multi-dimensional domain in which the first n1 electrons are integrated over Ω1 , the second n2 electrons over Ω2 , . . . , and the last nm electrons over Ωm , and N !Λ = N !/(n1 !n2 ! · · · nm !) is a combinatorial factor that accounts for electron indistinguishability. The 3D domains of these integrations can be arbitrary. However, some of these domains, such as the loge regions of Daudel [2], the atomic basins of the Quantum Theory of Atoms in Molecules (QTAM) [4], or the core, valence, and lone pair regions derived from the Electron Localization Function (ELF) [5,6], have been repeatedly used due to their relevant physical and chemical meanings. When using QTAM atomic basins, a partition of the N electrons of the molecule that assigns a given

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

623

number of electrons (including possibly 0) to each of these regions will be called a real space resonance structure (RSRS). There are NS = (N + m − 1)!/[N !(m − 1)!] of these for a given N, m pair. With the notation S(n1 , n2 , . . . , nm ) ≡ S({np }), or simply (n1 , n2 , . . . , nm ) ≡ {np }, we label the resonance structure having n1 electrons in Ω1 , n2 electrons in Ω2 , . . . , and nm electrons in Ωm . The importance of the computation of P (n1 , n2 , . . . , nm ) for all the RSRSs stems from the fact that they are the key quantities in a later calculation of the average electron population in a region ΩA [4], the delocalization index (DI) between two fragments A and B [7–11], the multi-center bond orders according to particular definitions [12–15], and other chemical bonding indices [16–19]. When Ψ (1, . . . , N) is a single Slater determinant, Eq. (2) may be worked out to give efficient formulas to compute P (n) [20–22]. However, the generalization of these formulas to any number of regions (m > 2) and/or correlated wave functions (i.e. Ψ (1, . . . , N ) is expanded as a linear combination of Slater determinants) is far from trivial. Recently [19], we have developed such a method to obtain the P (n1 , n2 , . . . , nm ) values simultaneously for the full set (NS ) of RSRSs, a task which has proven to be much easier than the computation of the probability of a single RSRS. In the following section we will present the basic equations of this method and extend it to include also the spin variables of the electrons, i.e. to split each P (n1 , n2 , . . . , nm ) into its spin-dependent probabilities. As we will see, all that is needed as input in the method are the overlap integrals between the molecular orbitals of the system within all the Ωi (i = 1, . . . , m) domains. The structure and input of the edf program implementing the method are described in Section 3. Finally, the results of a couple of test examples, for non-correlated and correlated wave functions, respectively, are given in Section 4. 2. Electron number probability distribution functions 2.1. Formal structure of Ψ Ψ We will assume a N -electron molecule with a multi-electron wave function Ψ (1, . . . , N) given by a linear combination of M Slater determinants Ψ (1, . . . , N) =

M 

Cr ψr (1, . . . , N) =

r=1

M 

Cr (N!)−1/2 det[ψ r ],

(5)

r=1

where ψ r is the (N × N ) matrix ⎡ χ r (1) . . . χ r (N ) ⎤ ⎢ ψr = ⎣

1

1

.. .

..

χNr (1)

.

...

.. .

⎥ ⎦,

(6)

χNr (N )

where χir are real orthogonal and normalized molecular spin-orbitals (MSO). The product Ψ Ψ = |Ψ |2 that appears in Eq. (4) is then given as |Ψ |2 =

M 

  Cr Cs (N!)−1 det ψ r ψ ts ,

(7)

r,s

where we have used the property det[ψ r ] det[ψ ts ] = det[ψ r ψ ts ]. The (i, j ) element of the (N × N ) product matrix ψ r ψ ts is given by     (ψ r )ik (ψ s )j k = χir (k)χjs (k). ψ r ψ ts ij = (8) k

k

k

k

 If a1 , . . . , aN and bk are column vectors, with ai = k αk bk , where αk are real or complex numbers, a useful property of determinants is     αk bk , . . . aN = αk det[a1 , . . . , bk , . . . , aN ]. det a1 , . . . , (9) If we apply successively Eq. (9) to columns 1, 2, . . . , N of ψ r ψ ts we have ⎡ χ r (k )χ s (k ) . . . χ r (k )χ s (k ) ⎤ 1 1 1 1 1 N N N N N     ⎥ ⎢ . .. t . .. .. det ψ r ψ s = ... det ⎣ ⎦. . k1 =1

kN =1

χNr (k1 )χ1s (k1 )

...

(10)

χNr (kN )χNs (kN )

In Eq. (10), if km = kn , the determinant can be put as χms (km )χns (kn ) times another determinant in which columns m and n are equal, and is thus equal to zero. This means that only the terms with all the kj (j = 1, . . . , n) different in Eq. (10) contribute to

624

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

det[ψ r ψ ts ]. The final expression for det[ψ r ψ ts ] is then      det χir (kj )χjs (kj ) , det ψ r ψ ts =

(11)

{kj }∈SN

where SN is the set of N! permutations of {1, 2, . . . , N }, {kj } is one of these permutations, and det[χir (kj )χjs (kj )] is an abbreviated form given to the determinant in the second member of Eq. (10), on the understanding that all the kj are different. 2.2. A linear system to compute P (n1 , n2 , . . . , nm ) Using Eq. (7) in Eq. (4) we have M      Cr Cs Prs {np } , P (n1 , . . . , nm ) ≡ P {np } =

(12)

r,s

where   Prs {np } = Λ



  det ψ r ψ ts dx1 · · · dxN ,

(13)

D

and using Eq. (11) in Eq. (13) we obtain      Prs {np } = Λ det χir (kj )χjs (kj ) dx1 · · · dxN .

(14)

{kj }∈SN D

Carrying out the summation over the spin variables and the spatial integration in Eq. (14) we have ⎡ Srs (k ) . . . Srs (k ) ⎤ 11 1 1N N    ⎢ ⎥ . .. . .. .. det ⎣ Prs {np } = Λ ⎦, . {kj }∈SN rs Srs N1 (k1 ) . . . SN N (kN ) where

Srs ij (kj ) =

⎧ rs (SΩ )ij ⎪ ⎪ ⎪ rs1 ⎪ ⎨ (SΩ2 )ij

if kj ∈ {1, . . . , n1 }, if kj ∈ {n1 + 1, . . . , n1 + n2 },

⎪ ⎪ ⎪ ⎪ ⎩

.. . if kj ∈ {N − nm + 1, . . . , N }

(Srs Ωm )ij

(15)

(16)

and (Srs Ωk )ij is the overlap integral between the ith MSO of the Slater determinant ψr and the j th MSO of the Slater determinant ψs within the Ωk domain, which can be computed as,  rs   r s SΩk ij = δ σi , σj (17) χir (r)χjs (r) dr, Ωk

by separating the spatial orbitals χir (r) and χjs (r), and integrating out the spin parts so that δ(α, α) = δ(β, β) = 1 and δ(α, β) = δ(β, α) = 0. Each Prs ({np }) is thus a sum of N ! determinants made of overlap integrals between the MSOs in the Slater determinants ψr and ψs . In each of these determinants, the n1 columns for which kj ∈ {1, . . . , n1 } are overlap integrals extended over Ω1 , the n2 columns for which kj ∈ {n1 + 1, . . . , n1 + n2 } are overlap integrals extended over Ω2 , . . . , and the nm columns for which kj ∈ {N − nm + 1, . . . , N } are overlap integrals extended over Ωm . To perform the summation that appears in the second member of Eq. (15) is not possible except for very small values of N . The interesting point is, however, that the following identity  m   n   nm rs 1 t1 × · · · × tm Prs {np } = det tk S Ωk (18) {np }

k=1

nm is valid for arbitrary values of t1 , t2 , . . . , tm [22]. Multiplying both members of Eq. (12) by t1n1 × · · · × tm and performing the summation over all possible {np } structures we arrive to  m M    n   nm t1 1 × · · · × tm P {np } = Cr Cs det tk Srs (19) Ωk . {np }

r,s

k=1

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

625

The left-hand side of Eq. (19) contains NS terms, i.e. as many as the number of {np } structures with n1 + · · · + nm = N . The general procedure to obtain the NS unknowns P (n1 , . . . , nm ) is now easy to grasp. Giving NS values to t1 , t2 , . . . , and tm−1 (remember that nm in the first member of Eq. (19) are computed for each of the NS structures, and the second tm = 1), the coefficients t1n1 t2n2 . . . tm member from the N × N overlap matrices Srs Ωk . A linear system of NS equations in the NS unknowns P (n1 , . . . , nm ) is obtained. Numerical instabilities aside, their solutions do not depend on the tk ’s, and it is immaterial which values we give to them from a strictly mathematical point of view. However, the linear system is ill conditioned and some care is necessary to choose the values of these parameters. In practice, by randomly selecting their values from the [0, 1] interval we avoid the linear dependence between the NS equations. Using this method of selection of the tk ’s, the same P (n1 , . . . , nm ) probabilities were obtained once and again in the two test examples given in Section 4 and also in other cases previously studied [19]. The rapid increase of NS with N and m prevents, however, that the linear system (19) may be exactly solved except for relatively small molecules [19]. When the exact solution is not possible an approximate method can still be used that consists in restricting the allowed {np } structures in Eq. (19) by avoiding those that do assign less (more) than a pre-assigned minimum (maximum) number of electrons to one or more of the m domains. This strategy reduces considerably the dimension of the linear system. Nevertheless, after solving it, we should test the consistency of the P (n1 , . . . , nm ) values by checking that they are all positive and add to a number only marginally different from 1.0. 2.3. Probabilities including spin variables The ordering of the MSOs in each Slater determinant ψr is arbitrary. However, since each MSO has spin α or β, it is always possible to reorder the set χir (i = 1, . . . , N ) within each ψr such that the Nα MSOs of type α are the first in the list and the Nβ MSOs of type β are the last ones (Nα + Nβ = N ). Moreover, Nα (and thus also Nβ ) is the same in all ψr ’s, so that all the overlap matrices Srs Ωk have the same block-diagonal structure:  rs,α S Ωk 0 rs S Ωk = (20) . rs,β 0 S Ωk This ordering process is computationally very cheap compared to the evaluation of the determinants in the second member of Eq. (19). It also allows us to generalize very easily Eq. (19) to include the spin variables of the electrons, i.e. to split each P (n1 , n2 , . . . , nm ) into its spin-dependent probabilities. Let us assume that we want to know the probability that exactly nα1 electrons with spin α are within Ω1 , nα2 electrons with spin α are within Ω2 , . . . , nαm electrons with spin α (nα1 + · · · + nαm = Nα ) β β β are within Ωm , and simultaneously, n1 electrons with spin β are within Ω1 , n2 electrons with spin β are within Ω2 , . . . , and nm β β rs,β rs,α electrons with spin β (n1 + · · · + nm = Nβ ) are within Ωm . Calling S˜ Ωk and S˜ Ωk the N × N matrices defined, respectively, as   rs,α   0 S 0 0 rs,β ˜ ˜Srs,α = Ωk and SΩk = (21) rs,β , Ωk 0 0 0 S Ωk β β rs ˜ rs,α ˜ rs,β and defining (i1 , . . . , im , im+1 , . . . , i2m ) ≡ {ip } = (nα1 , . . . , nαm , n1 , . . . , nm ), Trs k = SΩk , and Tm+k = SΩk (k = 1, . . . , m), these s probabilities, P ({ip }), may be obtained by solving the following linear system:  2m M  i i    i2m s  rs 1 2 (22) t1 t2 . . . t2m P {ip } = Cr Cr det tk Tk , {ip }

r,s

k=1 β

where t2m = 1 and t1 , . . . , t2m−1 are arbitrary real numbers. Eq. (22) is a linear system with NSα × NS equations and unknowns (NSσ = (Nσ + m − 1)!/[Nσ !(m − 1)!] (σ = α, β)) that can be solved as Eq. (19). As in that case, one can reduce the dimension of the above linear system by neglecting those structures that, based on chemical wisdom, would very likely have a negligible probability in case that the linear system were exactly solved. After solving Eq. (22), each P (n1 , . . . , nm ) spinless probability is easily obtained by adding those P s ({ip })’s for which i1 + im+1 = n1 , i2 + im+2 = n2 , . . . , and im + i2m = nm . 2.4. Chemical bonding indices from probabilities m-domain probabilities allow us to recover important chemical bonding indices. The most obvious of these indices is the average electronic population of each domain, ni , given by the standard definition of the average of a function over a discrete distribution:    ni  = (23) ni × P {np } , {np }

or

nσi ,

the average number of electrons in the domain i with spin σ :   σ   nσi = ni × P s {ip } .



{ip }

(24)

626

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

The delocalization index (DI) between two domains i and j (δ i,j ), which is a measure of the number of pairs of electrons shared between both domains [7–11,23,24], is usually obtained from γ 2 (x1 , x2 ) ≡ γ 2 (x1 , x2 ; x1 , x2 ), and the electron density ρ(x) as    γ 2 (x1 , x2 ) − ρ(x1 )ρ(x2 ) dx1 dx2 . δ i,j = −2 (25) σ1

σ2

Ωi Ωj

However, it can also be obtained from the P ({np }) probabilities as   δ i,j = −2 ni nj  − ni nj    σ σ    i,j = −2 (δ )σ1 σ2 , ni 1 nj 2  − nσi 1 nσj 2  = σ1

σ2

σ1

(26) (27)

σ2

where ni  and nσi  are given by Eqs. (23) and (24), respectively, and ni nj  and nσi 1 nσj 2  are obtained by using    ni nj  = ni nj × P {np }

(28)

{np }

and

  σ1 σ2   nσi 1 nσj 2 = ni nj × P {ip } .



(29)

{ip }

Eqs. (23) and (28) can be generalized to provide the average value of any function of the ni ’s.         f {np } = f {np } × P {np } .

(30)

{np }

In this way, the K-fragment DI’s, defined as [12–15] ! k  (−2)k−1  i,j,...,k δ nx − nx  , = (k − 1)!

(31)

x=i

can be expanded in terms of average electron populations, ni , and average values of the form (30). 3. Structure and input of the program 3.1. Description of the input The program edf reads the input from Fortran unit 5 (stdin in Unix), which can be redirected from a file (e.g., edf < edf.inp), and from two files whose names (variables filedat and wfnfile) are read from stdin. The output of edf is written in Fortran unit 6 (stdout in Unix), which can be redirected to a file (e.g., edf < edf.inp > edf.out). The input is read in records, which usually correspond to a line. They are enumerated in what follows (a means character format, r means any real format, i means integer format). CAPITAL letters denote fixed keywords in the input, which is case insensitive and free format. Structure of stdin file (Fortran unit 5) (1) ioverlap,ispin, (i,i): ioverlap controls the format of the file from which the atomic overlap matrices (AOM) within the QTAM atomic regions are read in. If ioverlap=0, this format is that of our own program promolden [25]. If ioverlap=0, this file is assumed to have the format of the proaimv Fortran code [1]. If ispin=0 the probabilities including spin variables are not computed; i.e. the linear system actually solved is Eq. (19). If ispin=0, probabilities both including and not including spin variables are computed. (2) filedat, (a): Name of the file containing the AOMs. (3) wfnfile, (a): Name of the file containing the names, coordinates, basis set, and molecular wave function of the current molecule. (4) NGROUP ngroup, (a,i): ngroup is the number of disjoint fragments in which the molecule is divided. After this record, ngroup records are read, each one with the following structure: nfugrp(i),(ifugrp(j,i),j=1,nfugrp(i)),elpol(1,i),elpol(2,i), (i,i,...): nfugrp(i) is the number of atoms of the current molecule assigned to fragment i. ifugrp(j,i),j=1,nfugrp(i) specify which atoms belong to fragment i according to the order given in file wfnfile. Each atom can only belong to a single fragment. elpol(1,i) and elpol(2,i) are the minimum and maximum number of electrons assigned to fragment i. The value of elpol(1,i) must be greater than 0 and the value of elpol(2,i) must be smaller than N , the total number of electrons of the molecule.

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

627

(5) EPSDET epsdet, (a,r): epsdet has the following meaning. The product Cr Cs in Eq. (7) neglected if |Cr Cs | < epsdet. The default value is epsdet=0.0. (6) PROBCUT probcut, (a,r): In the output, only probabilities greater than probcut are written. The default value is probcut=0.0. (7) NDETS nactdet, (a,r): When a multi-determinantal wave function is read, only the first nactdet determinants are used in Eq. (5). The maximum value of nactdet is ndets, the actual number of Slater determinants defined in file wfnfile. This is also the default value for nactdet. It is thus convenient that the determinants in file wfnfile are ordered by decreasing order of the Cr coefficients. (8) # Text line, (a): # This is a comment line in the input file. (9) END, (a): Line signaling the end of the input file. Records from (1) to (3) are mandatory and must be given in the above order. Records from (4) to (8) can be given in any order, and all of them are optional. In case that record (4) is not given or ngroup  0, edf takes ngroup equal to the number of atoms of the molecule, with each fragment being a single atom with a minimum and maximum electronic population equal to 0 and N , respectively. Structure of wfnfile file (Fortran unit 1) This is the aimpac-like wave function description (.wfn) file [1]. For monodeterminantal wave functions, it contains the names, nuclear charges, atomic coordinates, basis set, and natural MOs (with their occupation numbers and orbital energies) of the current molecule. For multi-determinantal wave functions we have generalized the original .wfn file in a backwards compatible manner by adding the canonical MOs used in the CI expansion, and the coefficients and descriptions of the determinants. All the new information may be obtained from the output of standard quantum chemical packages. In this work we have used a locally modified version of the gamess code [26] to write this kind of .wfn file. In the last line of the original wfnfile the character variable label is read. When label(1:5) is MCSCF, ALDET, or GENCI, the wfnfile comes from a multi-determinantal wave function. The generalized wfnfile contents after this line are as follows: (1) label, (a): Comment line. (2) label, (a): Comment line. (3) coef(i,j),j=1,nprims, (5E16.8): Expansion coefficients of each canonical MO in terms of the nprims not normalized Cartesian Gaussian functions of the basis set. Records (2) and (3) must be repeated nmo times from i=nmo+1 to nmo+nmo. The elements coef(i,j) from i=1,nmo and j=1,nprims, corresponding to the natural MOs, have already been read from this file and for monodeterminantal wave functions are those actually used to compute the AOMs read from the filedat file (see below). However, for multi-determinantal wave functions the AOMs stored in filedat are computed using the canonical MOs. Actually, the coef array is not used at all within the edf program. However, as it is well known, the .wfn file is also used as input of the aimpac package and of our own promolden code that computes the AOMs. All of them need this array. For compatibility reasons, it is also read in within edf. (4) check, (a): Line indicating the end of the input of the expansion coefficients for canonical MOs. This variable must necessarily be equal to the string ‘END DATA’ in order for edf to continue. (5) label, (a): Comment line. (6) nelact,ndets,ncore,nact (i,i,i,i): nelact is the number of active (or valence) electrons, ndets is the number of Slater determinants in the wave function expansion (M in Eq. (5)), ncore is the number of core MOs, and nact is the number of active (or valence) MOs. (7) label, (a): Comment line. (8) cidet(j),j=0,nelact, (r,i,i,...): cidet(0) is the coefficient of the Slater determinant in the wave function expansion (Cr in Eq. (5)). The absolute value of int(cidet(j)) (j=0) specifies the canonical MO assigned to the j th row of the active-space (non-core) part of the Slater determinant, and its sign determines the spin of this electron (positive for α and negative for β). Record (8) must be repeated ndets times. Structure of filedat file (Fortran unit 18) This file contains the AOMs between all the molecular orbitals (MO) χi within all the QTAM atomic basins, χi |χj a . The overlaps extended  to the Ωk (k=1,ngroup) regions in which the molecule has been divided are computed from the AOMs by using χi |χj Ωk = a∈Ωk χi |χj a . The structure of filedat depends on the value of the ioverlap parameter read from stdin. When ioverlap=0, this structure faithfully corresponds with that given in the output of the proaimv code of the aimpac package [1], and when ioverlap=0 it corresponds to that written by our promolden code [25]. ioverlap=0

(1) iaom, (i): Index of the atom for which the next record is read.

628

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

(2) ((aom(iaom,m,j),m=1,j),j=1,nmo), (6(1X,E16.10)): aom(iaom(i),m,j) is the AOM element χm |χj iaom , where nmo is the number of MOs. The above two records must be read ncenter times, where ncenter is the number of atoms in the molecule. Records (1) and (2) must be repeated ncenter times with the iaom variable taking all of its possible values (0, 1, . . . , ncenter). ioverlap=0

(1) iaom, (i): Index of the atom for which the next record is read. Each of these nmo records (m=1,nmo) has the following structure. (2) (aom(iaom,j,m),j=1,m), (8F10.6): aom(iaom(i),j,m) is the AOM element χj |χm iaom . Records (1) and (2) must be repeated ncenter times. 3.2. Main algorithm Following is a pseudo-code of the main algorithm: (1) (2) (3) (4) (5) (6) (7) (8)

(9) (10)

Read the first three records of stdin file. Read wfnfile file. Read filedat file.  centers Modify the AOM of the last center such that all χj |χm a = δj m . This is done to slightly improve the accuracy of a calculations. Read records from (4) to (9) of stdin file. Compute the χi |χj Ωk overlaps (i,j=1,nmo) from the AOMs. Determine all possible real space resonant structures ({np } when ispin=0 or {ip } when ispin=0). Solve the linear system (19) (ispin=0) or (22) (ispin=0). This is the core of edf. This task is divided in several steps. (a) Renormalize the wave function if necessary. i2m nm (b) Obtain the coefficients t1n1 × · · · × tm of Eq. (19) or t1i1 t2i2 . . . t2m of Eq. (22) by randomly selecting the ti values from the [−1, +1] interval. (c) Determine the second members of Eq. (19) or Eq. (22). (d) Solve Eq. (19) or Eq. (22). (e) Output results of probabilities. In case that ispin=0, compute the spinless P ({np }) probabilities from the P ({ip }) probabilities. Compute and output average electronic populations of the fragments and 2-fragment and 3-fragment delocalization indices.

3.3. Data structures The relevant data used in edf are stored in standard Fortran arrays: aom(m,i,j) contains the overlap integral extended to the atomic basin of atom m between the ith and the j th canonical MOs. sg(k,i,j) contains the overlap integral extended Ωk between the same pair of MOs. The electronic populations of the different real space resonant structures, {np }, are stored in resnc, resnc(j, i) being the number of electrons in fragment i for structure j . This array is obtained from the input data contained in the ifugrp, nfugrp, and elpol arrays, where nfugrp(i) is the number of atoms of fragment i, ifugrp(m, i) the order number of atom m of fragment i according to the order given in the input file wfnfile, and elpol(1, i) (elpol(2, i)) the minimum (maximum) number of electrons of fragment i. For each Slater determinant in the wave function expansion (Eq. (5)), cidet(0) is the expansion coefficient Cr , and int(cidet(j)) (j=0) specifies the MSO assigned to the (2 ncore+j)th row of this determinant: its absolute value identifies the active canonical MO corresponding to its spatial part, its sign describing the spin part (+ for α, − for β). The tpow array contains the ti values that appear in Eq. (19) (or Eq. (22) if the spin variables are included in the calculation), n n tpow(j,k) being the value corresponding to fragment k for structure j . From it and the resnc array, the products t1 1 × · · · × tmm are  rs,β rs,β rs,α rs,α  computed and stored in acoef. The ovea (oveb) and oveaa (ovebb) arrays contain the SΩk (SΩk ) and k tk SΩk ( k tk SΩk ) values. From oveaa, ovebb, and the cidet(0) coefficients, the second member of Eq. (19) or Eq. (22) is evaluated and stored in the prob array. After solving the linear system (19) or (22), prob also contains the P ({np }) or P s ({ip }) probabilities. Finally, from all the MSOs that participate in the wave function expansion (Eq. (5)), the integer arrays ordia (α MSOs) and ordib (β MSOs) indicate which ones appear in the Slater determinant ψ r and ordja and ordjb which ones appear in ψ s , in the product ψ r ψ ts (Eq. (7)).

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

629

3.4. Program routine blocks Aside from the main program, edf, there are several blocks of routines, described by groups in the following list: Wave function reading: The wave function is read in from the file wfnfile by using a modified version of the AIMPAC rdwfn routine. Resonant structures probabilities: The snbasins and nbasins routines compute the probabilities for all the real space resonant structures, either including (snbasins) or not including (nbasins) the spin variables, by solving the linear system (22) or (19), respectively. The electronic populations of these structures are previously determined in snprobs and nprobs routines. The spatialp routine obtains the P ({np }) probabilities, that do not include spin variables, from the P s ({ip }) values, in which the spin variables are explicitly considered. Finally, the colord routine, called by nbasins and snbasins, re-orders the MSOs in the Slater determinant ψs such that all the overlap matrices Srs Ωk have the block-diagonal structure of Eq. (20). Delocalization indices: The delocind and sdelocind routines obtain the average number of electrons in the Ωk regions, and the localization and inter-fragment delocalization indices from the previously computed probabilities. Linear algebra tools: The daxpy, dgeco, dscal, dgefa, dswap, dgesl, dgedi, ddot, idamax, and dasum routines of the Netlib Repository (http://www.netlib.org) were modified to work in real*16 precision and incorporated in a single file called netlib.f. The dgeco routine factors a real*16 matrix A by Gaussian elimination and estimates the condition of the matrix. After dgeco is called, the dgesl routine is used to solve a linear system AX = B, and the dgedi routine to compute det[A]. Random numbers generation: The random and semilla routines are used in the generation of random numbers. random is a slightly modified version of the ran1 routine of Numerical Recipes [27], and returns a pseudo-random number x ∈ [0, 1]. The semilla routine generates a eight-digits integer (num) by calling the fdate routine of the system. Calling random with the argument −num, the sequence of generation of random numbers is reinitialized. String manipulation: setword and setint read a character word or an integer value from a character variable; they are used to allow free format in some of the lines of the input, so that some variables can be omitted. The routine leng computes the length of a word within a character variable, not counting blank spaces at the end. The uppcase routine transform alphabetical characters in a word to uppercase. The atoi routine transforms a character variable to an integer, and the logical function isdigit(c) returns .true. if the character*(1) variable c is a digit and .false. otherwise. CPU time routines: timer and atimer routines are used to determine the CPU time used in different tasks of the program. timer accumulates and prints out the elapsed times of a series of up to mprocess tasks (mprocess is defined as a parameter in the routine; currently mprocess=100). The atimer routine is an interface to the timer routines of different operating systems. It returns machine seconds at the calling time. The error routine is called by timer when an error usage of the latter occurs; error prints an error message. 4. Compiling and testing the program 4.1. Compilation The linear systems (19) and (22) are very ill-conditioned. If the matrix of coefficients of the first member of these linear systems is written as its singular value decomposition, A = U · diag(w) · Vt , the ratio r = wj,min /wj,max (i.e. the inverse of the condition number of the matrix A) is very small. We have observed in some cases r values comparable to the machine’s floating point precision in double precision ( 10−12 ) or  smaller. When this happens  the uncertainty in the derived probabilities P ({np }) (or P ({ip })) can be very big and the condition {np } P ({np }) = 1 (or {ip } P ({ip }) = 1) is not very well satisfied. For this reason, edf is distributed in real*16 and must be compiled with a Fortran compiler supporting this precision. The test examples presented below have used an edf version obtained with the ifort intel Fortran 95 compiler. A disadvantage of this real*16 version is that it is slower than a double precision version. However, unless the above mentioned uncertainties in the probabilities were observed, a double precision version can safely be used to speed up the running of the program. To do this, three changes in the source code are necessary: (1) change the single line of the implicit.inc file from implicit real*16 (a-h,o-z) to implicit double precision (a-h,o-z); (2) change the values of the constants defined in the constants.inc file from real*16 to double precision; and (3) substitute the routines daxpy, dgeco, dscal, dgefa, dswap, dgesl, dgedi, ddot, idamax, and dasum contained in the distributed netlib.f file by the original ones (in double precision) of the Netlib Repository (http://www.netlib.org). A relevant fact also related to the uncertainties in the P ({np }) (or P ({ip })) values is the random generation of the ti numbers that define the matrix of coefficients (A) in Eqs. (19) and (22). We have tested many different ways of choosing these numbers. Some of them led to r values too small to be acceptable. The final choice (−1 < ti < +1) that appears in the distribution files has proven to be the best one in tens of examples. However, there is no warranty that this is always the best possibility. We invite the user of edf to try other ways of choosing these numbers in case that the condition {np } P ({np }) = 1 fails to be fulfilled with the present implementation. For this purpose, the matrix element tpow(i,k), currently computed in the line that reads tpow(i,k)=two*random(idum)-one within the nbasins.f and snbasins.f files, must be substituted by some other value.

630

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

4.2. Test input file, Example 1, stdin=h2cas.inp 0 1 h2cas.aom h2cas.wfn ngroup 2 1 1 0 2 1 2 0 2 end

4.3. Summary of the test output file, Example 1, stdout=h2cas.out # # # # # # # # # # # # # # # # # #

------------------- EDF ------------------------ELECTRON PROBABILITY DISTRIBUTION FUNCTIONS (c) Evelio Francisco 2007 University of Oviedo ------------------------------------------------Atomic overlaps matrix = h2cas.aom Wave Function file = h2cas.wfn Number of electrons = 2 Number of molecular orbitals = 2 Number of determinants = 2 Number of active electrons = 2 Number of core orbitals = 0 Number of active orbitals = 2 This is a Multideterminantal wavefunction Atomic overlap matrix read in with PROMOLDEN format Analysis including SPIN probabilities Reading atomic overlap matrix for atom 1 Reading atomic overlap matrix for atom 2

Cartesian coordinates (bohr) and nuclear charges H 1 0.00000000 0.00000000 -0.70055925 1.0 H 2 0.00000000 0.00000000 0.70055925 1.0 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ FRAGMENT 1 FORMED BY 1 ATOMS : 1 (MinElec,MaxElec) : 0 2 FRAGMENT 2 FORMED BY 1 ATOMS : 2 (MinElec,MaxElec) : 0 2 # Number of probabilities to be computed =

4

Linear System Solver: Netlib DGECO routine RECIPROCAL OF RCOND VALUE IS 0.619643004847E-01 # M-BASINS ELECTRON NUMBER PROBABILITY DISTRIBUTION INCLUDING SPIN # -----------------------------------------------------------------------# NUMBER OF GROUPS = 2 # TOTAL NUMBER OF PROBABILITIES = 4 # Gi(a) Gi(b) ARE THE NUMBER OF ALPHA AND BETA ELECTRONS IN GROUP i #-----------------------------------------------------------------------# Probability G1(a) G1(b) G2(a) G2(b) G3(a) G3(b) # 0.2083118319895712 0 0 1 1 # 0.2916881680104288 0 1 1 0 # 0.2916881680104288 1 0 0 1 # 0.2083118319895712 1 1 0 0 #-----------------------------------------------------------------------# 1.0000000000000000 <--- SUM, 4 PROBABILITIES > 0.0000000000E+00 # 1.0000000000000000 <--- TOTAL SUM #-----------------------------------------------------------------------Average populations and delocalization indices

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

delta_( 2 1)_{alpha,alpha) delta_( 2 1)_{alpha,beta ) delta_( 2 1)_{beta,alpha) delta_( 2 1)_{beta,beta )

= = = = = = = =

0.5000000000 0.5000000000 0.5000000000 0.5000000000 0.5000000000 -0.0833763360 -0.0833763360 0.5000000000

# M-BASINS ELECTRON NUMBER PROBABILITY DISTRIBUTION #-----------------------------------------------------------------------# NUMBER OF GROUPS = 2 # TOTAL NUMBER OF PROBABILITIES = 3 #-----------------------------------------------------------------------# Probability n1 n2 n3... # 0.2083118319895712 0 2 # 0.5833763360208577 1 1 # 0.2083118319895712 2 0 #-----------------------------------------------------------------------# 1.0000000000000000 <--- SUM, 3 PROBABILITIES > 0.0000000000E+00 # 1.0000000000000000 <--- TOTAL SUM #-----------------------------------------------------------------------Average populations and delocalization indices = 1.0000000000 = 1.0000000000 delta_( 2 1) = 0.8332473280

4.4. Test input file, Example 2, stdin=cas-lih.inp 0 0 cas-lih.aom cas-lih.wfn probcut 0.001 ngroup 2 1 1 0 4 1 2 0 4 end

4.5. Summary of the test output file, Example 2, stdout=cas-lih.out # # # # # # # # # # # # # # # # # #

------------------- EDF ------------------------ELECTRON PROBABILITY DISTRIBUTION FUNCTIONS (c) Evelio Francisco 2007 University of Oviedo ------------------------------------------------Atomic overlaps matrix = cas-lih.aom Wave Function file = cas-lih.wfn Number of electrons = 4 Number of molecular orbitals = 6 Number of determinants = 65 Number of active electrons = 4 Number of core orbitals = 0 Number of active orbitals = 6 This is a Multideterminantal wavefunction Atomic overlap matrix read in with PROMOLDEN format Analysis NOT including SPIN probabilities Reading atomic overlap matrix for atom 1 Reading atomic overlap matrix for atom 2

NEW PROBCUT VALUE = 0.100000E-02

631

632

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

Cartesian coordinates (bohr) and nuclear charges LI 1 0.00000000 0.00000000 0.24011857 3.0 H 2 0.00000000 0.00000000 3.25988143 1.0 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ FRAGMENT 1 FORMED BY 1 ATOMS : 1 (MinElec,MaxElec) : 0 4 FRAGMENT 2 FORMED BY 1 ATOMS : 2 (MinElec,MaxElec) : 0 4 # Number of probabilities to be computed =

5

Linear System Solver: Netlib DGECO routine RECIPROCAL OF RCOND VALUE IS 0.114533879175E-02 # M-BASINS ELECTRON NUMBER PROBABILITY DISTRIBUTION #-----------------------------------------------------------------------# NUMBER OF GROUPS = 2 # TOTAL NUMBER OF PROBABILITIES = 5 #-----------------------------------------------------------------------# Probability n1 n2 n3... # 0.0064374595987544 1 3 # 0.8906560161532174 2 2 # 0.1010556187347743 3 1 # 0.0018421988499153 4 0 #-----------------------------------------------------------------------# 0.9999912933366613 <--- SUM, 4 PROBABILITIES > 0.1000000000E-02 # 1.0000000000000000 <--- TOTAL SUM #-----------------------------------------------------------------------Average populations and delocalization indices = 2.0982851435 = 1.9017148565 delta_( 2 1) = 0.2104734619

4.6. Test input file, Example 3, stdin=hcn.inp 0 0 hcn.aom hcn.wfn probcut 0.001 ngroup 3 1 1 0 14 1 2 0 14 1 3 0 14 end

4.7. Summary of the test output file, Example 3, stdout=hcn.out # # # # # # # # # # # # # #

------------------- EDF ------------------------ELECTRON PROBABILITY DISTRIBUTION FUNCTIONS (c) Evelio Francisco 2007 University of Oviedo ------------------------------------------------Atomic overlaps matrix = hcn.aom Wave Function file = hcn.wfn Number of electrons = 14 Number of molecular orbitals = 7 This is a RHF wavefunction Atomic overlap matrix read in with PROMOLDEN format Analysis NOT including SPIN probabilities Reading atomic overlap matrix for atom 1 Reading atomic overlap matrix for atom 2

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

# Reading atomic overlap matrix for atom

3

NEW PROBCUT VALUE = 0.100000E-02 Cartesian coordinates (bohr) and nuclear charges H 1 0.00000000 0.00000000 -1.77539511 1.0 C 2 0.00000000 0.00000000 0.22357354 6.0 N 3 0.00000000 0.00000000 2.35182156 7.0 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ FRAGMENT 1 FORMED BY 1 ATOMS : 1 (MinElec,MaxElec) : 0 14 FRAGMENT 2 FORMED BY 1 ATOMS : 2 (MinElec,MaxElec) : 0 14 FRAGMENT 3 FORMED BY 1 ATOMS : 3 (MinElec,MaxElec) : 0 14 # Number of probabilities to be computed =

120

Linear System Solver: Netlib DGECO routine RECIPROCAL OF RCOND VALUE IS 0.136145722119E-10 # M-BASINS ELECTRON NUMBER PROBABILITY DISTRIBUTION #-----------------------------------------------------------------------# NUMBER OF GROUPS = 3 # TOTAL NUMBER OF PROBABILITIES = 120 #-----------------------------------------------------------------------# Probability n1 n2 n3... # 0.0023058610654656 0 3 11 # 0.0623842623916013 0 4 10 # 0.1291401792393889 0 5 9 # 0.1133328106266282 0 6 8 # 0.0526176225327567 0 7 7 # 0.0135370427960918 0 8 6 # 0.0018240431747734 0 9 5 # 0.0014421116993124 1 2 11 # 0.0720818003194062 1 3 10 # 0.1578190708901936 1 4 9 # 0.1446857315250721 1 5 8 # 0.0700556506683609 1 6 7 # 0.0188163798839127 1 7 6 # 0.0026538317927779 1 8 5 # 0.0207785507054511 2 2 10 # 0.0496276497870966 2 3 9 # 0.0487391232466751 2 4 8 # 0.0251327311132273 2 5 7 # 0.0071706398852352 2 6 6 # 0.0010737630625026 2 7 5 # 0.0016675699291450 3 3 8 # 0.0012355338543090 3 4 7 #-----------------------------------------------------------------------# 0.9981219601893837 <--- SUM, 22 PROBABILITIES > 0.1000000000E-02 # 1.0000000000000000 <--- TOTAL SUM #-----------------------------------------------------------------------Average populations and delocalization indices = 0.7861825349 = 4.8015879535 = 8.4122295116 delta_( 2 1) = 0.9158855687 delta_( 3 1) = 0.0840277492 delta_( 3 2) = 2.2922487320 delta_( 3 2 1) = 0.0354845620

633

634

E. Francisco et al. / Computer Physics Communications 178 (2008) 621–634

Acknowledgements Financial support from the Spanish Ministerio de Educación y Ciencia, Project No. CTQ2006-02976, is acknowledged. References [1] T.A. Keith, K.E. Laidig, P. Krug, J.R. Cheeseman, R.G.A. Bone, F.W. Biegler-König, J.A. Duke, T. Tang, R.F.W. Bader, The AIMPAC95 programs, the code is available at http://www.chemistry.mcmaster.ca/aimpac, 1995. [2] R. Daudel, The Fundamentals of Theoretical Chemistry, Pergamon Press, Oxford, 1968. [3] P. Ziesche, Int. J. Quantum Chem. 60 (1996) 149–162. [4] R.F.W. Bader, Atoms in Molecules, Oxford University Press, Oxford, 1990. [5] A.D. Becke, K.E. Edgecombe, J. Chem. Phys. 92 (1990) 5397. [6] B. Silvi, A. Savin, Nature 371 (1994) 683. [7] R.F.W. Bader, M.E. Stephens, J. Am. Chem. Soc. 97 (1975) 7391. [8] J.G. Ángyán, M. Loos, I. Mayer, J. Phys. Chem. 98 (1994) 5244–5248. [9] X. Fradera, M.A. Austen, R.F.W. Bader, J. Phys. Chem. A 103 (1999) 304–314. [10] X. Fradera, J. Poater, S. Simon, M. Duran, M. Solà, Theor. Chem. Acc. 108 (2002) 214–224. [11] J. Poater, M. Duran, M. Solà, B. Silvi, Chem. Rev. 105 (2005) 3911–3947. [12] K.C. Mundim, M. Giambiagi, M.S. de Giambiagi, J. Phys. Chem. 98 (1994) 6118–6119. [13] A.B. Sannigrahi, T. Kar, Chem. Phys. Lett. 299 (1999) 518–526. [14] A.B. Sannigrahi, T. Kar, J. Mol. Struct. (Theochem) 496 (2000) 1–17. [15] R.C. Bochicchio, R. Ponec, A. Torre, L. Lain, Theor. Chem. Acc. 105 (2001) 292–298. [16] A. Martín Pendás, E. Francisco, M.A. Blanco, Faraday Discuss. 135 (2007) 423–438. [17] A. Martín Pendás, E. Francisco, M.A. Blanco, Phys. Chem. Chem. Phys. 9 (2007) 1087–1092. [18] A. Martín Pendás, E. Francisco, M.A. Blanco, J. Phys. Chem. 111 (2007) 1084–1090. [19] E. Francisco, A. Martín Pendás, M.A. Blanco, J. Chem. Phys. 126 (2007) 094102. [20] E. Chamorro, P. Fuentealba, A. Savin, J. Comput. Chem. 24 (2002) 496–504. [21] A. Savin, in: K.D. Sen (Ed.), Reviews of Modern Quantum Chemistry, World Scientific Publishing, Singapore, 2002. [22] E. Cancès, R. Keriven, F. Lodier, A. Savin, Theor. Chem. Acc. 111 (2004) 373–380. [23] L. Rincón, J.E. Alvarellos, R. Almeida, J. Chem. Phys. 122 (2005) 214103. [24] L. Rincón, J.E. Alvarellos, R. Almeida, J. Chem. Phys. 122 (2005) 214104. [25] A. Martín Pendás, unpublished. [26] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.H. Jensen, S. Koseki, N. Matsunaga, K.A. Nguyen, S.J. Su, T.L. Windus, M. Dupuis, J.A. Montgomery, J. Comput. Chem. 14 (1993) 1347. [27] W.H. Press, B.P. Flannery, S.A. Teukolski, W.T. Vetterling, Numerical Recipes in Fortran 77: The Art of Scientific Computing, second ed., Cambridge University Press, Cambridge, 1992.