Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://beams.chem.msu.ru/sotrudniki/larin/IJQC83.pdf Äàòà èçìåíåíèÿ: Wed Jan 5 02:00:28 2000 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 20:29:46 2012 Êîäèðîâêà:

Cumulative Coordinates for Approximations of High- Order Atomic Multipole Moments in Aluminosilicate and Aluminophosphate Sieves
A. V. LARIN,1 D. P. VERCAUTEREN
1

2

Laboratory of Molecular Beams, Department of Chemistry, Moscow State University, Leninskie Gory, Moscow, B-234, 119899, Russia 2 Institute for Studies in Interface Sciences, Laboratoire de Physico-Chimie Informatique, FacultÈs Universitaires Notre Dame de la Paix, Rue de Bruxelles 61, B-5000, Namur, Belgium Received 14 April 2000; revised 14 February 2001; accepted 19 February 2001

ABSTRACT: A scheme to obtain approximate analytical functions for the atomic
distributed multipole moments of the crystallographically different atoms within aluminosilicate and aluminophosphate sieves is discussed. Respective atomic multipole moments are derived within the CRYSTAL95 ab initio periodic Hartree­Fock code with different basis sets, from minimal STO-3G to 6-21G . In order to illustrate the possible applications, distributed analyses are carried out for various structural models from all-siliceous zeolites and aluminophosphates with ratio Al/P = 1 to hydrogen forms of aluminosilicates. Simple approximate forms based on charge and geometry coordinates are proposed for the high-order moments of each atom, which are further required for the calculation of the electrostatic field within the structures. The possibility to use this analytical approach to evaluate the electrostatic field within embedded cluster models is also shortly discussed. c 2001 John Wiley & Sons, Inc. Int J Quantum Chem 83: 70­85, 2001

Key words: periodic Hartree­Fock; distributed multipole analysis (DMA);
electrostatic field; zeolites; ALPO sieves

Dedicated to the fruitful carrier of Dr. Evelyne Cohen de Lara. Correspondence to: D. P. Vercauteren; e-mail: Daniel. Vercauteren@fundp.ac.be. Contract grant sponsor: Interuniversity Research Program on Reduced Dimensionality Systems. Contract grant number: PAI/IUAP 4/10.



International Journal of Quantum Chemistry, Vol. 83, 70­85 (2001) c 2001 John Wiley & Sons, Inc.


HIGH-ORDER ATOMIC MULTIPOLE MOMENTS

Introduction
he evaluation of the electrostatic interactions remains one of the most important tasks within theoretical studies concerned with the interaction between biomolecules and solvent particles. The importance of the electrostatic effects between polar components of a cluster was confirmed by the good agreement between the theoretically predicted and experimentally recorded geometries, for example, for hydrogen bonding components [1], as well as for other types of molecular clusters [2]. Also, large values of the measured electrostatic field in zeolites [3, 4] suggest that the interactions with guest molecules in their adsorbed state are essentially of electrostatic character. The computation of the electrostatic interactions can a priori be achieved via any distributed multipole analysis (DMA) scheme [5 ­ 9] wherein atomic and bond centers can be considered as sites for the decomposition of the electron density. To our knowledge, the idea of distributed moments was first developed on the basis of the partition of the electron density, evaluated at the iterative extended HÝckel theory level [5] and CNDO level [5, 6], including the atomic centers only. Presently, it is, however, known that an approach taking into account both the atomic and bond centers, as developed in Stone's work [7], is usually preferable as it is more precise in providing convergence of the electrostatic interactions in a larger available space around the molecule under study [10]. Reaching the same precision for the electrostatic field values with atomic centers only would indeed necessitate to calculate high-order moments. For example, it was shown that the knowledge of the multipole moments (MM) up to the fourth order is required to obtain field values with a precision below 1% with the Saunders et al. scheme using atomic centers [9]. In this sense, it has also already been clearly shown that, as soon as the MMs depend on the geometry of the atom studied [7, 9 ­ 11], the knowledge of the dependences of the MMs with respect to the atomic internal coordinates is preferred to any known MM values only. The possibility to derive approximations of the MMs with respect to any structural parameter is indeed attractive because it allows the calculation of the respective moments from the geometrical data of the model only, such strategy being particularly useful in solid-state calculations. For zeolite crystals, for example, ab initio

T

methods indeed often fail either owing to the large number of atoms per elementary unit cell or to the unordered location of the substituting atoms in the framework, both preventing the wide application of the embedded cluster techniques to these materials. Appropriate analyses of the MM behavior with respect to the geometry could thus help, for example, to evaluate the long-range contributions in most of the zeolite catalysts of industrial interest. To be complete it should, however, be stressed that the derivation of such approximations within these materials requires some care, mainly because the very wide interval for the interatomic distances and respective conformational angles within most of the zeolites considered could correspond to different favored electronic states. Hence, few useful relations between MMs and structural parameters were proposed so far. Beginning of the 1980s, a linear correlation was observed between the quadrupole moment Q0 (C) at 2 the carbon atom, calculated with Stone's scheme [7], and the occupation of orbitals and number of nitrogen atoms in the rings of a series of azabenzenes [11]. Later on, it was demonstrated that the atomic MMs within some flexible molecules (ethanol, glycine, ...) vary strongly with the conformational angle by expansion of the MM as a series of the different angles [12, 13]. In these studies, it was shown that the MMs change in an appreciable manner even when the atom is remote from the two atoms positioned on the rotational axis. Evidently, such a behavior does not allow to obtain transferable MM values, which would be useful for any DMA estimations within a relevant class of large flexible molecules. The transferability of the lowest MM, i.e., the Mulliken atomic charges, on the contrary, has been clearly verified for a series of all-siliceous zeolites, which generally present a narrow interval for the Si­O bond length and a wide Si­O­Si angle variation [14 ­ 16]. Simple one- and twodimensional approximate analytical forms were indeed proposed for the Si and O atomic charges with respect to the internal geometry of the framework atoms. The Si and O charge values of the silicalite/STO-3G and all-siliceous mordenite/621G obtained via the respective analytical dependences [14 ­ 16] coincide with those computed by periodic Hartree­Fock (PHF) calculations [17, 18]. Similarly, through other low-order moments, the O dipole moments were approximated with respect to the bond distances and angles for a series

INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY

71


LARIN AND VERCAUTEREN of all-siliceous zeolites, aluminosilicates, and aluminophosphates [19, 20]. A simple sine function was found to be appropriate for a qualitative representation of the dipole dependences for the O atoms through all the systems at different basis set levels. Together with approximations of the Mulliken charges of the Si and O atoms already shown [14 ­ 16], these dependences of the MMs could thus provide an easy way to calculate the electrostatic field within any arbitrary zeolite framework even if it would be untractable through direct PHF computations. In this work, we will analyse the high-order MMs of the various framework atoms of several threedimensional solids such as zeolite models and other crystalline sieves. In the next section, we discuss some aspects of the MM computations and present the models chosen. In the third section, a general expression for the approximate representations of the MMs of any arbitrary order is given. It is based on Stone's relation [7] between the central moment of a molecule and the atomic MM of lower or equal order located on the atoms and on the bonds (although the latter are not considered in the scheme implemented in CRYSTAL95 [21]). In the fourth section, we discuss the results of the use of the approximations with respect to the internal coordinates for the Si and O atoms within all-siliceous zeolites, for the Al, P, and O within the ALPO sieves, and for the H, Al, Si, and O atoms within the H-form aluminosilicates. Deviations between the calculated and approximated moment values are analyzed in the last section. TABLE I
Symbol and number of atoms of different Si and O types, of atomic orbitals (AO) per unit cell (UC), and symmetry group of the all-siliceous zeolites [22]. Symbola AO/UC (6-21G ) 138 138 184 368 460 552 552 552 Symmetry group P32 21 P62 22 P42 21 2 I41 /amd P21 21 2 P43n C2/m R3c

Computational Aspects
The all-siliceous zeolite models considered herein were already studied earlier with various different basis set levels, from minimal STO-3G to 6-21G, ps-21G , and 6-21G (this last one, however, just only for the chabazite, CHA, and montesommaite, MON, frameworks) [14 ­ 16]. Herein, we included computations for the low-pressure form of cristoballite (CRI), - and -quartz, edingtonite (EDI), dachiardite (DAC), and sodalite (SOD), studied at the 6-21G level (Table I). This allowed to analyze the behavior of different order multipole moments (MMs) of Si and O versus their internal geometry in the respective frameworks. Four aluminophosphates (ALPOs), ATN, AST, ATO, CHA, were also considered at the 6-21G level to include the MMs of Al and P, and of other O atoms in a second type of framework (Table II). To test the influence of the geometry, nonoptimized and optimized sets of structural models were taken into account. The total geometry optimizations were performed for all-siliceous zeolites and ALPOs (with the exception of the AlPO4 -34 or CHA framework) using molecular mechanics based on the empirical Burchart­Kramer­van Santen (BKS) force field [25] as implemented in Cerius 2 [26]. To our knowledge, the only optimization reported so far is the one of the CHA aluminophosphate structure, using a plane wave method [24]. Five "small-size" hydrogen form zeolites were optimized at the STO-3G level [20] with the CRYSTAL95 version, in which we adopted the Polak­

Name

Atoms/UC 9 9 12 24 30 36 36 36

nSi /nO 1/1 1/1 1/1 1/3 3/5 2/1 4/9 1/4

-Quartz -Quartz
Cristoballite Montesommaite Edingtonite Sodalite Dachiardite Chabazite
a b

-QUAb -QUAb
CRIb MON EDI SOD DAC CHA

Ref. [23]. Own notation for convenience.

72

VOL. 83, NO. 2


HIGH-ORDER ATOMIC MULTIPOLE MOMENTS TABLE II
Symbol and number of atoms, of different Al, P (nP = nAl ), and O types of atomic orbitals (AO) per unit cell (UC), and symmetry group of the aluminophosphate sieves [22], all of them corresponding to the Al/P = 1. Name MAPO-39 AlPO4 -16 AlPO4 -31 AlPO4 -34b
a b

Symbola ATN AST ATO CHA

Atoms/UC 24 30 36 36

nAl /nO 1/4 2/3 1/4 1/4

AO/UC (6-21G ) 368 460 552 552

Symmetry group I4 F23 R3 R3

Ref. [23]. Ref. [24].

Ribiere algorithm [27]; then, single-point calculations were considered at the ps-21G level. Four of them were further successfully considered with the 6-21G one (Table III), with the exception of the HCHA framework which requires a strong variation of the Gaussian exponents at the 6-21G level. All MMs related to the atomic positions up to fourth order were determined via the distributed multipole analysis scheme proposed by Saunders et al. [9] available in the CRYSTAL95 [21] ab initio Hartree­Fock code for periodic systems. The Durant­Barthelat pseudopotential ps-21G basis set was considered with sp/d exponents as recommended in Ref. [21]. The used Gaussian sp/d exponents for the 6-21G basis are 0.9, 0.15/0.35, 0.15/0.4, and 0.42/0.72 a.u.-2 on the H, Al, Si, and O atoms within the H forms, 0.17/0.45 and 0.36/0.65 a.u.-2 on the Si and O atoms within the all-siliceous zeolites, and 0.14/0.35, 0.15/0.5, and

0.45/0.72 a.u.-2 on the Al, P, and O atoms within the ALPOs. All computations with the CRYSTAL95 code were carried out on an IBM 15-node (120-MHz) Scalable POWERparallel platform (with 1 Gb of memory/CPU). For all cases, the thresholds for the calculations were fixed to 10-5 for the overlap Coulomb, the penetration Coulomb, and overlap exchange, to 10-6 and 10-11 for the pseudo-overlap exchange, and to 10-5 for the pseudopotential series. A typical total geometry optimization of one bridged BrÆnsted center (12 variables) within STO3G takes 2­3 days on the above cited CPU. The single-point computations with the split-valence basis were executed directly without keeping the bielectronic integral files. The respective shortest convergence (7­8 cycles) took around 1.5­2 h in the case of the BrÆnsted centers of the HABW zeolite and around 20 min for any form of cristoballite with 6-21G .

TABLE III
Symbol and number of atoms, of different Al, Si, and O types (nH = nAl ) of atomic orbitals (AO) per unit cell (UC), and symmetry group of the H-form aluminosilicates [22]. Symbola HABW HCAN HCHA HEDI HNAT AO/UC (6-21G ) 388 582 517b 480 480 Symmetry group Pna21 P63 R3 P21 21 2 Fdd2

Name ABW Cancrinite Chabazite Edingtonite Natrolite
a b

Atoms/UC 28 42 39 34 34

nAl /nSi /nO 1/1/4 1/1/4 1/3/8 1/2/5 1/2/5

Ref. [23]. With ps-21G basis set.

INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY

73


LARIN AND VERCAUTEREN

Approximate Expressions for the Multipole Moments
A multipole integral expression for any arbitrary multipole moment (MM) of L-order (equation II.4.15 in Ref. [28]) has been determined for separate shells, defined by all quantum numbers of the atomic orbitals. It includes the summation over all electronic orbitals i (r)bothcentered onthe nuclei within the elementary cell 0 as well as shifted by the direct translational lattice vector g belonging to the shells of the respective atoms. For example, for a given atom A with coordinates (XA , YA , ZA ), after summation over all shells of the atom, the m component of the L order atomic MM may be written as: Qm (A) = L
1 2 g mm â NL XL (r,A) dr, (1)

absolute dipole moment values for the different groups of O atoms, i.e., Si­O­Si, Si­O­Al, and Si­ O(H)­Al. Following this simple idea of dependence between the moments on the different centers, we adopted an expression connecting the central and local moments similar to equation (11) in Ref. [7] developed by Stone for any atom or any given site different from an atom in the molecule:
N L S

Q (A) =
m L i = 1 S = 0 P =-S

L+m S+P

L-m S-P

-P â QP (i)Rm-S (A, i), (3) S L

P

g 12

0 1 (r)2 (r)

g

where Rm (A, i) corresponds to the respective LegenL dre polynomial whose argument is the vector between the considered atom A and site i, Qm (i) is the L m component of the L-order MM, the summation i running over all the N neighbor atoms of A. Even if only the very first term of the series in Eq. (3) is considered:
N L S

where P is the density matrix, coefficients N = (2 - m0 ) â (L -|m|)!/(L + |m|)!, m0 being the delta m function, and the real solid harmonics XL (r,A) are expressed as:
m L m XL (r,A) =

g 12

Qm (A) = L
i = 1 S = 0 P =-S N

a

LmSP

-P QP (i)Rm-S (A, i) S L

=
i=1

a

Lm00

Q0 (i)Rm (A, i) + ··· , 0 L

(4 )

DL (t, u, v)(x - XA )t (y - YA )u (Z - ZA )
(t,u,v)

v

(2)

the coefficients DL (t, u, v) being tabulated in Appendix A of Ref. [28], and r = (x, y, z) relates to the electron position. In previous studies, we derived approximations for the Mulliken atomic charges, i.e., the Q0 (A) 0 moments, for the A = Si and O atoms within all-siliceous zeolites [14 ­ 16], and for A = Al, P, and O within aluminophosphate sieves with ratio Al/P = 1 [19], based on the assumption that only the local geometry with respect to the nearest neighbors would be important for the resulting charge distribution, i.e., expressed through the moment on the atom studied. This idea was later developed for the local dipole moment Qm (O) of all types of crys1 tallographically different O atoms within series of aluminophosphates [19], as well as for the corresponding H forms as considered herein, plus the CHA framework at the ps-21G level [20]. More particularly, it was proposed that the local atomic moment should be proportional to the analogous "central" MM of the respective zeolite fragment. Simple geometrical relations between the local and central dipoles were shown to be precise enough for leading to satisfactory approximations of the

with aLmSP being dependent on the Rm (A, i) funcL tions used, a useful relation can be obtained for the behavior of the atomic MM with respect to the charge and geometry of the respective fragment. The term in the right-hand side of the series expansion (4) then corresponds to the nuclear contribution of the central L moment of the fragment including N neighbor atoms and centered on atom A. As an illustration in the case of a molecule, the relation between the total central moment and its "nuclear " part is nearly linear, for example, for dihydrogen in the gas state [29]. But going now to three-dimensional (3D) solids, it is evident that the choice of the closest neighbors within the crystal will not be limited by the formulas presented above and that the definition of the number of neighbors N requires testing. Unfortunately, the Rm (A, i) polynomials preL sented in Stone's work [7] distinguish from those applied by Saunders et al. to construct the MMs [9, 28]. In order to apply relation (3) to the MM calculated with CRYSTAL95, one should replace the normalized polynomials of Stone's work by the unm normalized functions XL (r, A) [Eq. (2)] considered in CRYSTAL. Then relation (3) provides a precise connection between the central and all local moments of low order of the fragment. Such adap-

74

VOL. 83, NO. 2


HIGH-ORDER ATOMIC MULTIPOLE MOMENTS tation is presently under progress in our group. Meanwhile, here we will consider the results of the approximations with the first term only, i.e., via formula (4). Replacing the angular Legendre polynomials Rm (A, i) by those expressed through the L nucleus Cartesian coordinates of atom A and its neighbor atom i (Xi - XA , ...) as developed in [28], one gets
m XL (A, i) = (t,u,v)

Dm (t, u, v)(Xi - XA ) L

t

â (Yi - YA )u (Zi - ZA )

v

(5)

in the same form as in expression (2). Then, one can deduce the coordinates for the charge and geometry dependences of the MMs from Eq. (4) as: Qm (A) = a L
Lm00

Rm (A) + bLm00 , L

(6 )

to the first neighbor ones. For most of the atoms within the structural models herein studied, i.e., zeolites and other aluminophosphates, an evident first choice for the neighbor atoms to be included in Eq. (6) can be proposed easily. For Si, Al, and P in all the frameworks, it is logical to consider four neighbors; for both H and O in the Si­O(H)­Al moiety, one considers three neighbors, and for O in the T­ O­T moieties, one considers the two first T atoms as neighbors. For all crystallographically different atoms A, we thus will compute the charge Q0 (A) 0 values with the Mulliken approach in CRYSTAL95 and include them in the calculations of the Qm (A) L via expressions (8) and (6) without any approximate formula of the MM versus the internal geometric parameters, as given in Ref. [14 ­ 16, 19, 20].

where aLm00 and bLm00 will be fitting parameters, and the Rm (A) functions correspond to the unnormalL ized functions (5) considered in CRYSTAL:
N

Approximations of the Multipole Moments for Aluminosilicate and Aluminophospate Sieves
H-FORM ALUMINOSILICATES We first present approximations via Eq. (6) of the MMs for the H-form zeolites as they have been partly optimized with the CRYSTAL95 code [20] and because several tendencies are more clearly evident within this series of models. An instructive illustration of the importance of the neighbor MM values and of the geometry of the closest atoms only (not including the MMs) is presented for the dipole components of the crystallographically different Si and Al atoms within the H-form aluminosilicates. Dipole and quadrupole moments of these atoms appear owing to a distortion of the tetrahedra TO4 ; they are absent for a straight Td symmetry. Our first attempt was to approximate the dipole moments with function (6) [Figs. 1(a) and 1(c)] replacing the individual different Mulliken charge values computed by CRYSTAL for the neighbor O atoms within the TO4 by averaged equal charges. In Figures 1(a) and 1(c), no correlation can be observed between the MM Qm (T) L and the charge and geometry-dependent coordinate Rm (T). But when we include the charge differL ences into expression (8), there is a clear correlation between Qm (T) and Rm (T) [Figs. 1(b) and 1(d)], L L which proves that this inclusion is extremely important for the approximation of the dipole moment for both the T = Si and Al atoms. Interestingly, one also notes that the highest the MM, the lowest the influence of the approximate coordinate function (8)

Rm (A) = L
i=1

m Q0 (i)XL (A, i) 0

(7 )

function [used in Eq. (3)] of instead of the Stone's method, the summation i running over all the N neighbors of A. If coordinate form (7) from the expansion series (4) includes only the charge Q0 (i) 0 values, then one can show that the parameters aLm00 and bLm00 in Eq. (6) should be independent on the m value and used simply as aL and bL . This undermeans that the slope of the approximation for the different m components of L-order MM should be similar. However, in order to emphasize the difference between the central moment described by relation (3) and a local one, which should be less dependent on the moments of the remoted atoms versus the closest ones, we decided to consider, instead of coordinate form (7), a modified form:
N

Rm L

Rm (A) = L
i=1

m Q0 (i)XL (A, i) d 0

-K iA

,

(8 )

which includes a term inversely proportional to the distance between the i neighbor and A atom, diA = ((Xi - XA )2 + (Yi - YA )2 + (Zi - ZA )2 )1/2 , K being an empirical value whose choice should be discussed (see discussion below). At this stage, let us just say that all the results we will present further for all the systems have been obtained with K = 2L + 1. Both last functions Rm (A) [Eqs. (7) and (8)] should L thus be analyzed in terms of a quantity dependent on the closest N atoms, not necessarily restricted

INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY

75


LARIN AND VERCAUTEREN

FIGURE 1. Dependence of the dipole Qm (T) components versus the Rm (T) coordinate [Eq. (8) with K = 3] for the 1 1

(a, b) T = Al and (c, d) T = Si atoms, considering (a, c) equal and (b, d) different charges on the first four neighbor O atoms of all types of Al and Si atoms in the H-form aluminosilicates calculated with the 6-21G basis set: m = 0 (circles), |m| = 1 (squares), positive m values (open signs), negative m values (filled signs). Approximation [Eq. (6)] is depicted by solid line.

of the lower moments on the resulting fit of the higher moments, as shown in Figure 2. Hence, the influence on the T octupole moment approximation with equal or different charges in expression (8) becomes already very minor and is thus not shown herein. The better correlation (6) for the MMs of Al as compared to Si could be related with the difference of second-order neighbors for Si, which can be given in usual notation nSi(4 - n)Al (n = 0, 1, ... ,4). Actually, the limited number of zeolite models considered here does not allow to completely verify the validity of correlation (6) while expanding the number of neighbors up to 8 for the T atoms in Eq. (8). The basis set dependence remains a problem in any DMA scheme. The comparison between the Qm (T) components, T = Al, Si, L = 2, 3, obtained L with the ps-21G and 6-21G basis sets in Figure 2, illustrates the usual smaller slope of function (6) with ps-21G . Also, as seen from Figure 2(c) for the Qm (Si) with both basis sets, a quadratic function of 3 Qm (Si) versus Rm (Si) would be more appropriate 3 3 than just the linear form of Eq. (6). This deserves

a further discussion (see discussion below). The respective decrease of the Mulliken charges calculated using CRYSTAL varying the basis set from ps-21G to 6-21G is more important for the T atoms (whose core electrons are described by pseudopotentials with ps-21G ), i.e., almost a factor of two, as compared to the H atoms, which almost do not vary. One may hope, however, that the final choice of a correct basis set level could be reached when a complete set of theoretical values for a chosen simple structural model, for example, relevant parameters out of the infrared (IR) vibrational spectra of several adsorbed species, would meet agreement with the experimental data. Other higher moments of the O atoms and dipoles of the H atoms for which a correlation was rather satisfactory are shown in Figures 3(a)­3(c) and 3(d), respectively. As the data for the crystallographically different O atoms are too numerous and are presented for the three different types of T­O­T moieties, we have separated the components QL of L = 1, 3, and 4, for these three L moieties. One can see that the "splitting" between

76

VOL. 83, NO. 2


HIGH-ORDER ATOMIC MULTIPOLE MOMENTS

FIGURE 2. Dependence of the (a, b) quadrupole Qm (T) and (c, d) octupole Qm (T) components versus the Rm (T) 2 3 L

coordinate [Eq. (8)] for the (a, c) T = Si and (b, d) T = Al atoms, in the H-form aluminosilicates: m = 0 (circles), |m| = 1 (squares), |m| = 2 (triangles up), |m| = 3 (diamonds), positive m values (open signs), negative m values (filled signs). 6-21G basis (black signs); ps-21G basis (gray signs). Respective approximations [Eq. (6)] are depicted for ps-21G basis (dotted line) and 6-21G basis (solid line).

the three groups demonstrates clear differences between the groups. If most components QL (O), L = L 1, 3, and 4, behave nearly linearly [Eq. (6)] for all three Si­O(H)­Al, Si­O­Al, and Si­O­Si types, the dipole component Q-1 (O) for Si­O(H)­Al does 1 not correlate with R-1 (O). On the contrary, the 1 two last moieties, Si­O­Al and Si­O­Si, reveal very little variation between the Q-1 (O) dependences 1 with R-1 (O). Despite the different geometries of the 1 bridged BrÆnsted centers obtained by optimization with STO-3G [20], the respective H-dipole components Qm (H) for all crystallographically different 1 H atoms present very close slope values aLm00 in Eq. (6). This case [Fig. 3(d)] demonstrates the largest differences between the approximations (6) for different m components of the same L-order moment throughout all MMs considered in this work. In all other cases, approximations (6) nearly coincide, which is in accordance with our conclusion about the slope values aLm00 for different m components as drawn in the discussion on approximate expressions above. Although not presented here for all

different H atoms, we should also mention a satisfactory correlation for the MMs of L = 2, 3, and 4th order. The lower m = 1 components are usually better approximated than the upper ones for higher L > 1. In the case of the hydrogen atoms, we also tested the expansion of the N neighbors in expression (8) from 1(O) atom to 3(O, Si, Al) atoms; no improvement in the precision of the aLm00 and bLm00 parameters of function (6) was obtained. The insensitivity here observed for the MMs of the bridged hydrogens can be explained by a minor charge transfer between Si or Al and H as a result of a long T­H distance, T = Al, Si, which is in agreement with that determined experimentally from nuclear magnetic resonance (NMR) measurements [30]. However, such an influence of the number of neighbors should also be taken into account for the other atoms, in particular for the Si atoms presenting different second-order neighbors of the nSi(4-n)Al type.

INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY

77


LARIN AND VERCAUTEREN

FIGURE 3. Dependence of the (a­c) multipole moment Q-L (O) and (d) dipole Qm (H) components versus the Rm (X) 1 L L
coordinate [Eq. (8)] in the H-form aluminosilicates calculated with the 6-21G basis set. Notations for (a­c) cases: Si­O(H)­Al (filled circles), Si­O­Al (open circles), Si­O­Si (triangles up); approximations [Eq. (6)] are depicted for the Si­O(H)­Al (dotted line), Si­O­Al (solid line), Si­O­Si (dashed line). Notations for (d) case: m = 0 (circles), |m| = 1 (squares), positive m values (open signs), negative m values (filled signs); approximations [Eq. (6)] are depicted for m = 0 (solid line), m = 1 (dotted line), m = -1 (dashed line).

ALL-SILICEOUS ZEOLITES All-siliceous zeolites are less convenient systems for using approximation (6). A fairly good agreement was obtained only for two lowest allowed MMs, i.e., dipole moment of the O atom [Figs. 4(a)­ 4(b)] and for the octupole moment of Si [Fig. 4(d)], while just a qualitative correlation was observed for the hexadecapole moment of O [Fig. 4(c)]. Comparison of the dipole components of the oxygens calculated at the STO-3G or 6-21G levels illustrates the variation of the aLm00 value with the shift of basis set [Fig. 4(a)]. Further improvements of the basis set level should thus be checked on the basis of the closeness between other calculated values and experimental data. In this sense in another work [31], we have shown a good agreement between the experimental (for "large"-size zeolites) and calculated (for "small"-size zeolites) quadrupole coupling constants Cqcc of the 17 O, 2 H, and 27 Al atoms, which confirms the correct values of the gradient field on

these atoms. Further similar steps should be considered to reach a better coincidence with the measured electrostatic field too [3, 4]. One of the reasons for the deviations between the calculated MM values and approximation (6) could come from the use of nonoptimized structures. For example, the points near R-3 (O) 0 a.u. 4 in Figure 4(c) [filled black diamond, above linear dependence (6) for R-3 (O) for the MON structure, and 4 filled black triangle down, below dependence (6) for R-4 (O) for DAC], corresponding to an Si­O­Si an4 gle near 180 , increase the average deviation. The existence of zeolite structures with Si­O­Si angles near 180 was clearly questioned [32]. Very large changes within the framework geometry could lead to a respective shift in the electron density distribution. And evidently, all possible density variations cannot reasonably be presented with a single correlation function (6). We could thus think that such function would only be precise in a relatively narrow range of bond distances and angles. The ques-

78

VOL. 83, NO. 2


HIGH-ORDER ATOMIC MULTIPOLE MOMENTS

FIGURE 4. Dependence of the multipole moments Qm (X) components versus the Rm (X) coordinate [Eq. (8)] for the L L
(a­c) X = O and (d) X = Si atoms, in the all-siliceous zeolites calculated with the (a) STO-3G and (b­d) 6-21G basis sets: m = 0 (circles), |m| = 1 (squares), |m| = 2 (triangles up), |m| = 3 (diamonds), |m| = 4 (triangles down), positive m values (open signs), negative m values (filled signs). Nonoptimized models (black signs); BKS optimized models (gray signs). Approximation [Eq. (6)] is depicted by solid line.

tion then is: Could it be easily recommended for the systems under study? To our opinion, a pertinent choice of the type of optimization procedure could therefore also present a problem to tackle within the further studies. The good coincidence already shown between the geometries of some crystalline solids calculated by ab initio [33] and optimized using the GULP approach [34] suggests to try in a first step the purely empirical schemes. That is why we considered optimizations of the all-siliceous frameworks with the empirical BKS force field [25] as implemented in Cerius 2 [26]. The optimized values with the BKS force field, i.e., 163.1, 153.5, 147.1, and 158.7 for the Si­O­Si angle within the SOD, CHA, QUA, and -QUA, respectively, correlate well with the ab initio values obtained with CRYSTAL at the 6-21G level, i.e., 159.2, 151.9, 144.7, and 157 , respectively [33]. Moreover, a predictable conclusion is that the Si­O­Si angle near 180 observed within the DAC or MON frameworks cannot be simply altered just by optimization with the BKS force field.

The optimization "tries" to equalize all the Si­O distances to an average one. The Si­O­Si angles are usually larger and the average Si­O bond distances are smaller for all structures optimized with the BKS force field than the ones obtained by ab initio methods at different levels [33] (the same tendency holds for the results obtained at the 6-31G basis set level). It allows to hope that better force fields, which would consider Si­O­Si angles near 180 as less stable ones, would be developed in the future. Question of the MMs, the important result is that the optimization of the structure has almost no influence. The behavior of the MMs for the optimized geometry [gray signs in Figs. 4(b)­4(d)] is similar to the one for the nonoptimized initial models [black signs in Figs. 4(b)­4(d)]. This agrees well with the coincidence already noted between the charge [14] and dipole [19] values computed either for the optimized models, or for the ones on the basis of X-ray data which were given in [22]. And even more important is that all the structural changes due to optimization did not lead to a better agreement us-

INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY

79


LARIN AND VERCAUTEREN

FIGURE 5. Dependence of the (a) dipole Qm (O), (b) quadrupole Qm (O), and octupole Qm (T) components versus the 1 2 3

Rm (T) coordinate [Eq. (8)] for the (c) T = Al, and (d) T = P atoms, in the aluminophosphate sieves calculated with L 6-21G basis set: m = 0 (circles), |m| = 1 (squares), |m| = 2 (triangles up), |m| = 3 (diamonds), positive m values (open signs), negative m values (filled signs). Nonoptimized models (black signs); the BKS optimized models (gray signs). Approximation [Eq. (6)] is depicted by solid line.

ing approximation (6) neither for the dipole of the Si atoms, nor for their quadrupole moments. After the BKS optimization, the number of the points that increase the average deviation from approximation (6) becomes larger for the Qm (O) values [gray 1 signs in Fig. 4(b)] and lower for the Qm (Si) values 1 (not shown here). That is why one supposes that the role of an efficient geometry optimizations is not yet clear, but it rather could improve the resulting fitting of the MMs via Eq. (6). ALUMINOPHOSPHATE SIEVES In general, the approximation of the MMs of the O atom in the aluminophosphates (ALPOs) is better than the one for O within the all-siliceous models. The components of both the dipole [Fig. 5(a)] and quadrupole [Fig. 5(b)] moments of O can be satisfactorily presented by the respective approximations of linear type (6). The correlation coefficients for the lowest allowed octupole moment of the tetrahedral atoms Al and P, using approx-

imation (6), are 0.984 and 0.990, respectively. The question of the choice of an appropriate optimization procedure discussed above for the all-siliceous zeolites, however, remains relevant for the ALPOs. The same type of optimization of the AST, ATN, and ATO aluminophosphate sieves with the BKS force field does not show large improvements of the MMs for the Al, P, and O atoms for the optimized models (gray symbols in Fig. 5) as compared to those for the nonoptimized ones (black symbols in Fig. 5).

Discussion
Before discussing in more detail the reasons for some deflections from the suggested approximation (6), let us return first to the choice of the power K in the charge and geometry-dependent coordinate function (8). Surprisingly, this simple idea to introduce a term inversely proportional on the distance diA seems to be very useful to provide, for

80

VOL. 83, NO. 2


HIGH-ORDER ATOMIC MULTIPOLE MOMENTS

FIGURE 6. Dependence of the dipole Qm (O) components versus the Rm (O) coordinate [Eq. (8)] with (a) K = 0, 1 1

(b) K = 1, (c) K = 3, and (d) K = 4 in all-siliceous zeolites calculated with the 6-21G basis set: m = 0 (circles), |m| = 1 (squares), positive m values (open signs), negative m values (filled signs). Approximations [Eq. (6)] are depicted by solid line.

most cases, a smaller average deviation between the MM approximated via Eq. (6) and calculated directly by the CRYSTAL. Different behaviors with different K power values in coordinate form (8) can be illustrated by passing from K = 0 in Eq. (8) for the Qm (Al) and Qm (Si) components in the H-form 1 1 aluminosilicates, where no correlation is observed, to K = 3 [Figs. 1(b) and 1(d), where L = 1] which shows some correlation. The problem of the choice of a precise K value is not a simple task because it requires the best possible set of MM values for testing Eq. (8). It is thus complicated by the presence of the "wrong" points that are present in the chosen data, as for the Qm (O) within all-siliceous zeolites (Fig. 6), and 1 which increase the average deviation. Having no strict criterion for disregarding the unappropriate points, owing, for example, to the strongly distorted geometry, we estimated the K power value empirically. Looking to Figure 6, one notes that K = 2L + 1 generally provides a better precision in approximation (6) as compared to the other analyzed values up to K = 3L + 1. This arbitrary (empirical) introduction

of K in Eq. (8) can, however, also lead to distorted contributions to the central moment, obtained using a relation analogous to Eq. (4) with polynomials (5), from the local MMs of the neighbor atoms located within the different "shells" (first, second, ...). That is why we suggest that a "final" choice of the K value could be done analyzing the obtained MMs in relation with the analyses of the contributions of the second-order neighbors, i.e., N = 8 for the O atoms within the T­O­T moieties and all Al, Si, and P tetrahedral atoms, and N = 9 for the O atoms within the Al­O(H)­Si moieties. With the exception of a larger number of neighbors for H, we did not consider the contributions from the different shells here and thus have used K = 2L + 1 as estimated above. The linear approximation (6) of the high-order multipole moments (MMs) logically extends the earlier findings deduced for the Mulliken charges [14 ­ 16, 19, 20]. The charges cannot evidently be described in terms of the charge and geometrydependent coordinate (7) or (8). The right-hand side in Eq. (7) then becomes simply the total charge of

INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY

81


LARIN AND VERCAUTEREN all the N neighbor atoms for the Q0 (A) moment. 0 Despite the fact that a dependence clearly exists between the atomic charge of a central atom and those of the neighbor ones, we prefer to consider as a first step an approximate function based only on the geometrical coordinates in order to have the possibility to construct a consequent scheme for all MMs, i.e., L = 0 L = 1 ··· . Hence, all the dependences derived for the Si and O charges [14 ­ 16, 19, 20] are still important for any further development. So far, all moments cannot be obtained with a satisfactory precision but future studies with higher terms in expansion series (4) should complete the required MM approximations. We emphasize, however, in accordance with our proposition when discussing approximate expressions above, that most of the cases of different m components of L-order MM can be approximated with Eq. (6) with nearly the same slope value aL as can be seen from Figures 1­5. The validity of dependence (6) for all systems considered is also confirmed by the small free term value bL approaching zero, i.e., being usually of the same order value as the error of the bL estimation. The precision of the expansion series (4) is closely connected with the role of the atomic centers to effectively "represent" the total electronic density. On one hand, application of more sophisticated basis sets of higher level should lead to the appearence of important areas of enhanced electron density, which could need the necessary introduction of bond centers (or any additional points) for a correct representation of the density. Hence, the MMs calculated with a higher basis set quality could deviate from our approximation (6) limited by the atomic positions only. On the other hand, we should consider the basis set of the highest level for the most precise representation of the electron density. We would like to show herein that a compromise between these two tendencies can be reached for the materials studied. As we have shown above, the most important MM, i.e., the dipoles on O atoms and octupoles on the T atoms, can be well presented with the linear form at the 6-21G level. This basis set level should also be sufficient for the computation of the electrostatic field gradients on framework 17 O, 2 H, 27 Al atoms [31], as it has been shown that the close 6-31G basis set provides correct values of NMR shielding constants [35]. Two evident problems, however, appeared herein with the a priori most suitable minimal STO-3G basis set [its suitability is suggested, without showing it here, for example, by the nice linear behavior (6) of the O quadrupole Qm (O) and Si hexadecapole 2 Qm (Si) components]. The dipole and quadrupole 4 moments for the T atoms, T = Al and Si, within both the ALPOs and all-siliceous structures (including the BKS-optimized models) cannot be approximated with Eq. (6); only a satisfactory correlation was found for the Qm (T), T = Si and Al, L = 1­3, L within the H forms. One could suggest that the sophisticated behavior of the dipole and quadrupole moment components requires a more delicate geometry optimization than the one we have used so far. Another interesting relevant aspect is the "intrinsic" cumulativity of the Saunders et al. scheme [9]. This scheme, as also Stone's scheme [7], does not include any relation between the atomic MM of different orders on the same center, in the manner formulated by Sokalski et al. within the cumulative atomic multipole moments (CAMM) method [8]. On the contrary, dependence (6) between the MMs of different orders related to the chemically bonded atoms in zeolites and other sieves could be considered as analogous to the "cumulativity" concept. Such property means that with the hierarchical construction of the lower order moments as we suggest, we could succeed to approximate the higher order moments on the neighbor atoms. Approximation (6) demonstrates the capabilities of a scheme based on the atomic centers only. However, it has been shown that a decomposition including the positions of bond centers provides a more powerful convergence that can be crucial for an isolated molecule like N2 [10]. The bond positions depend on the wave functions applied (if they are not chosen arbitrarily), while the atomic positions can be determined experimentally and hence serve easily as decomposition sites. This last convenience might be sufficient to prefer the atomic positions only for the simulation of the electric field provided that transferable functions for the MMs versus the atomic geometry are estimated from another source, using, for example, expression (6) presented herein. Dependences (6) and (8) could clearly be useful in further modeling of solid structures for which most known approaches neglect the charge variation with geometry. Scheme (6) derived herein allows also to consider the respective variation of the higher moments. It has been shown that their contributions to the resulting electrostatic potential (EP), in the important zeolite spaces available for small adsorbed molecules, decrease the latter by a factor of two as compared to the values produced by the atomic charges only, both for aluminosilicates [20] and aluminophosphates [19]. It should

82

VOL. 83, NO. 2


HIGH-ORDER ATOMIC MULTIPOLE MOMENTS

FIGURE 7. (a) Electrostatic potential values EP (a.u.) with respect to the Al­O(1)­P plane of the ATN framework

calculated with the 6-21G basis sets; (b­d) EP differences presented as (1 - EP/EP(3)) â 100(%) obtained using (b, c) Mulliken charges on all atoms and high-order dipole and quadrupole moments on O and octupoles on Al and P and (d) Mulliken charges only relative to the potential representation EP(3) allowing all the moments up to third order on all atoms. EP differences include values of high-order MMs from (b) CRYSTAL calculation and (c) approximation via Eq. (6).

thus essentially improve the empirical simulation of adsorbed species. Therefore, we have thus verified the precision of approximation (6) for constructing an EP map. For this, we computed the EP within

the section passing through the three Al, P, and O(1) atoms of the MAPO-39 sieve with the ATN framework [Fig. 7(a)]. Previously, it was shown that using all MMs up to the octupoles on all atoms

INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY

83


LARIN AND VERCAUTEREN provided a very reliable EP map [19]. Here, the EP was built considering MM values for which we suggested the use of approximations, namely dipoles and quadrupoles on O and octupole moments on Al and P (Fig. 5). In Figures 7(b) and 7(c), we present difference (1 - EP/EP(3)) â 100% maps using Mulliken charges taken from the CRYSTAL calculation as well as dipoles and quadrupoles on O and octupole moments on Al and P. The difference between both figures is that in Figure 7(b) high-order moments were calculated directly with CRYSTAL, while in Figure 7(c) they were obtained from the respective approximation (6). In both cases, the dipole and quadrupole moments on Al and P were neglected. It is very instructive that disregarding these lower order MMs on Al and P allows nevertheless a correct presentation of the EP [Fig. 7(b)] in the space available for small adsorbed molecules (upper right corner). Very large differences (1 - EP/EP(3)) are usually observed along the line EP(3) = 0 [the line closest to zero corresponds to the EP of 0.03 a.u. value in Fig. 7(a)], but such a behavior is only a minor disadvantage of the EP differences. The differences illustrate more clearly the EP variation with different methods. Both maps provide better agreement with more precise EP calculation (when MM up to third order over all atoms are included) as compared to the Mulliken charges only [Fig. 7(d)]. Introducing the simple analytical dependences already obtained [14 ­ 16, 19, 20] for the charges versus the internal geometry in formula (6) or in more advanced relations based on expression (3) (equation (11) in Ref. [7]), one could develop a consequent approximation of all required MMs to reach a high precise electrostatic field representation. These calculations of the MMs and, hence, of the electrostatic field could be included in the more and more common embedded techniques wherein the electronic charge distribution of the most important atoms is usually calculated precisely via ab initio cluster calculation (corresponding to a "bielectronic" zone of the Saunders scheme [9]) and the long-range interactions ("monoelectronic" zone) would be approximated by the MM set. Evidently, a different partition between these zones as accepted in Ref. [9] and in the embedded approaches should be discussed more deeply. tems with the periodic Hartree­Fock CRYSTAL95 code using the 6-21G basis sets, i.e., "small-size" all-siliceous zeolites, aluminophosphates, and an equivalent series of hydrogen forms of aluminosilicates. Approximated functions for the MMs were first proposed in terms of a simplified expression [Eq. (8)] being a consequence of a general ratio between the central and local moments related to the different neighbor sites as developed by Stone [7] for an isolated molecule. Reasonable correlations were found for all the components of the dipole and quadrupole moments for the O atoms (not for the O quadrupole within the SiO2 systems) and for the octupole moments of H, Al, P, and Si atoms for all studied systems. It was shown that the inclusion of the differences between the neighbor O charges is important for constructing the approximate charge and geometry-dependent coordinate (8) for the high-order moments of the T = Si and Al atoms in the same series of H-form aluminosilicates. This influence, however, decreases with the order of the multipole moment. The larger number of neighbor atoms, i.e., O, Si, and Al, considered in the case of the H atoms does not lead to a better precision as compared to the approximation of the MMs of the H atoms regarding only one neighbor O. This insensitivity to the neighbor 's "expansion" within the derived approximation for H can be considered as a confirmation of the local character of the charge distribution on the bridged atoms. The possibility to improve the quality of the approximation of the electrostatic potential is illustrated for the aluminophosphate ATN framework. The proposed scheme allows to construct an hierarchical or cumulative (in terms of Ref. [8]) coordinate sequence based on the expansion series (4) including known geometry and approximations of the lower order MMs of the neighbor atoms to consequently describe the higher MMs. ACKNOWLEDGMENTS The authors wish to thank the FUNDP for the use of the Namur Scientific Computing Facility (SCF) Centre, a common project between the FNRS, IBMBelgium, and FUNDP as well as MSI for the use of their data in the framework of the "Catalysis and Sorption" consortium. We are grateful for the partial support of the Interuniversity Research Program on Reduced Dimensionality Systems (PAI/IUAP 4/10) initiated by the Belgian Government. All the authors thank Prof. R. Dovesi and Prof. W. J. Mortier for fruitful discussions.

Conclusions
Atomic multipole moments (MMs) are calculated for three types of three-dimensional model sys-

84

VOL. 83, NO. 2


HIGH-ORDER ATOMIC MULTIPOLE MOMENTS

References
1. Buckingham, A. D.; Fowler, P. W. J Chem Phys 1983, 79, 6426; Ibid. Can J Chem 1985, 63, 2018. 2. Price, S. L.; Stone, A. J. J Chem Phys 1987, 86, 2859. 3. Cohen de Lara, E.; Delaval, Y. J Chem Soc Faraday Trans II 1978, 74, 790. 4. Marra, G. L.; Fitch, A. N.; Zecchina, A.; Ricchiardi, G.; Salvalaggio, M.; Bordiga, S.; Lamberti, C. J Phys Chem B 1997, 101, 10653. 5. Rein, R. Adv Quantum Chem 1973, 7, 335. 6. Dovesi, R.; Pisani, C.; Ricca, F.; Roetti, C. J Chem Soc Faraday Trans II 1974, 70, 1381. 7. Stone, A. J. Chem Phys Lett 1981, 81, 233. 8. Sokalski, W. A.; Sawaryn, A. J Chem Phys 1987, 87, 526. 9. Saunders, V. R.; Freyria-Fava, C.; Dovesi, R.; Salasco, L.; Roetti, C. Mol Phys 1992, 77, 629. 10. Price, S. L. Mol Phys 1986, 58, 651. 11. Price, S. L.; Stone, A. J. Chem Phys Lett 1983, 98, 419. 12. Koch, U.; Popelier, P. L. A.; Stone, A. J. Chem Phys Lett 1995, 238, 256. 13. Koch, U.; Stone, A. J. J Chem Soc Faraday Trans II 1996, 92, 1702. 14. Larin, A. V.; Leherte, L.; Vercauteren, D. P. Chem Phys Lett 1998, 287, 169. 15. Larin, A. V.; Vercauteren, D. P. Int J Quantum Chem 1998, 70, 993. 16. Larin, A. V.; Vercauteren, D. P. Int J Inorg Mater 1999, 1, 201. 17. White, J. C.; Hess, A. C. J Phys Chem 1993, 97, 6398. 18. White, J. C.; Hess, A. C. J Phys Chem 1993, 97, 8703. 19. Larin, A. V.; Vercauteren, D. P. J Mol Catal A 2001, 106, 73.

20. Larin, A. V.; Vercauteren, D. P. J Mol Catal A 2001, 168, 123. 21. Dovesi, R.; Saunders, V. R.; Roetti, C.; CausÞ, M.; Harrison, N. M.; Orlando, R.; AprÞ, E. CRYSTAL95 1.0, User 's Manual, University of Torino, Italy, 1996. 22. (a) Newsam, J. M.; Treacy, M. M. J. ZeoFile--A Stack of Zeolite Structure Types (catalysis and sorption software and databases from Molecular Simulations Inc., San Diego CA, 1995); (b) Newsam, J. M.; Treacy, M. M. J. Zeolites 1993, 13, 183. 23. Barrer, R. M. Pure Appl Chem 1979, 51, 1091. 24. Jeanvoine, Y.; àngyÀn, J. G.; Kresse, G.; Hafner, J. J Phys Chem B 1998, 102, 5573. 25. van Beest, B. W. H.; Kramer, G. J.; van Santen, R. A. Phys Rev Lett 1990, 64, 1955. 26. Cerius 2, User 's Guide, Version 4.0.0; MSI: San Diego, 1997. 27. Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Wetterling, W. T. Numerical Recipes; Cambridge University Press: New York, 1988. 28. Pisani, C.; Dovesi, R.; Roetti, C. Hartree­Fock Ab Initio Treatment of Crystalline Systems; Springer: New York, 1988. 29. Poll, J. D.; Wolniewicz, L. J Chem Phys 1978, 68, 3053. 30. Fenzke, D.; Hunger, M.; Pfeifer, H. J Magn Reson 1991, 95, 477. 31. Larin, A. V.; Vercauteren, D. P. Int J Quantum Chem 2001, 82, 182. 32. Mortier, W. J., private communication. 33. Civalleri, B.; Zicovich-Wilson, C. M.; Ugliengo, P.; Saunders, V. R.; Dovesi, R. Chem Phys Lett 1998, 292, 394. 34. Gale, J. D. GULP (The General Utility Lattice Program); Royal Institution/Imperial Colledge: London, UK, 1992/1994. 35. McMichael Rohlfing, C.; Allen, L. C.; Ditchfield, R. Chem Phys 1981, 63, 185.

INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY

85