Electrostatic interaction energies with overlap effects from a localized approach

Chemical Physics Letters 445 (2007) 315–320 www.elsevier.com/locate/cplett

Electrostatic interaction energies with overlap effects from a localized approach Fazle Rob, Rafał Podeszwa *, Krzysztof Szalewicz Department of Physics and Astronomy, University of Delaware, Newark, DE 19716, United States Received 15 June 2007; in final form 20 July 2007 Available online 1 August 2007

Abstract An ab initio method is proposed for calculations of intermolecular forces which should be applicable to interactions of monomers containing several hundreds atoms and scale linearly with system size. This approach unifies asymptotic approximation to intermolecular forces with symmetry-adapted perturbation theory (SAPT) treatment. Inexpensive asymptotic formulas are applied to interactions of fragments of monomers which are far apart, whereas fairly expensive SAPT formulas are applied to fragments which are near each other, where overlap effects cannot be neglected. This Letter implements the method for the electrostatic interactions. The electrostatic energies are reproduced to 1% with about dozen or so atoms in near-range regions. Ó 2007 Elsevier B.V. All rights reserved.

1. Introduction In the past decade, symmetry-adapted perturbation theory (SAPT) has become a popular approach for ab initio calculations of intermolecular interaction energies [1,2]. SAPT is accurate but expensive in terms of computer resources and therefore it cannot be used for monomers with more than about 10 atoms. The more recently developed SAPT(DFT) method [3], i.e., SAPT based on the density functional theory (DFT) description of monomers, in its current form [4] can be applied routinely to compute interaction potentials for systems with monomers containing up to about 30 atoms each. Since SAPT(DFT) codes have already been fairly well optimized, applications of SAPT(DFT) to still larger systems require a completely novel approach. We propose a method that will take advantage of the fast, exponential decay of the overlap effects and of the simplicity of the asymptotic calculations of interaction energies. The method is implemented here

*

Corresponding author. E-mail address: [email protected] (R. Podeszwa).

0009-2614/$ - see front matter Ó 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2007.07.065

for the electrostatic interactions, but similar ideas involving separation of the exact and asymptotic components can be applied to the first-order exchange energies and to the second-order terms. If the electron densities of monomers are known, the exact (for a given density) electrostatic energy can be calculated from a simple formula involving two-electron Coulomb integrals. Such energy includes all the overlap (charge penetration) effects [1,2,5]. For large separations R between centers-of-mass (COMs) of the monomers, i.e., when R is several times larger than the size of either monomer, the overlap effects can be neglected and the electrostatic energy can be computed from an asymptotic expansion which can be expressed in terms of COM-located multipole moments and R [1,5]. For larger monomers, one usually applies the so-called distributed expansions, utilizing multipoles located at several centers on each monomer, typically all atoms. Several methods exist for calculating multipole moments assigned to each center [5–12]. One of the most popular approaches is the distributed multipole analysis (DMA) method developed by Stone and collaborators [5–9]. There exist also methods that combine the exact and asymptotic formulas [13–18]. The latter methods are somewhat similar to the one pro-

316

F. Rob et al. / Chemical Physics Letters 445 (2007) 315–320

posed in the present Letter and will be discussed in more detail later on. 2. Localized approach to theory of intermolecular interactions When two molecules are very far apart, not only the electrostatic energy, but the whole interaction energy is very well represented by the asymptotic expansion in inverse powers of R. This expansion, determined by monomers’ multipole moments and polarizabilities, is never exact, but can be arbitrarily close to the true interaction energy provided that R is large enough [5,1]. The asymptotic expansion neglects two effects: the physical interactions resulting from tunneling (exchanges) of electrons between interacting systems and the contributions dependent on the overlap of wave functions of the interacting systems in the electrostatic, induction, and dispersion energies. For large molecules, similarly as in the case of electrostatic energies, asymptotic expansions with distributed monomer properties have to be used [9]. However, these expansions are still not applicable at the most relevant conformations of van der Waals minima which typically ˚. involve closest approach interatomic distances of 2–3 A Thus, despite the huge effort invested in the development of distributed expansions [9], the practical usefulness of this approach is limited due to the neglect of the overlap and exchange effects. The approach proposed here may change this picture. The basic idea is as follows. For approximately spherical large monomers near their van der Waals minimum configuration in a dimer, only a relatively small number of atoms ˚ , i.e., within the range are separated by less than, say, 6 A where the overlap effects can be of appreciable importance. Therefore, the majority of atom–atom interactions can be treated accurately by the asymptotic expansion and only a small fraction has to be computed with the overlap effects. We can then distinguish the ‘near-range’ (NR) and ‘far-range’ (FR) fragments of each molecule, as illustrated in Fig. 1. Only the NR–NR interactions, i.e., the interaction involving the atoms (or more precisely compo-

Fig. 1. Acetaminophen dimer.

nents of wave functions) located in the shaded regions, need to be treated by SAPT(DFT), all the remaining ones can be computed from the asymptotic expansion. Due to this feature, the method will be denoted as LSAPT(DFT), where L stands for ‘localized’. Although intermolecular interactions cannot be reduced to atom–atom interactions, it is possible to separate them into interactions between regions of molecules by fitting the electron densities (related to multipole moments) and the linear response functions (related to polarizabilities) with properly chosen atom-centered basis sets [9]. The effectiveness of the method will depend on the shapes of monomers and their relative orientations. The best case scenario for LSAPT(DFT) is the case of two polymers in a collinear orientation, where the time of calculations of the interaction energies (not counting the DFT calculations for monomers) will be practically independent of the monomer size (beyond some minimal size). The most difficult case will be calculations for two planar molecules in the ‘sandwich’ configuration. In this case, one has to specify several NR/FR regions (domains) in each monomer. Even in this case, for sufficiently large monomers the calculations should scale linearly. To our knowledge, the idea outlined above has not been explored in theory of intermolecular interactions except for the electrostatic interactions. In this case, one can naturally approximate the electron density of a monomer by a sum of atom-centered densities. Fitting of densities [19] has recently been very popular in electronic structure methods. Although the standard density fitting techniques have not yet been applied in calculations of electrostatic energies, some related methods have been proposed. Wheatley and Mitchell [13] used an approach named by them ‘Gaussian multipoles’ (GMUL) method which is similar to DMA in that each product of primitive Gaussian functions appearing in the electron density is moved to the nearest nucleus using a truncated Taylor series. The total electron density assigned in this way to each atom can be used to compute atom-centered multipole moments. Depending on the distances between atoms, the electrostatic interactions can be computed using the atom-centered densities or multipoles. Another method of this type has been developed by Coppens et al. [14–18]. This approach, named the exact potential with multipole moments (EP/MM) method, also represents electron densities as a sum of atom-centered densities. The latter densities are universal for a class of molecules and are collected into a ‘databank of transferable aspherical pseudoatoms’. The electrostatic energies are then computed similarly as in the GMUL method. The method developed in Section 3 differs from GMUL and EP/MM mainly by the fact that we do not approximate the densities in the NR region. Ideas somewhat similar to those of Refs. [13–18] and of the present work are also used in the continuous fast multipole method (CFMM) [20]. However, to our knowledge, CFMM has never been applied to calculations of electrostatic intermolecular interaction energies.

F. Rob et al. / Chemical Physics Letters 445 (2007) 315–320

317

3. Localized electrostatic energy

where

The energy of electrostatic interaction between two unperturbed monomers with electron charge distributions qA and qB can be defined as Z Z XZ qA ðr1 ÞqB ðr2 Þ 3 3 ð1Þ Eelst ¼ qA ðrÞ d r1 d r2  jr1  r2 j b XZ Zb Za 3 qB ðrÞ  d r d3 r jr  Rb j jr  Ra j a X X ZaZb ð1Þ þ jRa  Rb j a b

~vaabb ¼

where Za and Zb are the nuclear charges, Ra and Rb are the positions of the nuclei, and the indices a and b run over the NA and NB atoms in monomer A and B, respectively. Eq. (1) can be rewritten in two equivalent forms. The first one uses a generalized potential Z Z ð1Þ Eelst ¼ qA ðr1 Þ~vðr1 ; r2 ÞqB ðr2 Þd3 r1 d3 r2 ; ð2Þ

~N ~FB Þ þ Eelst ð~ ~FB Þ þ Eelst ð~ ~N Eelst  Eelst ð~ qN qFA q qN qFA q Aq B Þ þ E elst ð~ Aq B Þ;

0 0

~vðr1 ; r2 Þ ¼ r1 12  þ

X Z a =N el X Z b =N el B A  jr1  Rb j jr2  Ra j a b

X Z a Z b =ðN el N el Þ A B ; jR  R a bj a;b

ð3Þ

where r12 = jr1  r2j and N elA and N elB are the numbers of electrons in monomers A and B, respectively. The other form uses generalized densities: Z Z ð1Þ ~A ðr1 Þr1 ~B ðr2 Þd3 r1 d3 r2 ; q ð4Þ Eelst ¼ 12 q where ~A ðrÞ ¼ qA ðrÞ þ q

X

Z a dðr  Ra Þ;

ð5Þ

a

and similarly for the monomer B. The Dirac function d(r  R) represents the unit point charge of the nuclei. If qX is obtained from a Hartree–Fock (HF) or a Kohn– Sham (KS) determinant, it has the form X

qX ðrÞ ¼

N bas X

P Xll0 vXl ðrÞvXl0 ðrÞ

ð6Þ

ll0

where vXl are the basis functions (atomic orbitals) for monomer X, X = A or B, N Xbas is the number of these functions, and P Xll0 are the density matrix elements defined for closed-shell systems by N occ molecular orbital expansion X coefficients cXli occ

P Xll0

¼2

NX X

cXli cXl0 i :

ð7Þ

i

The electrostatic energy can now be written as A

ð1Þ Eelst

¼

B

N bas N bas X X aa0

bb0

0 0

B ab PA vab ; aa0 P bb0 ~

ð8Þ

Z Z

A vA vðr1 ; r2 ÞvBb ðr2 ÞvBb0 ðr2 Þd3 r1 d3 r2 : a ðr1 Þva0 ðr1 Þ~

ð9Þ

2

B Eq. (9) requires calculating 14 ðN A bas N bas Þ two-electron integrals over the basis functions (unless prescreening is used [21]). This scaling makes such calculations time consuming and not possible for large monomers (consisting of about 100 or more atoms at the current levels of computer resources). To develop a fast and accurate method for calculations of the electrostatic energy, let us express this energy (exactly) as a sum of contributions ð1Þ

ð1Þ

ð1Þ

ð1Þ

ð1Þ

ð10Þ 0 ð1Þ ~R where Eelst ð~ qR B Þ are Aq N ~A þ q ~FA ties such that q

obtained from the generalized densi~A , and similarly for B. Any split¼q ting of the densities can formally be used, but we want the splittings to represent the densities in the NR and FR regions, as the labels suggest. It is possible to split the densities into two non-overlapping parts using a spacial division of each monomer as illustrated in Fig. 1, but if the splitting uses a basis set division, as it will be the case in our approach, the densities have to overlap. However, the major~N ~FX will be located in the NR and FR regions, ity of q X and q ð1Þ ~N respectively. Then, only the Eelst ð~ qN Aq B Þ contribution will have to be computed with the expensive approach of Eq. (8), the three other contributions from Eq. (10) can be accurately approximated with a multipole expansion. To define the splitting of density, let the indices of P be ordered in such a way that l, m = 1, . . ., NN correspond to orbitals centered in the NR region and l, m = NN + 1, . . ., NN + NF – in the FR region, where NN and NF are the numbers of basis functions centered in the NR and FR regions, respectively. The density matrix assumes then a block form ! PNN PFN ; ð11Þ PNF PFF

where PNN represent the submatrix with l 6 NN and m 6 NN, and so on [PFN = (PNF)T since the density matrix is symmetric]. Blocks PNN and PFF contain the majority of the density in the NR and FR regions, respectively, since any product of Gaussian functions centered in a certain region is also centered in the same region (except when the regions are non-convex but the exponential decay of the density matrix elements with distance [22] ensures that such contributions are small). In contrast, PFN contributes to densities in both region. Nevertheless, we assign the whole of PFN to the FR region. Although this choice may appear too crude, it works (as it will be shown below) due to the exponential decay mentioned above. The nuclear charges are also split into NR and FR regions. Thus, we can write explicitly the generalized density for the NR region as

318

F. Rob et al. / Chemical Physics Letters 445 (2007) 315–320

~N q A ðrÞ ¼ 

NN X

A A PA lm vl ðrÞvm ðrÞ þ

l;m¼1

X

Z a dðr  Ra Þ;

ð12Þ

a2NR

and similarly for monomer B. The generalized density in ~A  q ~N ~FA ¼ q the FR region is then given by q A . The NR– NR contribution to the electrostatic energy is computed ~N ~N using q A and q B in Eq. (4). These two densities are then used to compute the distributed multipoles of the DMA method [6]. Such moments are also computed separately ~FA and q ~FB . The remaining three terms in for the densities q Eq. (10) are then calculated using the multipole expansion. We used the GDMA program described in Ref. [8] to compute distributed multipoles (up to L = 10 except for L = 6 in 6-31G basis) with the traditional algorithm [7] (the parameter ‘switch’ in the program equal to 0). 4. Results To define the NR and FR regions, we assume a ‘cutoff’ distance. An atom of monomer A belongs to the NR region if the minimum of its distances to the atoms of monomer B is smaller than the cutoff distance. This simple definition is sufficient for topologies such as presented in Fig. 1, but not for cases like two parallel polymers, where several domains have to be defined. For testing our method, we have selected the acetaminophen (C8H9NO2) dimer shown in Fig. 1. The monomer geometry has been extracted from the crystal structure [23]. For the dimer configuration, we used a structure with ˚ , the second shortest R in the crystal. This R = 5.894 A dimer structure is centrosymmetric. In this configuration of the dimer, the overlap effects are fairly large. We performed calculation with three basis sets: 6-31G [24], cc-pVDZ [25], and aug-cc-pVDZ [26]. The results for different cutoff lengths are plotted in Fig. 2. The monomer density was obtained with the PBE0 DFT method [27,28] using the DALTON [29] program. The GRLB asymp-

10000

6—31G cc—pVDZ aug—cc—pVDZ mixed

1000

Error (%)

100 10 1 0.1 0.01 0.001 1e—04 3.5

4

4.5 5 5.5 Cut—off length (Å)

6

6.5

Fig. 2. Error in localized electrostatic energies for acetaminophen dimer with increasing cutoff length. For the description of the ‘mixed’ basis set, see text.

totic correction [30] was used with the experimental ionization potential of 0.278 hartree [31]. Fig. 2 shows that when one increases the size of the NR region, the error with respect to the exact result in a given basis set decreases very quickly for non-augmented basis sets and becomes practically negligible (below 1%) already ˚ , when only six atoms per monowith the threshold of 3.7 A mer are included in the NR region. On the other hand, for the augmented basis (aug-cc-pVDZ) the errors of the localized method are very large unless the cutoff is larger than ˚ (12 or more atoms per monomer in the NR region). 4.8 A Even for such cutoffs, the errors are several percent, reach˚ (this value of ing as much as 33% for the cutoff of 5.743 A the cutoff has to be chosen to increase the number of atoms per monomer in the NR region from 15 to 16). What we observe here is the well known difficulty of DMA appearing for basis sets with diffuse functions, see Ref. [8] for an extended discussion of this problem. If such functions are used, the distributed multipole moments are generally large and show significant oscillations with changes of the basis set [8]. Table 1 shows that the large errors of the electrostatic energy computed in the aug-cc-pVDZ basis are connected with the excessive charge splittings between the NR and FR regions. Similarly, we found that the higher multipole moments are large, distributed or not. In fact, the large errors of the localized approach are almost entirely due to the slow convergence (or even divergence) of the asymptotic part of our calculations since the NR–NR interactions are calculated exactly [from Eq. (4)] despite the excessive charge splittings. The divergence is indicated by the result with the expansion terminated at L = 6 equal to 2.043 kcal/mol being more accurate than the L = 10 result ˚ . No good of 3.199 kcal/mol, both for the cutoff of 5.743 A solution was found for this problem so far. Stone has recently recommended [8] to abandon the function-spacebased distribution for diffuse orbitals and use instead physical-space-based partitions (which requires numerical integrations over the density). This approach would not help in LSAPT due to the excessive charges introduced already by the splitting of the density matrix. One might think that problems described above would be diminished if the density matrix were expressed in terms of primitive Gaussian functions. Then, the elements of PFN could be uniquely assigned to either NR or FR region since the product of two Gaussian functions is another Gaussian function centered between them (analogously to the first step of the DMA analysis [6]). However, we have checked that such a procedure leads to as excessive charge splittings as the present one. In fact, the problem is not related to the partition of the monomer. We have performed a regular DMA calculation for the whole acetaminophen molecule and the atomic charges were very large, with the largest one amounting to 22.4 electron charges. Such charges are much larger than found for the molecules studied in Ref. [8]. This may be due to the size of the acetaminophen or to the presence of the aromatic ring as the largest charges are on the aromatic carbons.

F. Rob et al. / Chemical Physics Letters 445 (2007) 315–320

319

Table 1 Electrostatic energies of the acetaminophen dimer in the aug-cc-pVDZ and mixed (aug-cc-pVDZ in NR region and cc-pVDZ in FR region) bases Cutoff

NRA

4.6 4.8 5.0 5.3 5.4 5.743 5.8

11 12 13 14 15 16 17

Aug-cc-pVDZ

Mixed

LSAPT

Error (%)

Charge

SAPT

LSAPT

Error (%)

Charge

13.029 2.511 2.526 2.497 2.459 3.199 2.386

444.74 5.00 5.61 4.39 2.80 33.75 0.26

15.07 9.98 9.75 9.22 9.02 18.91 0.57

2.286 2.301 2.302 2.358 2.359 2.386 2.388

2.345 2.306 2.300 2.358 2.359 2.387 2.388

2.56 0.21 0.09 0.02 0.02 0.04 0.01

0.032 1.390 1.246 0.940 0.620 0.848 1.225

The column NRA gives the number of atoms in the NR region of monomer A (same as B). The ‘charge’ column gives the magnitude of charge on the NR (same as on FR) region of monomer A (or B). The exact electrostatic energy in aug-cc-pVDZ is 2.392 kcal/mol.

As a way around the problem introduced by the use of diffuse functions, we propose to apply ‘mixed’ basis sets – augmented in the NR region and non-augmented in the FR region. This choice is reasonable since the diffuse functions are needed chiefly in the NR region to represent well the tails of density which are critical for an accurate description of overlap effects. At the same time, since the diffuse functions are located only in the NR region, the excessive charge splitting should be eliminated. This is indeed the case, as shown by the results in Table 1 obtained with the aug-cc-pVDZ basis in the NR region and the cc-pVDZ basis in the FR region. The charge splitting is now always small, never exceeding 1.4 electron charges. Accordingly, the localized method converges well, see Table 1 and Fig. 2. The errors fall below 1% for ˚ . The decrease of the basis set by the cutoff of about 4.8 A the omission of some diffuse functions leads, of course, to a deterioration of the exact electrostatic energy. For the ˚ , this energy in the mixed basis is cutoff of 4.8 A 2.30 kcal/mol compared to the result in aug-cc-pVDZ basis equal to 2.39 kcal/mol, a 3.8% difference, ˚ cutoff. The exact which decreases to 1.4% for the 5.3 A result in cc-pVDZ constitutes only 56% of the aug-ccpVDZ value. To find out how the size of the NR region for an assumed accuracy changes with the size of molecule, we calculated the electrostatic interaction energy between two a-tocopherol (C29H50O2, vitamin E) molecules. The monomer geometry is taken from Ref. [32]. The dimer was arranged in a head-to-head orientation, shown in ˚ and the distance between the closFig. 3, with R = 20.5 A ˚ . The calculation was perest-approach atoms of 2.226 A formed in the 6-31G basis set using PBE0 and the ionization potential of 0.265 hartree [33]. As one can see in Fig. 3, the errors in the localized approach decrease very fast. The size of the NR region needed to reach 1% accuracy is virtually the same as in the case of the acetaminophen dimer. Also the computational effort beyond the DFT calculations for the monomers is almost the same for both systems. Thus, in this (favorable) case the scaling of the method with the system size is better than linear.

Fig. 3. Errors in a-tocopherol calculations with 6-31G basis set.

5. Conclusions We have shown that the localized approach to calculations of intermolecular interaction energies works well for the electrostatic interactions in the acetaminophen dimer near its van der Waals minimum. With about a dozen of atoms included in the NR regions where the overlap effects are accounted for, the errors of the localized approach are around 1% compared to the exact values in the same basis set. To achieve a stable method, one has to use mixed basis sets, augmented only in the NR region, which leads to an additional error of the order of 1%. These errors should be compared to those given by the asymptotic expansions (for complete monomers). The DMA method gives for the same configuration the electrostatic energy of 0.709 kcal/mol, only 30% of the exact result. The COM asymptotic expansion diverges and gives a completely unphysical energy of +246 kcal/mol. The localized procedure works equally well for the much larger a-tocopherol dimer. Acknowledgement This work was supported by ARO DEPSCoR and NSF CHE-0555979 Grants.

320

F. Rob et al. / Chemical Physics Letters 445 (2007) 315–320

References [1] B. Jeziorski, R. Moszyn´ski, K. Szalewicz, Chem. Rev. 94 (1994) 1887. [2] K. Szalewicz, K. Patkowski, B. Jeziorski, in: D.J. Wales (Ed.), Intermolecular Forces and Clusters, Structure and Bonding, Springer, 2005, p. 43. [3] A.J. Misquitta, R. Podeszwa, B. Jeziorski, K. Szalewicz, J. Chem. Phys. 123 (2005) 214103. [4] R. Podeszwa, R. Bukowski, K. Szalewicz, J. Chem. Theory Comput. 2 (2006) 400. [5] A.J. Stone, The Theory of Intermolecular Forces, Clarendon Press, Oxford, 1996. [6] A.J. Stone, Chem. Phys. Lett. 83 (1981) 233. [7] A.J. Stone, M. Alderton, Mol. Phys. 56 (1985) 1047. [8] A.J. Stone, J. Chem. Theory Comput. 1 (2005) 1128. [9] A.J. Stone, A.J. Misquitta, Int. Rev. Phys. Chem. 26 (2007) 193. [10] W.A. Sokalski, R.A. Poirier, Chem. Phys. Lett. 98 (1983) 86. [11] G. Jansen, C. Ha¨ttig, B.A. Hess, J.G. Angyan, Mol. Phys. 88 (1996) 69. [12] D.S. Kosov, P.L.A. Popelier, J. Phys. Chem. A 104 (2000) 7339. [13] R.J. Wheatley, J.B.O. Mitchell, J. Comp. Chem. 15 (1994) 1187. [14] A. Volkov, X. Li, T. Koritsanszky, P. Coppens, J. Phys. Chem. A 108 (2004) 4283. [15] A. Volkov, T. Koritsanszky, P. Coppens, Chem. Phys. Lett. 391 (2004) 170. [16] A. Volkov, P. Coppens, J. Comput. Chem. 921 (2004) 25. [17] A. Volkov, P. Coppens, H.F. King, J. Chem. Theory Comput. 2 (2006) 81.

[18] P.M. Dominiak, A. Volkov, X. Li, M. Messerschmidt, P. Coppens, J. Chem. Theory Comput. 3 (2007) 232. [19] B.I. Dunlap, J.W.D. Connolly, J.R. Sabin, J. Chem. Phys. 71 (1979) 4993. [20] C.A. While, B.G. Johnson, P.M.W. Gill, M. Head-Gordon, Chem. Phys. Lett. 230 (1994) 8. [21] D.S. Lambrecht, C. Ochsenfeld, J. Chem. Phys. 123 (2005) 184101. [22] C. Ochsenfeld, M. Head-Gordon, Chem. Phys. Lett. 270 (1997) 399. [23] D.Y. Naumov, M.A. Vasilchenko, J.A.K. Howard, Acta Crystallogr. C 54 (1998) 653. [24] W.J. Hehre, R. Ditchfield, J.A. Pople, J. Chem. Phys. 56 (1972) 2257. [25] J.T.H. Dunning, J. Chem. Phys. 90 (1989) 1007. [26] R.A. Kendall, T.H. Dunning, R.J. Harrison, J. Chem. Phys. 96 (1992) 6796. [27] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [28] C. Adamo, M. Cossi, V. Barone, J. Mol. Struct. (Theochem) 493 (1999) 145. [29] DALTON, a molecular electronic structure program, Release 2.0, 2005, see http://www.kjemi.uio.no/software/dalton/dalton.html. [30] M. Gru¨ning, O.V. Gritsenko, S.J.A. van Gisbergen, E.J. Baerends, J. Chem. Phys. 114 (2001) 652. [31] S.G. Lias, Ionization Energy Evaluation, in NIST Chemistry WebBook, NIST Standard Reference Database Number 69, http:// webbook.nist.gov. [32] DrugBank, a comprehensive resource for in silico drug discovery and exploration, see http://redpoll.pharmacy.ualberta.ca/drugbank. [33] S. Nagaoka, K. Mukai, T. Itoh, S. Katsumata, J. Phys. Chem. 96 (1992) 8184.