Convergence of electric field and electric field gradient versus atomic ...

Report 2 Downloads 82 Views
Convergence of Electric Field and Electric Field Gradient Versus Atomic Basis Sets in All-Siliceous and Mg Substituted Phillipsites A. V. LARIN, I. K. SAKODYNSKAYA, D. N. TRUBNIKOV

Department of Chemistry, MSU, Moscow, GSP-2, 119992, Russia Received 23 August 2007; Revised 8 February 2008; Accepted 10 February 2008 DOI 10.1002/jcc.20973 Published online 28 April 2008 in Wiley InterScience (www.interscience.wiley.com).

Abstract: Electrostatic potential (EP), electric field (EF), and electric field gradient (EFG) values are calculated in periodic models of magnesium substituted phillipsite (MgPHI) zeolite forms using periodic DFT (PDFT) hybrid B3LYP level with fourteen different basis sets. Relative root mean square differences between the EP, EF, or EFG values estimated between different basis sets are evaluated in several spatial domains available for adsorbate molecules in the zeolite. In these areas, the EF increase in MgPHI is evaluated relative to all-siliceous PHI types. The EP is interpreted in terms of framework ionicity for MgPHI and all-siliceous PHI models. Angular Si O Si dependence of the EFG asymmetry at 17O atoms in all-siliceous zeolites is discussed. q 2008 Wiley Periodicals, Inc.

J Comput Chem 29: 2344–2358, 2008

Key words: electrostatic potential; electric field; electric field gradient; periodic boundary conditions; hybrid B3LYP; phillipsite; cationic form

Introduction Enhanced catalytic activity of high silica aluminosilicate zeolites exchanged with divalent cations with respect to H2, CH4 was discussed both experimentally1,2 and theoretically.2 Spectroscopic evidence of H H and C H bond cleavage in H2 and CH4, respectively, were demonstrated already at moderate temperatures.1,2 The reactions of H2 and Zn-cation for some specific locations of two replacing Al atoms in the ZSM-5 fragment were considered using isolated cluster approach.3 The importance of electrostatic energy for cationic (GaZSM-5) forms is approved by emphasized intensities in vibrational spectra of ethane versus respective intensities in HZSM-5 at nearly the same coverage.2 This increase of intensities is the consequence of elevated electric field (EF) in the cationic system. The relevance of electrostatic contributions to the reaction barrier is very probable for the system with two remote Al atoms and one compensating divalent cation. This electrostatic factor cannot, however, be involved into such a consideration at the isolated cluster level. Together with different schemes already proposed to consider the influence of electrostatic potential (EP) on the barriers,4–11 cumulative coordinate technique was developed.12,13 But the level of converged basis set for its application to model accurate EP and its derivatives such as EF should be known. However, the EP or the EF convergence versus the basis set has not yet been analyzed for cationic form zeolites when considering quantum mechanical (QM) calculations.

Recently, importance of accurate basis sets with rich polarization functions was emphasized to avoid over-polarization of crystalline urea while calculating integrated Laplacian, atomic volumes, and charges calculated with AIM approach and using experimental electron densities.14 The situation remains less clear for ionic crystals, like zeolites. In this respect, important work devoted to all-siliceous zeolites15 ascertains a minor variation of their geometry starting from 6-21G basis set within framework of periodic computations. Sufficiency of 6-31G* basis set was demonstrated for optimization of isolated all-siliceous clusters capped by hydrogen atoms.16 This contradiction between the required basis set levels for accurate optimization of isolated or periodic structures should be assigned to spatial restrictions imposed by periodic conditions. Moreover, it was shown that periodic optimization using some empirical force fields (FFs) provides nearly experimental geometry for all-siliceous systems.15 The Sauer-Sierka FF17 was also shown to provide an accurate geometry relative to experiment for all-siliceous zeolites15 while Catlow FF18 is the most accurate one for aluminophosphates (AlPOs).18 On the basis of this advantage of empirical FFs, behavior of the EP, EF values was studied versus a series of eight basis sets for all-siliceous zeolites.19 In this work, one geometry optimized with the GULP software20 and Catlow FF21 was used observing minor variations of the framework geometry of all-siliceous zeolites versus the basis sets15 for a wide series of basis sets for Correspondence to: A.V. Larin; e-mail: [email protected]

q 2008 Wiley Periodicals, Inc.

Atomic Basis Sets in All-Siliceous and Mg Substituted Phillipsites

each studied zeolite. The EP and EF values were calculated in the porous space being available for adsorbate molecules as limited by negative EP values. As a result, regarding the EF, relative root mean square (RRMS) deviations of not more than 10% (TON type) and 16% (CHA type) were observed between the highest level of basis set and the lowest ones starting at the 621G*(Si)/6-31G*(O) level.19 More expanded series of thirteen basis sets was used for Hforms.22 The coordinates were fixed after empirical optimization for framework atoms with all basis sets and partially optimized for protons with some basis sets. For both all-siliceous and Hform zeolites, a linear correlation between the EP and the ionicity or the average charge of oxygen atom was observed with all basis sets with exception to the minimal STO-3G one. This correlation permits re-scaling of the EP grid for hybrid QM/molecular mechanical (MM) computations provided that the ionicity could be evaluated, for example, as proposed in ref. 23. In this work, we address to the MgPHI model to reveal the EP, EF, and electric field gradient (EFG) behavior versus basis sets of improving quality. It seems that this task suits for periodic periodic Hartree-Fock (PHF) or PDFT computations only. The size of zeolite cluster for accurate EF calculations corresponds to 7 a-units of the NaA zeolite24 or to 14 sodalite units of all-siliceous FAU type.25 Such fragment is evidently too large for high level basis set computations at the level of an isolated cluster. Hence, periodic models for zeolites are mandatory for EF calculations. Local models for EFG calculations were proposed,26 but these models cannot, however, explain the noncorrelation between Cqq(17O) and respective T O T0 angles as was approved at different computational levels of an isolated cluster27 or PHF computations.28 The periodic ab initio calculations are safer for final EFG interpretation particularly for cationic form zeolites. The proposition of one unique fixed geometry upon basis set shift as applied for all siliceous zeolite systems is, however, not reliable for cationic forms wherein Mg interaction with zeolite has essentially electrostatic character while EP varies with basis set as already shown in all-siliceous19 and H-forms.22 Hence, Mg position can vary in a greater extent upon basis set shift as compared to the Si O and Al O bond lengths. In this article, we first will present the results of the optimization at empirical level of all framework atomic positions in Mg-form PHI zeolites. Considering the optimized zeolite model, we will control the convergences of the EP, EF, and EFG using the B3LYP hybrid functional with a series of fourteen basis sets. Partial optimization of the Mg coordinates at the ab initio QM level with four basis sets will be also considered. This will allow confirming the convergence of the EF values suggesting partial ab initio optimizations. In the conclusion, we propose a simple way of transformation of the EP grid values if the EP grid was once calculated to estimate EP data corresponding to another basis set and respective zeolite ionicities.

Computational details First, we fully optimized a series of MgPHI zeolite models considering periodic conditions using GULP20 with the Catlow FF18,21 for framework atoms and Bush FF29 for Mg cation (Tables 1 and 2). Electrostatic properties, i.e., EP, EF, and EFG,

2345

˚ , degrees) and Fractional Coordinates Table 1. Lattice Parameters (A from Full Optimization Using the GULP Software with the Catlow FF for MgPHI (Si/Al 5 3) and all Siliceous PHI Models Versus Experimental Ones of Hydrated K/Ca- and Ba/Ca-Cationic Forms [30] all of them Possessing P21/m Symmetry Group. Zeolite

Parameters a, a, a, a, a,

Relaxed PHI

Mg Al Si1 Si2 Si3 O1 O2 O3 O4 O5 O6 O7 O8 O9 a, b, c, b

K2Ca1.648 (H2O)12PHIa Ba2Ca0.6 (H2O)12PHIb

b, b, b, b, b,

c, c, c, c, c,

b b b b b

MgPHI21 MgPHI14 MgPHI11 MgPHI15 MgPHI32

Values

Si1 Si2 Si3 Si4 O1 O2 O3 O4 O5 O6 O7 O8 O9 a, b, c, b

10.03491 12.77449 8.88849 130.92 10.03474 12.77404 7.92719 122.09 10.26508 13.00136 8.75381 128.48 10.45118 12.74722 8.80482 130.26 9.925933 12.197272 8.209284 129.432085 0.008935 0.250000 0.707529 0.776805 0.076440 0.269340 0.395669 0.131743 0.965599 0.060772 0.977559 0.280063 0.225721 0.125062 0.153716 0.209163 0.029889 0.276063 0.712127 0.557973 0.237120 0.594072 0.149995 0.080873 0.954015 0.875519 0.099872 0.942674 0.072786 0.253715 0.379310 0.405924 0.148338 0.872071 0.386214 0.534678 0.721075 0.750000 0.121892 0.238114 0.250000 0.232535 9.861495 14.093136 8.666541 124.680745 0.731996 0.018613 0.284970 0.420169 0.138929 0.023841 0.053759 0.017810 0.285158 0.105265 0.138852 0.025526 0.095701 0.107540 0.200087 0.637028 0.576522 0.162660 0.605028 0.108331 0.197812 0.024989 0.923099 0.161292 0.891295 0.044381 0.281640 0.292805 0.379160 0.084507 0.789840 0.498949 0.499961 0.581314 0.750000 0.028360 0.054538 0.250000 0.974806 9.865 14.300 8.668 124.20

a, b, c, b

9.879 14.139 8.693 124.81

Si/Al 5 2.019 [30]. Si/Al 5 3.0 [30].

a

b

of optimized Mg-form were then compared at the hybrid PDFT/ B3LYP (30% of Hartree-Fock exchange) level with the basis sets available in CRYSTAL98.31 The 30%-weight for HartreeFock exchange was chosen because it corresponds to intermediate ionicities expressed via Mulliken charges corresponding to B3LYP and Hartree-Fock solutions as demonstrated for a series of perovskites.32 Details of the BS1 base for Si and O are given in the CRYSTAL manual.31 For 6-21G*, the Al and Si exponents of the outer sp and d components are 0.15/0.35 and 0.11/ 0.4 au22, respectively; the O exponents are 0.42/0.72 au22. For 86-1G* (Mg) and 85-11G*(Mg), the sp and d exponents are

Journal of Computational Chemistry

DOI 10.1002/jcc

Larin, Sakodynskaya, and Trubnikov • Vol. 29, No. 14 • Journal of Computational Chemistry

2346

Table 2. Types of MgPHI Models, T-Position of Si Atom Replaced by Al in the Mg-Forms, Framework

Density (FD, g 3 cm23), PDFT Energies U 1 DU (au) Calculated at the PDFT Levels (shift by DU 5 27253 au for 85-11G* Basis Set and DU 5 22862 au for ps-21G* Basis set) Compared to the GULP ˚ ), and Mg. . .O Energies U 1 DU (au, Shift by DU 5 271 au), the Shortest Mg. . .Mg Distance (A ˚ ). Distances (A Type

T

FD

2U/85-11G*

2U/ps-21G*

2U/GULP

jMg. . .Mgj

PHI32 PHI21 PHI14 PHI11 PHI15

T(1) T(3) T(1) T(3) T(1)

2.175 1.940 1.940 1.826 1.866

0.986362 0.959437 0.959431 0.908497 0.896305

0.281782 0.260058 0.260056 0.199334 0.194209

0.417272 0.440791 0.440791 0.399265 0.395590

6.933 7.709 7.710 10.265 7.849

0.4/0.54 au22 and 0.688/0.28/0.54 au22, respectively. For 8411G*, the O sp and d exponents are 0.5/0.191/0.47 au22. For 85-11G*, the Al sp and d exponents are 0.59/0.35/0.51 au22. The 88-31G* or 8-31G* exponents for Si, Al, the 6-31G* and 6311G* exponents for O correspond to the recommended values.31,33,34 Two additional basis sets were used to characterize the relative stabililty of all five empirical MgPHI models. Pseudopotential Durand-Barthelat basis sets on T 5 Si and Al atoms was considered together with full electron basis sets for other 8511G(Mg)/ps-21G* (T)/6-21G*(O) denoted below as ps-21G*. The second basis set is 85-11G(Mg)/85-11G*(Al)/66-21G*(Si)/8411G*(O) denoted below as 85-11G*. The latter corresponds to the BS6 without d-polarization functions at Mg atom. Standard parameters have been applied for the first basis set, while the parameters of the second basis set have been described in previous publication.35 Owing to minor importance of d(Mg)-polarization functions,35 the basis set 85-11G* is very similar to the BS6 one. Standard tolerance criterions for the treatment of the two-electron integrals of 5, 5, 5, 7, 10 (TOL1-TOL5) and DFT tolerance criterions of 9, 9, 14 were applied. The SCF FMIXING convergence parameter was fixed to 30%. As standard integration grid, we used 44 points per interval as developed in CRYSTAL98, while Becke procedure was chosen as weights of the grid points.36 For the EP and EF evaluations, we choose three planes in MgPHI (Fig. 1a and 1b): Si18-O49-Si15y (grid of 90 3 107 points), O47-O31-O45 (84 3 107), and O19-Si11-O35 (95 3 100), the first one being nearly perpendicular to the vector which connects fifth coordinated O48 atom and Mg1 cation. Two other planes pass through one (third plane) or two (second plane) Mg cations. The PHI framework possesses channels in all three dimen˚ along with sions with largest effective diameters of 3.8 3 3.8 A the OX axis (Fig. 1b) that should be available for small molecules. The Mg cation can be thus in close contact with adsorbate. Two all-siliceous PHI forms were constructed. In the first one, i.e., nonrelaxed model, Mg cations were deleted and Al atoms replaced by Si at fixed framework. The second one, i.e., relaxed model, was obtained from the first PHI model by using full relaxation including cell parameters at the GULP/Catlow FF level (Table 1). The first one was only used for the analysis of

jMg. . .Oj 2*2.042, 2*2.076, 2*2.076, 2*1.924 2*2.001,

2*2.059, 2.092 2*2.112, 2.221 2*2.112, 2.221 2*2.198, 2.312

cation influence on EFG properties, while EP and EF behavior was discussed using the relaxed PHI model. When considering all-siliceous PHI models, the same atoms already selected in MgPHI are chosen as corners of the EP grids (Fig. 1c). Their atomic numbers in all-siliceous PHI are shifted by two because two Mg1 and Mg2 cations are absent. So, for the relaxed and non relaxed PHI we have three grids: Si16-O47-Si13 (89 3 111), O45-O29-O43 (84 3 112), and O17-Si9-O33 (95 3 101). Accordingly to respective three planes in MgPHI, they are also named below as first, second and third planes. The planes pass through pores wherein EP and EF values are important for modeling chemical reactivity, adsorption, and diffusion. A qualitative criterion to select areas available for adsorbate molecules can be defined by the EP values themselves.20,22 Namely, the domains corresponding to negative EP values limit the space remote from the atomic centers and hence serve to define the volume available for adsorbed species (Figs. 2a–2i). The ‘‘negative’’ EP space of the plane is estimated through the DM ratio of the number of M points with negative EP to the total number of M0 points in each plane. The M0 is slightly different for MgPHI and relaxed PHI models, i.e., 9630, 8988, 9500, and go to 9879, 9408, 9595 in the three respective planes. P Average EP (or absolute EF) values were also calculated as ( i ViBSN )/M in the considered planes. For comparison of the charge distributions obtained using different basis sets, the Mulliken atomic charges were expressed via the q0 ionicity or average oxygen charge Q00 ðOaver Þ versus the basis set levels. The dependence between the EP surfaces and q0 obtained with the different basis sets in the same plane was expressed according to the equation20,22: q0 ¼ a 3 EP þ b

(1)

where a and b are fitting parameters. To define closeness between the EP (or absolute EF) values calculated with different basis sets in the same plane, for example, between basis set BSN (N 5 1– 13) and BS14 (the largest basis), we used RRMS differences:

y

In-line notations, such as O49 in text or Figure 1, are related to the label of the atom, while subscript notations, such as in O1 in Table 1, are related to the independent crystallographic type of atom.

Journal of Computational Chemistry

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uM uP BSN u ðVi  ViBS14 Þ2 ui¼1 RRMSðBSNÞ ¼u u M P t ðViBS14 Þ2 i¼1

DOI 10.1002/jcc

(2)

Atomic Basis Sets in All-Siliceous and Mg Substituted Phillipsites

2347

1 jFyj2 1 jFzj2)1/2, Fj being the EF component along the jdirection, j 5 x, y, z). The closeness of the a slope coefficient in eq. (1) for various planes of the same zeolite19,22 suggests a possibility of recalculation of the EP surface basing on the already known EP of relatively accurate quality. This signifies a kind of parallel transfer of the EPBS14 surface not changing its curvature. The more simical lar EP shapes are, the less respective errors are. The EPBSN i;app values corresponding to the basis set BSN at any point i were shifted from the EPBSL i;cal values calculated with basis set BSL according to expression (1) via the formula 1 BSN BSL EPBSN  qBSL i;app ðxi; yi Þ ¼ EPi;cal ðxi; yi Þ þ ðq0 0 Þ a

(3)

where the qBSN and qBSL ionicities expressed via Mulliken 0 0 charges correspond to the values obtained at the BSN and BSL BSN levels. Then the RRMS difference between EPBSN app and EPcal surfaces were calculated as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uM uP BSN 2 u ðEPBSN i;app  EPi;cal Þ ui¼1 BSN u RRMSðEPcal Þ ¼ u (4) M P t 2 ðEPBSN Þ i;cal i¼1

with i from 1 to M, M being the total number of negative EP values. In eq. (3), L was chosen as 14 (for the largest basis set). We thus checked the differences between the shapes of the and other EPBSN EPBS14 cal cal surfaces. The nuclear quadrupole coupling constant Cqq characterises the quadrupole interactions of an asymmetric nuclei. The Cqq values (in MHz) within the Mg-forms were obtained here using the EFG at their respective positions37 Cqq ¼ 2:3496 3 102 qrEzz

(5)

where the coefficient in the right hand side corresponds to !Ezz expressed in e3au23, and the nuclear quadrupole moments q are 0.140238 and 20.02639 barn (1 barn 5 10228 m2) for 27Al and 17 O, respectively. Another parameter which influences the spectra of a nucleus having a quadrupolar nuclear moment is the EFG anisotropy of electric field gradient (AFG) g ¼ ðrExx  rEyy Þ=rEzz

Figure 1. (a) Mg cation local coordinates in elementary unit cell, (b) zeolite fragment versus OX axis of MgPHI and (c) elementary unit cell of PHI. Atoms chosen to define the planes wherein the EP and EF values were calculated are connected by line of the same type. The atomic labels are related to the position in the cell. The Mg O bonds are omitted in (a) and shown by dotted lines. OABC correspond to OXYZ axes.

where ViBSN are the EP (or EF) values calculated at the BSN level, with i from 1 to M, M being the total number of negative EP values in each plane at the upper BS14 level (or EF 5 (jFxj2

(6)

wherein all EFG elements are related to the EFG tensor principal axes. For the Cqq(17O) and g(17O) comparison between MgPHI and two all-siliceous PHI forms the most complete BS13 and BS14 basis sets for MgPHI were used which correspond to the BS10 for all-siliceous forms. No 25Mg analysis was performed. No experimental EFG properties for the cation in zeolites are known.

Results Geometries of Optimized Mg-Form Phillipsites

The favored Mg position was chosen according to the minimal cell energy of the MgPHI32 model being the densest one as

Journal of Computational Chemistry

DOI 10.1002/jcc

2348

Larin, Sakodynskaya, and Trubnikov • Vol. 29, No. 14 • Journal of Computational Chemistry

Figure 2. Electrostatic potential (EP, au) calculated at the PDFT/B3LYP/BS14 level in the (a) Si18O49-Si15, (b) O47-O31-O45, and (c) O19-Si11-O35 planes of MgPHI and (d, g) Si16-O47-Si13, (e, h) O45-O29-O43, and (f, i) O17-Si9-O33 planes of relaxed (g–i) and nonrelaxed (d–f) PHI. Iso-contour lines start from the EP 5 0 au level. Internal parts (in white) correspond to the domain of negative EP values wherein the EP and EF convergences were tested.

compared to other four models, all of them computed with the GULP/Catlow-Bush FF (Table 2). Presence of small Mg cation perturbs strongly the PHI framework resulting in decrease of the dehydrated cell parameters versus experimental ones. The latter

correspond to hydrated forms with different cations, which all of them are larger than Mg (Table 1). The primary Mg role for the cell contraction and not of parallel ‘‘dehydration’’ is approved by the cell parameters of relaxed all siliceous PHI form. The lat-

Journal of Computational Chemistry

DOI 10.1002/jcc

Atomic Basis Sets in All-Siliceous and Mg Substituted Phillipsites

ter are much more similar to the experimental ones30 (Table 1) while the PHI form was optimized without water molecules with the same Catlow FF.17 Relative stability of the models is similar with both ps-21G* and 85-11G* basis sets applied (Table 2). Earlier, accuracy of the GULP code for simulation of dehydrated Mg and Zn cationic forms of goosecreekite (GOO) type was confirmed by high correlation between calculated relative GULP and B3LYP energies, i.e., with correlation coefficients of 0.881 and 0.901 for series of MgGOO and ZnGOO, respectively.35 Herein, the order of GULP energies for the MgPHI32 and two isoenergetic MgPHI21 or MgPHI14 models is only inversed versus respective B3LYP energies (Table 2). Two last models differ by the positions of Mg and of the Si/Al replacement. Regarding the correlation between GULP and B3LYP energies,35 this inversion between stabilities of MgPHI32 and MgPHI21 or MgPHI14 looks to be exclusion. The deepest minimum of the MgPHI32 cell energy corresponds to five-coordinated Mg (Fig. 1a) whose coordination number follows from the analysis of relative shortest Mg O distances (Table 2). However, we will demonstrate that four nearest O atoms belong to first coordination sphere of Mg versus fifth oxygen atom on the basis of the analysis of AFG of coordinated oxygen atoms. The Mg. . .Mg distance of ˚ between the nearest Mg sites in the MgPHI32 model is 6.93 A large enough to allow the consideration of isolated Mg cations. One should mention that this is the shortest one through the series of the MgPHI models (Table 2). The geometry differences lead to the order of ionicities which does not change for all basis sets of the quality better than that of BS4 (Table 3): q0(relaxed PHI) [ q0(non-relaxed PHI) [ q0(MgPHI). Charge anisotropy of oxygen sublattice can be expressed via ratio fi ¼ Q00 ðOi Þ=Q00 ðOaver Þ (Table 4). Their qualitative analysis relative to 1, i.e., if fi [ 1 or if fi \ 1, reveals a closer similarity between both the PHI models versus MgPHI32. Only two fi values, i 5 4 and 9, from nine differ between relaxed and nonrelaxed PHI models, the first O4ss type is coordinated to Mg, while five fi values, i 5 3, 5, 6, 7, and 8, from nine are different in relaxed PHI and MgPHI32 models. Even if drastic change of the

Table 3. Average Mulliken Atomic Q00 (Oaver) Charges (e) Calculated with the Different Basis Sets at the PDFT/B3LYP Level for MgPHI and all Siliceous PHI Models.

MgPHI32 BS1 BS2 BS3 BS4 BS5 BS6 BS7 BS8 20.607 20.989 21.215 21.102 21.099 20.972 21.461 21.037 BS9 BS10 BS11 BS12 BS13 BS14 20.672 21.022 20.958 20.940 21.194 21.179 Non relaxed PHIa BS1 BS4 BS5 BS6 BS7 BS8 BS9 BS10 20.567 21.013 20.995 20.837 21.326 20.584 20.900 21.059 Relaxed PHI BS1 BS2 BS3 BS4 BS5 BS6 BS7 BS8 20.608 20.915 21.008 21.114 20.995 20.828 21.321 20.536 BS9 BS10 20.887 20.747 a

no SCF convergence is achieved for BS2 and BS3.

2349

2 Table 4. Mulliken Atomic Q00 Charges (e ) Calculated with the Upper Quality Basis Sets at the PDFT/B3LYP Level for MgPHI and all Siliceous PHI Models.

Atoms

MgPHI/BS14

PHI/BS10a

PHI/BS10b

Mg Al/Si Si1 Si2 Si3 O1 O2 O3 O4 O5 O6 O7 O8 O9 Oaver

1.685 2.119 2.152 2.184 2.134 21.052/0.898 21.228/1.042 21.199/1.017 21.200/1.018 21.199/1.017 21.072/0.909 21.319/1.119 21.255/1.064 21.067/0.905 21.179/1.000

– 2.130 2.096 2.123 2.125 21.056/0.997 21.064/1.005 21.052/0.993 21.058/0.999 21.057/0.998 21.088/1.027 21.030/0.973 21.057/0.998 21.082/1.022 21.059/1.000

– 1.492 1.495 1.492 1.494 20.7435/0.996 20.7550/1.011 20.7438/0.996 20.7547/1.011 20.7462/0.999 20.7558/1.012 20.7348/0.984 20.7399/0.991 20.7392/0.990 20.7466/1.000

For oxygens, charge ratio fi 5 Q00 /(Oi)/Q00 (Oaver) are also presented after slash. a non relaxed PHI model. b relaxed PHI model.

Si O Si angles is observed upon relaxation for the most O atoms (Table 5), the oxygen charge anisotropy varies very little between the all-siliceous PHI forms (Table 4). Models with Mg Positions Fixed Irrespective of Basis Set Level EP and EF Convergence

As the first test of the EP or EF convergence, we considered the models obtained by full optimization at the empirical level for MgPHI, without further improvement of the atomic coordinates versus the different basis sets. Considering the basis sets chosen (Table 6), one difference versus usual series for all-siliceous and H-forms should be noted. No SCF convergence was observed regarding 3-21G and 6-21G* levels, that is why the basis sets of poor quality (BS2, BS4) were obtained by omitting some polarization functions of the BS9 basis set. It is instructive to mention the advantage of the 88-31G(Si) basis set without polarization functions versus 66-21G*(Si) addressing to the energies obtained using BS6 and BS7 basis sets (Table 6). Higher importance of rich O basis set versus basis set at Si atom15 can be approved versus basis sets at the Mg and Al cations as well. Such conclusion is facilitated looking on the energies corresponding to the BS8 with best available 85-11G*(Mg)/8-31G*(Al, Si) set at the cations versus BS9 with more simple 86-1G*(Mg)/6-21G*(Al, Si) cationic basis set. Nevertheless, the corresponding 6311G*(O) basis set in the second BS9 case determines its final energetic advantage relative to 6-31G*(O) at the BS8 level. Regarding sufficient basis sets at each atomic type, the Si and O basis sets were already compared with respect to the geometry15 or versus EF value,19,22 while the basis sets for Mg and Al are not. Lower quality of the 86-1G*(Mg) than that of 85-11G*(Mg)

Journal of Computational Chemistry

DOI 10.1002/jcc

2350

Larin, Sakodynskaya, and Trubnikov • Vol. 29, No. 14 • Journal of Computational Chemistry

Table 5. T O T0 Angles (8), Quadrupole Coupling Constants Cqq(17O) (MHz), and Asymmetry g(17O) Values for MgPHI, all Siliceous Nonrelaxed and Relaxed PHI Models Calculated at the B3LYP/BS14 and B3LYP/BS10 Levels, Respectively.

Atoms

T-O-T0

Cqq

g

BS14/MgPHI O1ss O2as O3as O4ssMg O5as O6ss O7asMg O8ssMg O9ss

138.49 130.29 132.52 132.15 132.31 126.69 122.73 116.35 138.54

5.996 3.839 4.645 4.488 4.365 5.212 3.723 4.863 6.145

0.237 0.377 0.186 0.919 0.103 0.353 0.901 0.383 0.171

BS10/PHI nonrelaxed O1ss 138.49 O2ss 130.29 132.52 O3ss O4ss 132.15 132.31 O5ss O6ss 126.69 122.73 O7ss O8ss 116.35 O9ss 138.54

5.967 5.313 6.537 5.652 6.317 5.694 6.163 4.972 6.310

0.182 0.592 0.314 0.537 0.242 0.177 0.789 0.771 0.176

BS10/PHI relaxed O1ss O2ss O3ss O4ss O5ss O6ss O7ss O8ss O9ss

6.972 6.870 6.976 6.874 7.111 6.848 7.161 7.206 7.213

0.182 0.182 0.183 0.181 0.136 0.193 0.134 0.174 0.171

142.48 146.30 142.52 146.41 153.76 144.62 145.23 148.31 148.65

was demonstrated by comparison of the respective BS11 and BS12 energies (Table 6). For Al, the advantage of the 88-31G* versus 85-11G* can be clearly seen from BS10 and BS14 energies or between BS5 and BS8 (Table 6). The part of the porous space limited by negative EP values is shown for MgPHI and PHI models in all the planes (Figs. 2a–2i). The drastic variation of ionicity or average O charge with the basis level (Table 3) corroborates the absence of convergence of the EP in the series of bases (Table 7). One can, however, observe several basis set levels which provide a sufficiently accurate behavior of the RRMS values of the EF, i.e., nearly 12% starting from BS4 (Table 7), relative to the one at the BS14 level (the largest deviations are given in bold in Table 7 for several of the BS4-BS14 basis sets). The largest EF value is obtained in the Si18-O49-Si15 plane which is nearly perpendicular to the Mg1. . .O48 direction (Fig. 1a). The two other planes contain Mg cation and the nearest oxygen atoms. This space is thus screened by neighboring O atoms in higher extent than in the first plane. The upper EF error value of 12% is even smaller than 16 and 15% observed for all-siliceous (CHA)19 and H-form (HCHA, HBRE) zeolites,22 respectively. Smaller error value of

10% was obtained for all-siliceous TON framework.19 A slightly worse 20% convergence of EF was observed herein for relaxed PHI model (see Fig. 3). Unique source of experimental data on EF in zeolites is the intensities of vibrational transitions. To evaluate the EF values from infrared intensities, application of molecular gas properties, i.e., polarizability or dipole values etc., is usually adopted even if they are perturbed in adsorbed state. For molecular dihydrogen, it was shown that its polarizabilies and other molecular properties are slightly perturbed in adsorbed state so that the parameters of linear approximations of its properties are valid in the adsorbed state within 1–2%-error in a limited interval of internuclear distances.40 In the case of CO, the methods of the EF evaluations are well tabulated,41 while care should be taken regarding large monovalent cations like Rb11 whose weak EF can result in the CO localization closer to zeolite wall due to stronger dispersive interaction with framework oxygens.42 So, these gas properties cannot be a reason of the difference between the calculated and experimental EF (Table 8). A principal reason for such difference is that the experimental EF values relate to the molecular position and do not correspond to EF averaged over negative EP domain where the diatomic probe should also be. The aim of comparison of the experimental and calculated EF values is to illustrate a close order of values that cannot result in their coincidence. That is why some discrepancies versus experimental EF values are natural. Nevertheless, the average absolute EF values obtained for Mg-form (Table 8) are close to the ‘‘experimental’’ order of magnitude, i.e., from 28 3 109 V/m (or 5.5 3 1022 au) for O2 adsorbed in zeolite A exchanged with Ca21 cations43 to 2.4 3 109 V/m (0.47 3 1022 au) for N2/NaRbY,44 and repeat in general the convergence shown by the RRMSs of the EF (Table 7) when the basis set quality is better than BS3. We should also note some discrepancies between experimental EF measurements. More particularly, the smaller Mg cation45 produces a smaller EF than the larger Ca cation does.43 Even if the data are related to different zeolites and molecules, it seems to be less probable considering the possible cationic sites in the A and b type. They do not correspond to shielded or fully coordinated Ca position in beta type. We assign these contradictions to difficulties of the measurements of total adsorbate coverage which is proportional to infrared intensity value. The ionicities of both relaxed and non relaxed all-siliceous forms of the PHI type (Table 3) are larger than that of MgPHI starting from BS5. At this level EF remains within 12%-variation versus the best available one, while EF(MgPHI)/EF(PHI) ratio varies from 2.0 to 2.8 times in the same porous space (see Fig. 4). It looks like that the EF increase due to the Si/Al replacement and the addition of Mg is not local and shifts strongly the average EF values in the pores. This observation is in a drastic difference with the situation in the H-forms. We remind that the EF(HCHA)/ EF(CHA) ratio approaches 1.25 with the largest BS13 used.22 The latter emphasizes the role of cation for the EF increase. As a resume of this part, we remark that convergence of the EF value in MgPHI does not overhead 12% versus to the best available value with the BS14 basis set starting from BS4. Average EF value in relaxed PHI model reveals a worse convergence within 20% starting from BS5.

Journal of Computational Chemistry

DOI 10.1002/jcc

Atomic Basis Sets in All-Siliceous and Mg Substituted Phillipsites

2351

Table 6. Total Energies (au) Calculated with the Different Basis Sets at the PDFT/B3LYP Level for MgPHI and all Siliceous PHI Models.

Basis set notationa

Energyb

MgPHI32 BS1 BS2 BS3 BS4 BS5 BS6 BS7 BS8 BS9 BS10 BS11 BS12 BS13 BS14

STO-3G 86-1G(Mg)/6-21G(T)c/6-311G(O) 86-1G(Mg)/85-11G(Al)/6-21G(Si)/8-411G(O) 86-1G(Mg)/6-21G(T)/6-311G*(O) 85-11G*(Mg, Al)/8-31G*(Si)/6-31G*(O) 85-11G*(Mg, Al)/6-21G*(Si)/8-411G*(O) 85-11G(Mg, Al)/8-31G(Si)/8-411G*(O) 85-11G*(Mg)/8-31G*(T)/6-31G*(O) 86-1G*(Mg)/6-21G*(T)/6-311G*(O) 85-11G*(Mg, Al)/8-31G*(Si)/8-411G*(O) 86-1G*(Mg)/8-31G*(T)/6-311G*(O) 85-11G*(Mg)/8-31G*(T)/6-311G*(O) 86-1G*(Mg)/85-11G*(Al)/8-31G*(Si)/8-411G*(O) 85-11G*(Mg)/8-31G*(T)/8-411G*(O)

27162.650835 27252.845954 27253.297667 27253.492173 27253.969293/27253.969295d 27253.977082 27254.023335 27254.070583 27254.303460/27254.303471e 27254.652030 27254.741868 27254.769227/27254.769260f 27254.879366 27254.887924/27254.887940g

PHI BS1 BS2 BS3 BS4 BS5 BS6 BS7 BS8 BS9 BS10

STO-3G 6-21G(Si)/6-311G(O) 8-31G(Si)/8-411G(O) 6-21G(Si)/6-311G*(O) 8-31G*(Si)/6-31G*(O) 6-21G*(Si)/8-411G*(O) 8-31G(Si)/8-411G*(O) 6-21G*(Si)/6-311G*(O) 8-31G*(Si)/6-311G*(O) 8-31G*(Si)/8-411G*(O)

26953.129356/26952.879129 27040.456428/–h 27041.075272/–h 27041.047628/-7040.789478 27041.843270/-7041.456332 27041.690368/-7041.330555 27041.715867/-7041.395712 27042.077396/-7041.714696 27042.521530/-7042.143823 27042.338732/-7042.317448

Basis

a

The notations of core orbitals for 66-21G* and 88-31G* basis sets at Si and Al atoms are shortened as 6-21G* and 8-31G*. b ˚ optimized at the GULP/Catlow In the first (MgPHI) part, energy is given for the Mg2. . .O47 distance of 2.0918 A FF level, second values after slash correspond to jMg2. . .O47j distances optimized at the BS5, BS9, BS12, and BS14 levels; in the second (PHI) part, energy is given for relaxed/non relaxed all siliceous PHI models, respectively. c T 5 Si, Al. d ˚. jMg2. . .O47j 5 2.0847 A e ˚. jMg2. . .O47j 5 2.0854 A f ˚. jMg2. . .O47j 5 2.0829 A g ˚. jMg2. . .O47j 5 2.0984 A h no SCF convergence was achieved for non relaxed model.

Electrostatic Potential—Ionicity Correlation

The ‘‘EP – q0’’ correlation (1)19,22 allows a simple shift (3) of the EP values obtained with basis sets of relatively moderate quality to get similar type EP surfaces corresponding to more elaborated basis set22 that, a priori, could not be applied to systems because of their size. Such procedure allows easily to avoid or to minimize EP gap that appears in QM/MM schemes if two or more basis sets are applied to different parts of a complex system.7,9 For shifting calculated EP values via eq. (3) to the q0 value corresponding to another basis set for a structure under study, one should define a new ionicity value corresponding to the new basis set to be applied. To realize it, detailed charge dependencies on the atomic geometry23 could be calculated for smaller size zeolites and then be used to approximate atomic charges in upper size zeolites. The linear ionicity correlation versus the EP (1) was shown for both relaxed all-siliceous and Mg-forms (see Fig. 5). In the

both cases, the three BS1-BS3 basis sets of lower quality deviate from the linear function observed previously also for other all siliceous19 and H-form zeolites.22 Omitting BS1-BS3 data from the fitting set of the ‘‘EP – q0’’ pairs one can improve the correlation (1) which is better for cationic and worse for all-siliceous forms (Table 9). Larger ionicity value in the MgPHI type at similar basis set level starting from BS5 (Table 3) seems to be the main reason of higher correlation for MgPHI than for PHI. For all-siliceous PHI zeolite, the slope is slightly smaller than that for MgPHI. This is in coherence with smaller slope values for the CHA than that for HCHA zeolite.22 It was found earlier that the slopes are rather specific for each structure, for example, between 21.7 and 22.4 for HCHA and between 16.6 and 18.2 for HBRE.22 Herein, we can not confirm the specificity of slope for different framework with only PHI type. Another difference versus all-siliceous and H-forms is a higher EP ‘‘homogeneity’’ in Mg-forms. The average EP values are

Journal of Computational Chemistry

DOI 10.1002/jcc

Larin, Sakodynskaya, and Trubnikov • Vol. 29, No. 14 • Journal of Computational Chemistry

2352

a Table 7. RRMS Difference (%, According to eq. (2) Without Shift (3) or According to eq. (4) Using the Shift) Between the EF or EP Values Calculated with the Different Basis Sets Relative to BS14 at the PDFT/B3LYP Level Considering the Negative EP Values in Each Plane of MgPHI.

EF (2)b Plane Si18-O49-Si15 O47-O31-O45 O19-Si11-O35

EP (2)b Si18-O49-Si15 O47-O31-O45 O19-Si11-O35

EP (4) Si18-O49-Si15 O47-O31-O45 O19-Si11-O35

DM

BS1

0.552 0.510 0.622

46.35 47.31 47.30

0.552 0.510 0.622

0.552 0.510 0.622

BS2 13.18 16.56 14.57

BS3

BS4

BS5

BS6

BS7

14.69 17.93 16.26

3.89 4.89 4.50

9.55/9.51 9.68/9.71 9.83/9.86

11.35 11.48 11.10

3.36 3.11 3.25

BS8 8.22 8.49 8.64

BS9 7.43/7.36 7.90/7.86 7.83/7.88

BS10 10.31 10.60 10.67

BS11 9.50 9.98 10.03

BS12 9.19/9.14 9.49/9.52 9.58/9.62

BS13 0.22 0.45 0.50

BS1 30.82 29.95 27.74

BS2 0.25 1.92 0.88

BS3 43.81 42.66 45.40

BS4 23.80 24.57 25.58

BS5 15.37/15.11 15.85/15.56 16.32/15.29

BS6 28.92 29.57 31.14

BS12 25.95/25.75 27.09/26.77 29.41/29.43

BS13 0.57 0.54 0.53

BS8 19.42 20.35 22.25

BS9 46.72/46.57 48.24/48.01 50.27/50.31

BS10 20.78 21.41 22.15

BS11 24.77 26.26 28.62

BS1 52.43 49.74 50.21

BS2 25.05 22.61 23.23

BS3 40.97 39.19 41.01

BS4 15.34 15.91 16.66

BS5 9.59 9.82 9.72

BS6 13.84 13.94 13.74

BS10 9.92 9.89 9.42

BS11 9.94 9.93 9.48

BS12 9.38 9.21 8.84

BS13 0.90 1.05 1.20

BS8 8.89 9.03 9.38

BS9 6.64 6.38 5.95

BS7 26.72 27.05 27.41

BS7 2.59 2.90 3.66

DM is the ratio of the number M of negative EP values versus the total number of points (M0) in each plane. a maximal RRMS for each property, EP or EF, for the basis sets with the quality better than the one of BS3 is given in bold. b the values after slash correspond to Mg positions optimized at respective basis set level (see corrected Mg O distances in Table 3).

Figure 3. Convergence of average absolute electric field EF(PHI) in the 1st (circles), 2nd (triangles up), and 3rd (triangles down) planes of relaxed PHI models versus basis sets BS1-BS9 at the PDFT/ B3LYP level. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]

closer in different planes of MgPHI than that in PHI. It can be seen comparing the wider EP intervals between the average EP values for different planes for the PHI model than of MgPHI (see Fig. 5). For the latter, three (EP, q0) points for three planes are contracted to nearly one at the same basis set level. This inhomogeneity is even larger for all-siliceous CHA and H-forms of BRE and CHA.19,22 All the latter possess lower ionicity than MgPHI. One should note that a visual similarity was also emphasized between the EP iso-surface calculated at the HF/6-21G* level using the isolated cluster for 10-member window of all-siliceous silicalite and the EP iso-surface of silicalite calculated at the PHF/STO-3G level.46 To quantify EP similarity according to eq. (4) between the EP calculated with different basis sets after shift (3), the basis set BS14 or L 5 14 was chosen. It was found that the shift (3) decreases the RRMS error (4) between EPBSN cal and BSN EPBSN and app relative to the RRMS error (2) between EPcal (Table 8). The RRMS of EP values shifted via (3) are EPBS14 cal less than 20% for all basis sets starting from BS4 (Table 7). Only minor RRMS between the EP values for the models with very close q0 values, such as qBS13 (21.194 e see Table 3) and 0 qBS14 (21.179 e), slightly increases after shift (last line in Table 0

Journal of Computational Chemistry

DOI 10.1002/jcc

Atomic Basis Sets in All-Siliceous and Mg Substituted Phillipsites

2353

2 Table 8. Average EF Values (310 , au) Calculated with the Different Basis Sets at the PDFT/B3LYP Level Considering the Negative EP Values in Each Plane of PHI and MgPHI.

Type

Plane

BS1

BS2

BS3

BS4

BS5

BS6

BS7

BS8

BS9

BS10

PHIa

Si16-O47-Si13 O45-O29-O43 O17-Si9-O33

1.529 1.863 1.377

2.054 2.523 1.849

1.726 2.127 1.506

2.122 2.609 1.935

1.607 1.980 1.391

1.674 2.060 1.447

1.784 2.200 1.577

1.568 1.943 1.354

1.589 1.960 1.375

1.602 1.982 1.375

MgPHI

BS1 BS2 BS3 3.517 4.934 5.116 3.359 4.478 4.611 2.985 4.070 4.197 BS8 BS9 BS10 4.182 4.193 4.212 3.722 3.723 3.745 3.436 3.443 3.432 5.5 for O2/CaA [43], 0.47 for N2/NaRbY [44], 2.5 for CO/Mg-beta [45], 1.9 for CO/Ca-beta [45]

BS4 4.501 3.981 3.684 BS11 4.177 3.743 3.454

BS5 4.211 3.734 3.433 BS12 4.165 3.720 3.428

BS6 4.301 3.878 3.537 BS13 4.509 3.980 3.689

BS7 4.640 4.071 3.784 BS14 4.505 3.974 3.683

Si18-O49-Si15 O47-O31-O45 O19-Si11-O35

Experiment

Experimental EF values (3102, au) are related to center-of-mass position of the molecules. a BS1-BS10 basis sets as used for all-siliceous zeolite in Table 3.

7). It means that the shift (3) deletes essential part of the differBS14 ence between EPBSN only regarding different ioniccal and EPcal ities. We should mention that the RRMS of the EP (4) depends mainly on different EP derivatives or EF values at each basis set. As soon as the EF convergence within 12% was achieved starting from the BS4 level, we should thus obtain small RRMS differences for the EP at a level of basis set higher than BS4 due to convergent EF values. Comparing the first (EF (2)) and the third (EP (4)) lines in Table 7 we indeed observe a similarity between the RRMS for the basis sets of the quality better than that of BS4. The higher is the ionicity difference between two basis sets, the higher is the gain in the RRMS decrease of EP (4) due to the shift (3). For example, the differences are the largest between BS14 and BS7 (0.282 e) or between BS14 and BS9 (0.507 e). Respective gain in the RRMS of EP (4) can be crudely estimated as from 27 to 3% for BS7 and as from 47–50% to 6% for BS9. For other basis sets, the RRMS decrease due to allowance of the different ionicities is comparable with the RRMS value

Figure 4. Ratio of the average absolute electric field EF(MgPHI)/ EF(PHI) between the values in the 1st (circles), 2nd (triangles up), and 3rd (triangles down) planes of MgPHI and relaxed PHI models versus basis sets BS1-BS14 at the PDFT/B3LYP level.

due to different EF values at respective BSN and BS14 basis sets. After considering the simultaneous variation of the EF values upon the shift of basis set, this method could be farther improved for implementation within embedded cluster, ONIOM or similar QM/MM schemes in ionic solids. In these QM/MM models a QM part can be embedded into MM part while both parts are described with different basis sets. Hence, QM part is subjected to EP created by more or less ionic system than EP created by QM part itself. This could result in over- or underestimated stability of the structures with substantial charge transfer.

Figure 5. Ionicity q0(e2) defined as the average O atomic charge versus the EP (au) calculated at the PDFT/B3LYP level with different basis sets: in the 1st (circles), 2nd (triangles), and 3rd (diamonds) planes of MgPHI (filled) and relaxed PHI (open) models planes. Approximations are depicted by lines, respectively. Correlation r and slope parameter a of the function q0 5 a 3 EP 1 b (1) are given in Table 9. The data obtained with BS1-BS3 are marked by the ellipses.

Journal of Computational Chemistry

DOI 10.1002/jcc

Larin, Sakodynskaya, and Trubnikov • Vol. 29, No. 14 • Journal of Computational Chemistry

2354

Finalizing the EP analysis, we should denote that a similarity between the EF and EP behaviour versus basis set shift is found. The EP reveals also a kind of convergence versus basis sets starting from BS5, if one subtracts its part dependent on the difference of the Mulliken charges/ionicity between the chosen basis set and the one of upper quality (BS14 for MgPHI). For these basis sets, this subtraction via EP shift (3) minimizes the RRMS of EP (4) between the EP surfaces up to the deviations appearing due to the different EF values calculated with the two basis sets, i.e., the upper quality basis set studied and the chosen one. As a result, RRMS of EP (4) resembles RRMS of EF (2) for the BS5-BS14 basis sets (Table 7). Convergence of Electric Field Gradient

Different number of neighboring atomic shells for convergence of chemical shift on the one hand and EFG dependent properties, i.e., Cqq(17O) and g(17O), on the other hand was shown with Hartree-Fock and DFT approaches for isolated clusters of silica polymorphs.47 Even if the nearest atoms are shown to be sufficient for convergent Cqq(17O) and g(17O) values in the polymorphs, additional questions remain about a cationic influence. These could be checked via Cqq(17O) and g(17O) analysis in allsiliceous pseudo-structure which does not relax upon Al/Si replacement and Mg withdrawing. Such structure holds local geometry and neighbor types of all Oss atoms whose Cqq(17O) and g(17O) vary mainly due to cation deletion. Oas atoms are less convenient because neighboring Al atom is replaced by Si one so that additional factor comes into play. Keeping the geometry of cationic form zeolite upon Al/Si replacement, Cqq(17Oss) and g(17Oss) can be compared and provide a measure of cationic influence. Surely, this estimate can be realized at high enough basis set level because of a lower stability and usual poor SCF convergence of periodic solution for such non relaxed structure. 27

Al Atoms. The Cqq(27 Al) value at tetrahedral position is determined by its distortion, whose extent can be expressed through either the rows of irreducible representations of Td group using valence bonds or angles,28,48 or the deviations versus perfect tetrahedral angle.49 The extent of distortions of bond lengths of Al atom in cationic form is smaller relative to the one in the H˚ and forms. Fe, Al O bond lengths vary from 1.713 to 1.891 A Table 9. Parameters of the q0 5 a 3 EP 1 b Correlation and Coefficient r Between the EP Values and the Ionicity q0 Defined as the Average O Charge of the Framework Corresponding to the Fitting Without the BS1 and BS3 Basis Sets (r0 , a0 ) or Without the BS1-BS3 Basis Sets (r, a) at the PDFT/B3LYP Level for MgPHI Zeolite Considering the Negative EP Values in Each Plane.

Type

Plane

r0

a0

r

a

PHI

Si16-O47-Si13 O45-O29-O43 O17-Si9-O33 Si18-O49-Si15 O47-O31-O45 O19-Si11-O35

0.758 0.759 0.786 0.978 0.980 0.977

11.30 11.67 12.13 14.05 13.98 13.54

0.869 0.876 0.880 0.991 0.991 0.989

11.83 12.26 12.35 14.77 14.42 14.07

MgPHI

Figure 6. Quadrupole coupling constants (a) Cqq(27Al) (MHz) and (b) asymmetry g(27Al) values (both are shown by circles) in MgPHI versus basis sets BS1-BS14. Experimental values for different M counter ions, M 5 Al1x, H11, Na11,51 are depicted by lines.

˚ in HBRE and HCHA, respectively, while from 1.645 to 1.823 A the bond interval is narrower as much as 2 in MgPHI, spanning ˚ . As a result, more isotropic Al geometry from 1.710 to 1.797 A (g(27 Al) 5 0.231) relative to the one in H-forms leads to a smaller Cqq(27 Al) 5 6.923 MHz value (Fig. 6a) than that related to H-forms. The error of this estimation could be obtained from deviations of Cqq(27Al) which is less than 15% with any basis set of a quality better than that of BS1. This error is even smaller than was calculated for H-forms.50 Then the calculated Cqq(27Al) value of 6.923 6 1.039 MHz is between the experimental ones being 16 6 0.5 or 14 6 0.5, 14.5 6 1.0, and 5.5 6 0.5 MHz for different counter cation types H11, A1xy), and Na11, respectively (Fig. 6a), studied in dehydrated NaHY and AlNaY forms as well as in dealuminated NaHY.51 Its comparison is extremely relevant to judge on the over- or underestimated Cqq(27Al) or g(27Al) values because only one AlO4 geometry is present in MgPHI. A direct comparison of Cqq(27Al) and g(27Al) with experimental data is not possible because no NMR data for Mg-exchanged forms are available. Nevertheless, the AlO4 distortion produced by small Mg21 cation seems to have y

We omit for shortness a discussion of the 21 charge of extra-framework Al cations, showing Alx1 charge as in original work [51].

Journal of Computational Chemistry

DOI 10.1002/jcc

Atomic Basis Sets in All-Siliceous and Mg Substituted Phillipsites

2355

intermediate place between the ones produced by the H11 and Na11 cations. Mg21 cation is smaller than Na1152 and could distort dehydrated framework more pronouncedly than Na11 thus leading to higher Cqq(Al) value. For H11, it is well known experimentally that AlO4 distortions in H-forms due to covalent O H bonds are more essential than in cationic ones.51 The AlO4 distortions produced by Alx1 counter cations should be closer to the case of H11 influence, i.e., more emphasized, owing to high Alx1 charge. All three situations with three different counter ions, i.e., Na11, H11, and Alx1, were studied experimentally.51 The Cqq(27Al) comparison with the data51 does not reveal a drastic underestimation of calculated value. Hence, this AlO4 type in MgPHI can serve for the discussion of Cqq(27Al) convergence. On the contrary, the g(27Al) of 0.231 calculated herein is smaller than the values obtained with any counter cation,51 i.e., 0.8 (Na11) or 0.3 (Alx1, H11) (Fig. 6b). In the case of g(27Al), we could suspect an underestimated value. Hence, more AlO4 geometries should be considered to judge more reliably about the g(27Al) convergence. Summarizing our results on AFG, we note that the AFG deviations of both 27Al and 17O atoms versus basis set are larger than their EFG values or Cqq values (see below), i.e., g(27Al) rises up to 36% at the BS12 level (Fig. 6b). 17

O Atoms. Analysis of influence of zeolite structure on experimental NMR 17O parameters was reliable using local TownesDailey model,26 i.e., considering the influence on the parameters from only nearest neighbors. However, several experimental observations were not fully interpreted. For example, chemical 17 O shift depends on the T O T0 angle,26,53 while Cqq(17O) 27 does not at isolated cluster or periodic level.28 Minor slope value of the Cqq(17O) dependence on the T O T0 angle is 28 smaller than respective error. Even if the Townes-Dailey model predicts the Cqq(17O) increase with T O T0 angle, various functional forms were also proposed later [53 and refs. therein] to interpret experimental Cqq(17O) data. g(17O) decrease with the T O T0 angle was shown in agreement with the Townes-Dailey model.54 However, cationic influence was not discussed in details neither for Cqq(17O), nor for g(17O). That is why new study should be performed with a series of basis sets and within framework of periodic approach which considers all contributions to the Cqq(17O) and g(17O) values. We remind that earlier periodic study of H-forms addressed to two basis sets of similar quality.28 Relative variations of Cqq(17O) and g(17O) versus basis set (see Fig. 7) are similar to the case of H-forms.50 The Oss4Mg and Oas7Mg atoms coordinated to Mg possess maximum g(17O) values (Table 5) like protonated oxygens O7H in HBRE and O8H in HCHA.50 Similarly to H-forms, the largest g(17O) variations upon basis set shift are related to the smallest g(17O) values at the Oas3 and Oas5 positions (Fig. 7a). Contrary to the H-forms, the Cqq(17O) variation does not exceed the 30%- and 20%-levels starting from BS4 and BS8 (Fig. 7b), respectively, while for HBRE the 20%-level is not achieved with any basis set applied. The 30%-level for g(17O) in MgPHI is available starting from BS4 (Fig. 7a), while the 50%-deviation occurs using BS12 for HBRE. We remind, the latter was also optimized at the empirical FF level with GULP like MgPHI.

Figure 7. Variations of (a) asymmetry g(17O) and (b) quadrupole coupling constants Cqq(17O) (MHz): circles (O2as and O1ss), triangles up (O3as and O4ssMg), triangles up (O5as and O6ss), diamonds (O7asMg and O8ssMg), and crosses (O9ss) in MgPHI versus BS1-BS14 basis sets. Oas and Oss types are shown by open and closed symbols, respectively. Upper 20%- and lower 230%-limits are depicted by horizontal lines. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]

An extent of Mg cation influence on the Cqq(17O) and g(17O) values as well as of the geometry variation, was illustrated compared to those for MgPHI and nonrelaxed all-siliceous PHI forms. Essential Cqq growth at all Oas positions while going from MgPHI to the non relaxed all-siliceous PHI form points on importance of nature of neighboring atom (Si or Al) because of closeness of the Si and Al charges (Table 4) in MgPHI and non relaxed PHI. At the O4ssMg position, the increase is determined by the deletion of the Mg only (25.9%). Comparing the Cqq variations between MgPHI and nonrelaxed PHI forms for the O4ssMg (25.9%) and O7asMg (116.2%) atoms coordinated to the Mg cation with the Cqq variations for other not coordinated Oas types, i.e., O2as (38.4%), O3as (40.7%), or O5as (44.7%), one can conclude that the Al/Si replacement and Mg deletion influence on the Cqq in a similar extent for the O7asMg atom only. For other cases Al/Si replacement is much more important for Cqq than Mg withdrawing. On the contrary, Mg deletion influences on the g(O) values of all Oas and Oss types to higher extent with exception of O9ss (Table 5).

Journal of Computational Chemistry

DOI 10.1002/jcc

2356

Larin, Sakodynskaya, and Trubnikov • Vol. 29, No. 14 • Journal of Computational Chemistry

As we proposed in the beginning of part Convergence of Electric Field Gradient, the higher Mg influence on g(O) than on Cqq can be controlled more precisely with three Oss type atoms, i.e., O1ss, O6ss, and O9ss. The Mg. . .O1ss distances are 4.090 and 4.600 ˚ , while the O6ss and O9ss atoms are a bit closer to Mg, i.e., A ˚ , respectively. O9ss atom is placed at the same 3.651 and 3.360 A symmetry plane as Mg (y(O9ss) 5 0.25 in Table 1) so that a weak influence of Mg on the g(O9ss) is partly explained owing to the symmetry of the O9ss site despite of its closer position toward to Mg. The relative variations of the g(17O) at other two O1ss and O6ss positions between the MgPHI and non relaxed PHI, i.e., 23.2 and 50.0%, are in agreement with their Mg. . .O distances. The O atoms of coordination sphere of Mg can be partitioned into two groups via comparison of g(O4ssMg) and g(O7asMg) on the one hand and g(O8ssMg) on the other hand. The first two values are close to 0.9 versus 0.383 for O8ssMg which is slightly ˚ for O8ssMg vermore distant as compared to others, i.e., 2.092 A ˚ sus 2.042 and 2.059 A for O4ssMg and O7asMg (Fig. 1a). Negligible Mg role for the Cqq(O8ssMg) value and, moreover, for another Oss atom is approved by minor Cqq(O8ssMg) variation, i.e., 2.2%, between the MgPHI and nonrelaxed all-siliceous PHI form. Even if the bond length difference between O4ssMg and O8ssMg ˚ , we or between O7asMg and O8ssMg is only 0.05 and 0.033 A would consider that O8ssMg is out of the first coordination sphere of Mg. Relaxation of the PHI framework results in the Si O Si increase from 48 to 328, the decrease of g(17O), and higher Cqq(17O) values. The O atomic positions become more similar comparing the Si O Si angles, g and Cqq values (Table 5). A drastic difference between the g(17O) behavior versus the T O T0 angle either in MgPHI, or in PHI was manifested. A linear correlation with coefficient r 5 20.765 was observed between g(17O) and respective Si O Si angle unifying the calculated data for all-siliceous PHI (relaxed form), ATS,50 and TON19 zeolites (see Fig. 8). g(17O) of the O atoms in MgPHI do not obey a similar dependence. Analogous g(17O) function decreasing with Si O Si angle was observed in some silica polymorphs54 that follows from the Townes-Dailey model.26

Figure 8. Asymmetry g(17O) values versus T-O-T0 angles (8) in allsiliceous relaxed PHI (triangles up), ATS50 (triangles down), TON19 (diamonds), and MgPHI (circles). Linear approximation for all-siliceous zeolites only is depicted by dashed line.

Models with Mg positions ab initio Optimized at Chosen Basis Set Level

The possibility of easy cation exchange proves that the cations are less rigidly bonded to zeolite framework as compared to the Si and Al atoms. The interaction between cation and zeolite has essentially electrostatic character that was confirmed in terms of Laplacian sign55 at the bond critical point rc positions or through the ratio of electron potential energy V(rc) to electron kinetic G(rc) energy.56,57 Positive Laplacian !2q(rc) values allow to classify Mg O, Al O, and Si O bonds as closed shell ionic type interactions. More delicate partition into ionic (Mg O) and intermediate (Al O, Si O) types was obtained58 basing on Espinosa scheme57 in the jV(rc)j/G(rc) 2 H(rc) or jV(rc)j/ G(rc) 2 !2q(rc) coordinates, where H(rc) is the total electron energy, i.e., H(rc) 5 G(rc) 1 V(rc). High cationic charge (Table 4) makes the cation more sensitive to EP variations upon basis set shift (Table 5) so that its position can vary to a greater extent with basis. Hence, we performed ab initio optimization of the Mg site keeping the atomic positions of other framework atoms fixed using BS5, BS9, BS12, and BS14. The fixed coordinates are partly justified at relatively high level of basis set starting from 6-21G(Si, O) as found for all-siliceous zeolites.15 Special cationic position (y(Mg) 5 0.25, Table 1 or y(Ca/Ba) 5 0.2530) simplifies the optimization of the Mg site. This symmetric 5-coordinated Mg position allows an easy way to check the possible displacement while changing basis set. Displacements along with two coordinates from three selected ones to optimize the Mg position are forbidden. These two coordinates correspond to Mg2-O35 and Mg2-O50 distances (Fig. 1a). The Mg displacements toward O50 or O35 atoms are parallel to this plane containing the O31, O34, O44, and O45 and are not allowed within the spatial symmetry group P21/m. Hence, only one Mg coordinate could be changed in the direction being nearly perpendicular to the plane. The Mg2. . .O47 distance was chosen as the approximate coordinate which is allowed for the Mg motion (Fig. 1a). The Mg2-O47 vector toward fifth coordinated O47 atom is tilt with respect to the plane defined by the four closest O31, O34, O44, and O45 oxygens. The optimal Mg2. . .O47 length was determined at the minimum position of a parabolic curve fitted to 7–8 points separated by equidistant ˚ (resulting Mg2. . .O47 distances are given in steps of 0.01 A ˚ for footnotes of Table 6). The maximum difference of 0.0099 A the Mg position versus ‘‘empirical’’ minimum with GULP/Catlow FF was calculated with BS12 throughout all basis sets. For comparison, the jMg2. . .O47j differences of 0.0059 and 0.0007 ˚ were also calculated with two other basis sets are for ps-21G* A and 85-11G*, respectively, both applied for the MgPHI in Ref. 35. These small values ascertain the accurate fit of Mg location using GULP/Catlow FF. The RRMS differences for the models with optimized Mg site were recalculated for BS14 and three other basis sets, i.e., BS5, BS9, and BS12, whose RRMS of the EF values for the non-optimized models deviate less than 12% from the one obtained at the upper BS14 level (Table 7). The RRMS differences of the EP, or EF, for the partially optimized MgPHI32 are very close to the ones obtained for the nonoptimized models (Table 7). The minor variations of the RRMS values while

Journal of Computational Chemistry

DOI 10.1002/jcc

Atomic Basis Sets in All-Siliceous and Mg Substituted Phillipsites

optimizing the Mg coordinates only, which depend on the variation of basis set to a higher extent, validate the RRMS estimations upon fixed framework geometry.

Conclusions Electrostatic properties of the MgPHI zeolites were calculated at the hybrid PDFT/B3LYP level with a series of 14 basis sets. The most stable MgPHI form was chosen through the ones optimized with GULP code using Catlow-Bush force field. The EP and absolute value of EF were calculated in various planes, while EFG values were computed at atomic positions. The ratio between two sets of EP or EF values obtained with different basis sets in the same plane was expressed via RRMS differences. It was shown that the RRMS difference for the EF holds within 12% for a wide series of basis sets relative to the EF obtained with the largest one. Comparing previous deviations of 15% for EF convergence in all siliceous zeolites19 and 20% for all siliceous PHI model herein, we would recommend the last value as more reliable estimate of the EF convergence. Simultaneously, the increase of the EF due to Si replacement by Al atom was evaluated at the upper basis set level. The ratio between the average EF values in the same plane of the MgPHI and all-siliceous PHI vary between 2.0 and 2.8. We also concluded a similar EF convergence and EP variations with or without partial optimization of Mg location which is the most sensitive to the basis set shift. We have thus determined the level of accurate EF calculation, at least 6-31G* for O, 86-1G for Mg atoms and not lower than 6-21G for T 5 Al or Si. This requirement providing convergent EF values is harder than the quality required for accurate geometry optimization of all-siliceous zeolites, i.e., 6-21G(O, Si) basis set.15 EP and average oxygen charge do not demonstrate any convergence versus the basis sets in MgPHI. A simple shifting procedure was tested on the basis of already computed EP values for the EP evaluation in zeolites with a not yet used basis set. The proposed procedure is based on a linear correlation between the EP values and the q0 ionicity or average Mulliken charge of the framework oxygen. The shifting algorithm was tested by calculating the RRMS between the EP computed with the largest BS14 level and the one shifted from respective EP to another qBSN ionicity value related to another BSN basis set, N being 0 from 1 to 13. With such procedure, i.e., not changing the curvature of the EP surface, the accuracy of the EP approximation is better than 20% for the basis sets of a quality superior to the one of BS3. The remaining error is shown to be due to the difference between EF calculated either with BSN, or with BS14. Convergence of NMR parameters such as quadrupolar coupling constant and EFG asymmetry versus basis sets was tested in MgPHI. Cqq(17O) and g(17O) variations do not exceed the 30%-level starting from BS4, while for Cqq(17O) 20%-level is also available starting from BS8. A linear correlation was O Si angle in allobserved between g(17O) and respective Si siliceous zeolites. g(17O) values in MgPHI do not obey to such correlation. Partition of the O atoms coordinated to Mg was observed on the ground of the g(17O) analysis that allows to assign four O atoms to the first coordination sphere of Mg. Var-

2357

iations of Cqq(27Al) are less than 15% with basis set better than BS1. As for 17O, the g(27Al) converges worse than Cqq(27Al) raising up to 50% at the BS12 level. Earlier, we proposed a cumulative coordinate method which allows propagating the atomic multipoles up to the 4th order12,13 for the EF evaluation in zeolites and AlPOs without direct ab initio QM computations.59,60 The present work defines the level of basis set for which such a prediction can provide reliable atomic multipoles and hence EF values to be used, for example, for molecular dynamic computations for cationic form zeolites.

Acknowledgments The authors deeply acknowledge Prof. D.P. Vercauteren. The authors thank the FUNDP, FNRS-FRFC, and Lotterie Nationale (convention 2.4578.02) for the use of the Namur Interuniversity Scientific Computing Facility (ISCF) Centre.

References 1. Kazansky, V. B.; Serykh, A. Microporous Mesoporous Mater, 2004, 70, 151. 2. Kazansky, V. B.; Subbotina, I. R.; Rane, N.; Van Santen, R. A.; Hensen, E. J. M. Phys Chem Chem Phys 2005, 7, 3088. 3. Zhidomirov, G. M.; Shubin, A. A.; Kazansky, V. B.; Van Santen, R. A. Int J Quantum Chem 2004, 100, 489. 4. Greatbanks, S. P.; Sherwood, P.; Hillier, I. H.; Hall, R. J.; Burton, N. A.; Gold, I. R. Chem Phys Lett 1995, 234, 367. 5. Sierka, M.; Sauer, J. J Phys Chem B 2001, 105, 1603. 6. Day, P. N.; Pachter, R.; Gordon, M. S.; Merrill, G. N. J Chem Phys 2000, 2063, 112. 7. Dapprich, S.; Komaromi, I.; Suzie Byun, K. J Mol Struct 1999, 1, 461. 8. Zygmunt, S. A.; Curtiss, L. A.; Zapol, P.; Iton, L. E. J Phys Chem B 1944, 2000, 104. 9. Joshi, Y. V.; Thomson, K. T. J Catal 2005, 230, 440. 10. Tokmachev, A. M.; Tchougreeff, A. L. J Phys Chem A 2003, 107, 358. 11. Assfeld, X.; Rivail, J.-L. Chem Phys Lett 1996, 263, 100. 12. Larin, A. V.; Vercauteren, D. P. Int J Quantum Chem 2001, 83, 70. 13. Larin, A. V.; Parbuzin, V. S.; Vercauteren, D. P. Int J Quantum Chem 2005, 101, 807. 14. Birkedal, H.; Madsen, D.; Mathiesen, R. H.; Knudsen, K.; Weber, H.-P.; Pattisond, Ph.; Schwarzenbach, D. Acta Crystallogr A 2004, 60, 371. 15. Civalleri, B.; Zicovich-Wilson, C. M.; Ugliengo, P.; Saunders, V. R.; Dovesi, R. Chem Phys Lett 1998, 292, 394. 16. Nedelec, J. M.; Hench, L. L. J Non-Cryst Solids 2000, 277, 108. 17. Sierka, M.; Sauer, J. Faraday Discuss 1997, 106, 41.College: UK, 1992/1994. 18. Henson, N.J.; Cheetham, A.K.; Gale, J.D.; Chem. Mater., 1996, 8 664. 19. Larin, A.V.; Trubnikov, D.N.; Vercauteren, D.P.; J. Comp. Chem., 2008, 29, 130. 20. Gale, J.D.; GULP 1.3, Royal Institution/Imperial Colledge: UK, 1992/1994. 21. Schro¨der, K. P.; Sauer, J.; Leslie, M.; Catlow, C. R. A; Thomas, J. M. Chem Phys Lett 1992, 188, 320. 22. Larin, A. V.; Trubnikov, D. N.; Vercauteren, D. P. Int J Quantum Chem 2007, 107, 3137.

Journal of Computational Chemistry

DOI 10.1002/jcc

2358

Larin, Sakodynskaya, and Trubnikov • Vol. 29, No. 14 • Journal of Computational Chemistry

23. Larin, A. V.; Mortier, W. J.; Vercauteren, D. P. J Comput Chem 2007 28, 1695. 24. Cohen de Lara, E.; Nguen Tan, T. J Phys Chem 1917, 1976, 80. 25. Janssens, G. O. A.; Baekelandt, B. G.; Toufar, H.; Mortier, W. J.; Schoonheydt, R. A. J Phys Chem 1995, 99, 3251. 26. Townes, C. H.; Dailey, B. P. J Chem Phys 1949, 17, 782. 27. Bull, L. M.; Bussemer, B.; Anupold, T.; Reinhold, A.; Samoson, A.; Sauer, J.; Cheetham, A. K.; Dupree, R. J Am Chem Soc 2000, 122, 4948. 28. Larin, A. V.; Vercauteren, D. P. Int J Quantum Chem 2001, 82, 182. 29. Bush, T. S.; Gale, J. D.; Catlow, C. R. A.; Battle, P. D. J Mater Chem 1994, 4, 831. 30. Rinaldi, R.; Pluth, J.J.; Smith, J.V.; Acta Cryst. B, 1974, 30, 2426. 31. Saunders, V. R.; Dovesi, R.; Roetti, C.; Causa`, M.; Harrison, N. M.; Orlando, R.; Zicovich-Wilson, C. M. CRYSTAL98 1.0, User’s Manual, Torino, 1999. 32. Cora`, F.; Alfredsson, M.; Mallia, G.; Middlemiss, D. S.; Mackrodt, W. C.; Dovesi, R.; Orlando, R. Struct Bond 2004, 113, 171. 33. Zicovich-Wilson, C. M.; Dovesi, R. J Phys Chem B 1998, 102, 1411. 34. Civalleri, B.; Ugliengo, P. J Phys Chem 2000, 104, 9491. 35. Larin, A. V.; Semenyuk, R. A.; Trubnikov, D. N.; Vercauteren, D. P. Russ J Phys Chem A 2007, 81, 495 ( transl. from Zur. Fiz. Khim. A, 2007, 81, 583), available on site http://dx.doi.org/10.1134/ S0036024407040012. 36. Becke, A. D. J Chem Phys 1988, 88, 1053. 37. Engelhardt, G. Stud Surf Sci Catal 1989, 52, 151. 38. Sundholm, D.; Olsen, J. Phys Rev Lett 1992, 68, 927. 39. Ajzenberg-Selone, F. Nucl Phys A 1970, 281, 1. 40. Larin, A. V.; Jousse, F.; Leherte, L.; Vercauteren, D. P. Chem Phys Lett 1997, 274, 345. 41. Pacchioni, G.; Cogliandro, G.; Bagus, P. Int J Quantum Chem 1992, 42, 1115.

42. Larin, A. V.; Vercauteren, D. P.; Lamberti, C.; Bordiga, S.; Zecchina, A. Phys Chem Chem Phys 2002, 4, 2424. 43. Jousse, F.; Cohen de Lara, E. J Phys Chem 1996, 100, 233. 44. Marra, G. L.; Fitch, A. N.; Zecchina, A.; Ricchiardi, G.; Salvalaggio, M.; Bordiga, S.; Lamberti, C. J Phys Chem B 1997, 101, 10653. 45. Li, P.; Xiang, Y.; Grassian, V. H.; Larsen, S. C. J Phys Chem B 1999, 103, 5058. 46. White, J. C.; Hess, A. C. J Phys Chem 1993, 97, 8703. 47. Xianyu Xue, Kanzaki;, M. Solid State NMR 2000, 15, 245. 48. Jeanvoine, Y.; Angyan, J. G.; Kresse, G.; Hafner, J. J Phys Chem B 1998, 102, 5573. 49. Ghose, S.; Tsang, T. Am Mineral 1973, 58, 748. 50. Larin, A. V.; Trubnikov, D. N.; Vercauteren, D. P. in preparation. 51. Jian Jiao, Kanellopoulos, J.; Wei Wang, Ray, S. S.; Foerster, H.; Freude, D.; Hunger, M. Phys Chem Chem Phys 2005, 7, 3221. 52. Shannon, R. D. Acta Cryst 1976, A32, 751. 53. Loeser, T.; Freude, D.; Mabande, G. T. P; Schwieger, W. Chem Phys Lett 2003, 370, 32. 54. Clark, T. M.; Grandinetti, P. J.; Florian, P.; Stebbins, J. F. J Phys Chem B 2001, 105, 12257. 55. Bader, R. F. W. Atoms in Molecules; Oxford Science Publications: Oxford, UK, 1990. 56. Cremer, D.; Kraka, E. Croat Chem Acta 1984, 57, 1259. 57. Espinosa, E.; Alkorta, I.; Elguero, J.; Molins, E. J Chem Phys 2002, 117, 5529. 58. Gibbs, G. V.; Cox, D. F.; Crawford, T. D.; Rosso, K. M.; Ross, N. L. J Chem Phys 2006, 124, 084704. 59. Larin, A. V.; Hansenne, C.; Trubnikov, D. N.; Vercauteren, D. P. Int J Quantum Chem 2005, 105, 839. 60. Larin, A. V.; Trubnikov, D. N.; Vercauteren, D. P. Rus J Phys Chem 2005, 79, 1709. (transl. from Zur. Fiz. Khim., 1927, 2005, 79).

Journal of Computational Chemistry

DOI 10.1002/jcc