- Email: [email protected]
Application of the DFT-based embedding scheme using an explicit functional of the kinetic energy to determine the spin density of Mg+ embedded in Ne and Ar matrices
17 September 1999
Chemical Physics Letters 311 Ž1999. 87–92 www.elsevier.nlrlocatercplett
Application of the DFT-based embedding scheme using an explicit functional of the kinetic energy to determine the spin density of Mgq embedded in Ne and Ar matrices Tomasz Adam Wesołowski
)
UniÕersite´ de GeneÕe, de Chimie Physique, 30 quai Ernest-Ansermet, CH-1211 GeneÕe 4, Switzerland ` Departement ´ Received 12 January 1999; in final form 24 June 1999
Abstract The formalism of the Kohn–Sham equations with constrained electron density is extended to the spin-polarized case. The isotropic hyperfine coupling constants Ž A isoŽMg.. of Mgq embedded in a Ne or Ar matrix represented using a cluster are calculated and compared to that of free Mgq. For the largest basis set used, the calculated values Ž222.9 and 210.4 gauss for Ar and Ne, respectively. agree with experimental measurements Ž222.4 and 211.6.. The shifts of A isoŽMg. relative to the values for free Mgq are rather basis-set-independent. q 1999 Elsevier Science B.V. All rights reserved.
1. Introduction The determination of the hyperfine structure of an isolated radical by means of quantum-chemistry methods usually involves computer-time-expensive calculations. Within the framework of conventional ab initio approaches, reliable results can be obtained provided the applied theoretical approach assures a good description of the electronic correlation. Moreover, since the isotropic hyperfine coupling constant is proportional to the spin density at the nuclear cusp, accurate determination of the hyperfine structure implies more severe criteria for the basis set completeness than the corresponding calculations of the energy and other global observables. Therefore, the application of conventional ab initio methods is limited to relatively small systems w1x. ) Fax: q41-22-702-6518; e-mail: [email protected]
Density functional theory ŽDFT. in its Kohn– Sham ŽKS. formulation w2x has been successfully applied for deriving the hyperfine structure parameters Žfor review, see Refs. w1,3,4x.. Compared to conventional post-SCF ab initio methods, the KS calculations make it possible to include electron correlation effects on the spin density at lower computational cost and to converge with the basis set more rapidly. As a consequence, much larger systems can be studied. In the case of a radical forming a complex with other molecules, both the KS formalism as well as the conventional ab initio methods involve the quantum-mechanical treatment of the whole system comprising the radical and the molecules Žor atoms. representing the environment. In this work, an alternative strategy for deriving the spin density of an embedded molecule is proposed. Instead of applying KS theory to the whole system, the spinrelectron densities of subsystems are treated separately. The subsystems correspond to
0009-2614r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 Ž 9 9 . 0 0 7 4 5 - 9
88
T.A. Wesołowskir Chemical Physics Letters 311 (1999) 87–92
the radical and to the embedding molecules, respectively. For a given electron density of the embedding molecules, the spin density of the radical is determined using the spin-polarized version of the KS equations with constrained electron density ŽKSCED. w5,6x. Similarly to the spin non-polarized case, the KSCED method makes it possible to derive the spin density of the embedded radical by means of the KS-like equations in which the effective potential includes one part which describes the isolated radical and the other one arising from the interactions with the embedding molecules. The interaction part of the one-electron KSCED effective potential includes the non-additive kinetic energy component which is approximated using an explicit functional of the kinetic energy. This is an essential feature of KSCED. Our previous studies on molecular complexes formed by closed-shell molecules showed that the KSCED energies and KSCED electron densities depend critically on the approximate functional used to derive the non-additive kinetic energy w7,8x. The best agreement between the supermolecular KS energies and the KSCED ones was obtained with the gradient-dependent kinetic energy functional ŽPW91K. which uses the same gradient dependency as the exchange functional of Perdew and Wang ŽPW91. w9x. Indeed, the KSCED calculations using the PW91K kinetic energy functional combined with the PW91 exchange-correlation functional led to the interaction energies of other weak complexes such as the ones formed by benzene and small diatomics ŽN2 , O 2 , and CO. w10x agreeing within 0.1 kcalrmol with experimental measurements. Subsequently, the KSCED equations using the PW91K kinetic energy were applied to study the properties of CO physisorbed on metal oxide surfaces ŽMgOŽ100. and ZnOŽ1010.. w11,12x. In these studies, the KSCED equations were used as a first-principle embedding scheme in which the surface was represented by means of a cluster with the electron density calculated with the absence of CO and kept frozen for all considered geometries of the adsorbed molecule. The shifts of the CO stretching frequencies agreed within a few cmy1 with the experimental measurements. The geometry of CO and the stretching frequency were not affected if the electron density was allowed to relax. In this work, the KSCED effective potential which applies the PW91K kinetic energy functional is ap-
plied for deriving the spin density of an embedded radical. For the studied complexes formed by the Mgq cation and rare-gas atoms, the effect of the surrounding atoms on the radical’s spin density arise mainly from the Pauli repulsion between the overlapping electron densities. The supermolecule KS calculations by Eriksson w13x showed that the hyperfine structure of Mgq involved in complexes with rare-gas atoms can be well described at this level of theory.
2. Methods 2.1. The formalism The derivation of the KSCED equations in the spin-non-polarized case can be found elsewhere w5,6x. The spin-polarized version is outlined below. The spin density of the whole complex is represented as a sum:
r s s r 1s q r 2s ,
s s a or b ,
Ž 1.
where the densities r 1 and r 2 integrate to the total number of electrons N1 and N2 , respectively. Similarly to the KS theory, the electron density r 1s is represented by means of one-electron functions: N1 s < s <2 r 1s s Ý nŽ1.i cŽ1.i ,
s nŽ1.i s 0 or 1 .
Ž 2.
i
For a given r 2s, r 1s which minimizes the total energy, the functional Ew r a , r b x is derived from the following one-electron equations: s ,eff s y 12 = 2 q V KSCED Ž r. Ž r 1a , r 1b , r 2a , r 2b . cŽ1.i s s s ´Ž1.i cŽ1.i Ž r. ,
is1, N1 ,
Ž 3.
s,eff with the one-electron effective potential V KSCED taking the form: s ,eff V KSCED s
ZA
Ý y < ryR A
A
<
q
H
r 1 Ž rX . q r 2 Ž rX . d rX < rX y r <
q Vxcs Ž r 1a q r 2a , r 1b q r 2b . q
d Tsnadd r 1a , r 1b ; r 2a , r 2b dr 1s
,
Ž 4.
T.A. Wesołowskir Chemical Physics Letters 311 (1999) 87–92
where: Ø r i s r ia q r ib , i s 1,2; Ø Tsnadd w r 1a , r 1b ; r 2a , r 2b x s Ts w r 1a q r 2a , r 1b q r 2b x y Ts w r 1a , r 1b x y Ts w r 2a , r 2b x where Ts w r s x denotes the kinetic energy of non-interacting electrons in a reference system as defined in the KS theory w2x; Ø Vxcs is conventional exchange-correlation potential Vxcs s Ž d Exc . r Ž dr s . ; Ø the index A runs through the nuclei of subsystems 1 and 2. Eqs. Ž3. and Ž4. provide the spin-polarized version of equation Ž14. of Ref. w5x and equation Ž20. of Ref. w6x which where derived for overlapping closeshell systems. s,eff . The KSCED effective potential Ž V KSCED given by Eq. Ž4. can be formally partitioned into two components s ,eff s ,eff ,0 s ,eff ,int V KSCED s V KSCED q V KSCED ,
Ž 5.
where s,eff,0 Ø VKSCED does not depend on r 2 and it can be identified with the conventional KS effective potential of the isolated subsystem 1: s ,eff ,0 VKSCED s VKsS,eff Ž r 1 . s Ý y
A
H
q
ZA < r yRA <
r 1 Ž rX . d rX q Vxcs Ž r 1a , r 1b . ; < rX y r <
Ž 6.
s,eff,int Ø VKSCED collects all the terms arising from the interactions between subsystems:
s ,eff ,int VKSCED s
ZB
Ý y < ryR B
B
<
q
H
r 2 Ž rX . d rX < rX y r <
q Vxcs Ž r 1a q r 2a , r 1b q r 2b . yVxcs Ž r 1a , r 1b q
d Tsnadd r 1a , r 1b ; r 2a , r 2b dr 1s
89
For a given spin density Ž r 1s ., derived from the KSCED calculations the isotropic hyperfine coupling constant Ž A iso ŽMg.. is calculated as: 4p A iso Ž Mg . s gb g Mg b Mg² Zs :y1 3 = r 1a Ž r Mg . y r 1b Ž r Mg . , Ž 8. where Ø g and g Mg are electron and nuclear magnetons; Ø b and b Mg are electron and nuclear gyromagnetic ratios. 2.2. The embedding electron density: r 2 Since the electron density of the whole complex is represented using the electron densities of subsystems, it is possible to apply different levels of theory Žor computational effort. to derive each of them, depending on the a priori knowledge of the system under investigation. This offers the possibility of introducing additional approximations or simplifications. The embedding electron density r 2 is used for expressing analytically the embedding effective potential. Following Eq. Ž3., any change of r 2 involves a change of r 1. In this work, two ways for obtaining r 2 were applied. In the first, r 2 was taken from the KS calculations for the isolated embedding subsystem. In the second, it was obtained from the ‘freezeand-thaw’ cycle w14x in which r 1 and r 2 exchange their roles in subsequent KSCED calculations until they converge at values minimizing the total KSCED energy. The first way implies that the relaxation of the electron density of the embedding subsystems resulting from the interactions is neglected. In this work, it was assumed that r 2a s r 2b , i.e., the embedding electron density is not spin-polarized. 2.3. Localization of the electron density of subsystems
;
Ž 7.
Ø the indexes A and B run through the nuclei of subsystems 1 and 2, respectively. The terms in Eq. Ž7. represent the nuclear attraction, electron–electron repulsion, exchange-correlation, and the Pauli repulsion between the overlapping electron densities, respectively.
As has been discussed elsewhere w8x, for each subsystem the electron density Ž r i . can be localized near its nuclei or it can be allowed to spread over the whole complex. The former case offers substantial computational savings owing to the fact that a smaller number of basis set functions is needed to expand the electron density of the selected subsystem than is needed for the whole complex. In such a case,
90
T.A. Wesołowskir Chemical Physics Letters 311 (1999) 87–92
solving Eqs. Ž3. and Ž4. requires a computational effort comparable to that needed to obtain the spin density of the isolated Ži.e., not interacting with other molecules. subsystem. Since the basis sets used for expanding the fragments’ electron densities are not complete, the use of the localized basis sets is not well suited for describing such effects as the charge transfer between the subsystems or the formation of covalent bonds between them. To take into account such effects, the basis sets must be complete or at least include all functions describing the complex correctly by means of the supermolecular calculations. Such KSCED calculations involve, however, a computational effort exceeding that of the KS calculations for the whole complex. It is therefore not recommended for practical use. It is applied here to illustrate the effect that the incompleteness of the basis sets used to expand the fragment’s electron density has on the KSCED results. 2.4. Approximate functionals The PW86 w15x exchange and the P86 w16x correlation functionals were used to approximate the exchange-correlation functional in both the KS and in the KSCED calculations. Previous supermolecular KS studies of the same systems by Eriksson w13x showed that this approximation of the exchange-correlation functional is sufficiently accurate to describe the spin density. The non-additive kinetic energy part of the effective potential wŽ d Tsnadd w r 1a , r 1b ; r 2a , r 2b x.rŽ dr 1s .x was derived analytically w7x using the gradient-dependent kinetic energy functional. The used approximate kinetic energy functional was constructed following the Lee, Lee and Parr scheme w17x in which the same gradient-dependency as the one of the PW91 exchange functional w9x was used. 2.5. Numerical implementation The calculations were performed using a modified version of the deMon w18x program into which the KSCED formalism was implemented. A fine grid was used with 128 radial points in all calculations. The calculations were made using DZVP-type functions basis functions applied previously by Eriksson
in studies of the same system w13x. DZVP basis sets for neon Žcontraction pattern: Ž621r41r1) .. and TZVP for argon Žcontraction pattern: Ž73111r 6111r1) . were used in all calculations whereas three different basis sets were used for magnesium: Ø basis A ŽDZVP with contraction pattern 6321r 411)r1) ; Ø basis B ŽDZVP with decontracted valence orbitals and the contraction pattern 411111111r21111r1; Ø basis C Žfully decontracted DZVP basis sets: 111111111111r111111r1. The auxiliary basis functions for fitting the charge density and Vxc were derived from the original functions used by Eriksson w13x for the same systems Ž5,2;5,2. for Ne, and Ž5,4;5,4. for Mg and Ar. The modified functions, which can be labeled using the same convention as: Ž5,2;10,4. for Ne, and Ž5,4;10,8. for Mg and Ar, retained the original charge density part but the exchange-correlation part was supplemented by the functions fitting Ž d Tsnadd w r 1a , r 1b ; r 2a , r 2b x.rŽ dr 1s .. The functions used for fitting the kinetic energy part were obtained by a simple procedure in which all the coefficients of the original exchange-correlation functions were multiplied by two. The multiplication factor is justified by the r 2r3 dominant power dependence of the kinetic energy part of the effective potential. For each complex, the magnesium isotropic hyperfine coupling constants were calculated using a cluster comprising eight atoms to represent the matrix. The cluster atoms occupy the corners of the cube and the Mg–X ŽX s Ne or X s Ar. distance corresponds to the minimum at the supermolecular ˚ R MgAr KS potential energy surface Ž R MgNe s 2.97 A, ˚ s 3.39 A. w13x. 3. Results Table 1 shows the isotropic hyperfine coupling constants for Mgq embedded in the cluster of eight rare-gas atoms ŽNe or Ar.. The electron density Ž r 2 . taken from the KS calculations for the isolated X 8 ŽX s Ne or Ar. cluster was used to derive the ems,eff,int . bedding effective potential Ž V KSCED . The absolute values of the isotropic hyperfine coupling constant for free Mgq and embedded Mgq depend strongly on the basis set used for magnesium. The effect of the matrix on A iso ŽMg., as measured
T.A. Wesołowskir Chemical Physics Letters 311 (1999) 87–92
91
Table 1 Isotropic hyperfine coupling constants Ž A iso ŽMg. in gauss. of free Mgq ŽKS. and embedded Mgq ŽKSCED: non-relaxed electron density of the embedding atoms and the localized expansion of the electron density of subsystems.. The basis sets A, B, and C, correspond to different decontractions of the DZVP basis for magnesium Basis set
A iso ŽMgq .
A iso ŽMg8Neq .
A iso ŽMg8Arq .
A iso ŽMg8Arq . –A iso ŽMg8Neq .
A B C Exp w19x.
y197.3 y216.1 y209.1
y209.4 y229.9 y222.9 y222.4
y197.2 y216.8 y210.4 y211.6
12.2 13.1 12.5 10.8
by the difference between A iso ŽMg.Žembedded. and A iso ŽMg.Žfree., is not affected by decontracting the basis set. As a consequence, the difference between the isotropic hyperfine coupling constants of Mgq embedded in either a Ne or Ar matrix – Ne Ž Ž D AAr Mg. s A iso ŽMg8Arq–A iso ŽMg8Neq. does iso not vary with the change of the basis sets Ž12.2–13.1 gauss. and it agrees to within about 2 gauss with the experimental value Ž10.8 gauss w19x.. For the largest basis set Žfully decontracted DZVP., the absolute values of A iso ŽMg. are in excellent agreement the experimental measurements. 3.1. The effect of the relaxation of r 2 on the Mg q spin density The calculations described in the previous section do not take into account the modification of the electron density of the embedding molecules on the
calculated isotropic hyperfine coupling constant. To include such deformations, the spin density of Mgq was calculated using the r 2 obtained by means of the ‘freeze-and-thaw’ cycle of KSCED w14x. Such a cycle converges quickly for the Ne and Ar complexes studied. The values of A iso ŽMg. obtained in the third iteration differ by less than 0.1 gauss from the ones derived from the fully converged ‘freezeand-thaw’ cycle. The values of A iso ŽMg. derived using relaxed r 2 are shown in Table 2. It can be seen that the r 2 relaxation involves a very small Žusually - 1 gauss. change of the calculated values of A iso ŽMg.. 3.2. The effect of the charge transfer between the subsystems on the Mg q spin density Table 2 shows the isotropic hyperfine coupling constants derived using the supermolecular expansion of both r 1 and r 2 . It can be seen that allowing
Table 2 Isotropic hyperfine coupling constants Ž A iso in gauss. of embedded Mgq derived from the KSCED calculations using Ža. relaxed embedding electron density Ž r 2 ., and Žb. supermolecular expansion of the spin density of the subsystems. The basis sets A, B, and C correspond to different decontractions of the DZVP basis for magnesium. Values in parentheses denote the A iso changes with respect to the KSCED results derived using non-relaxed r 2 and localized expansion of r s1 and r s2 Basis set
A iso ŽMg8Neq .
Relaxation of r 2 : A y209.0 Žq0.4. B y229.5 Žq0.4. C y222.7 Žq0.2. Supermolecular expansion of r 1 and r 2 : A y207.2 Žq2.2. B y226.1 Žq3.8. C y219.2 Žq3.7. Exp w19x.
y222.4
A iso ŽMg8Arq .
A iso ŽMg8Arq . –A iso ŽMg8Neq .
y196.9 Žq0.3. y215.6 Žq1.2. y209.8 Žq0.6.
12.1 Žq0.1. 13.7 Žq0.8. 12.9 Žq0.4.
y197.0 Žq0.2. y214.9 Žq1.9. y207.7 Žq2.7.
10.2 Žy2.0. 11.2 Žy1.9. 11.5 Žy1.0.
y211.6
10.8
92
T.A. Wesołowskir Chemical Physics Letters 311 (1999) 87–92
the spin density to spread over the whole complex reduces the magnitude of A iso ŽMg. in all cases. The effect is larger for the neon complex than for that of – Ne Ž argon leading to the reduction of D AAr Mg.. iso
4. Conclusions The difference between the magnesium isotropic hyperfine coupling constants measured in Ne or Ar – Ne Ž matrices Ž D AAr Mg.. is well reproduced by iso means of the KSCED effective potential which is expressed as a function of the electron density of the atoms surrounding the radical. The KSCED calculations which use a non-relaxed electron density of embedding atoms and localized spin density of the radical lead to the values of the isotropic hyperfine coupling constants which agree to within 1 gauss with experiment for the largest basis set used. Such calculations involve a substantial reduction of computational effort compared to the supermolecular KS calculations. The absolute values of A iso ŽMg. depend strongly on the size of the basis orbital basis set used for magnesium. The relative value obtained as the difference between the A iso ŽMg. calculated for the complex and for the isolated Mgq which measures the effect of the matrix is rather basis set insensitive. – Ne Ž Similarly, the difference D AAr Mg. does not sigiso nificantly vary with the change of the magnesium basis set. The effects of the relaxation of the embedding electron density and the charge transfer between the interacting subsystems were also determined within the framework of the KSCED calculations. For the studied systems, they are not significant amounting to less than 20% of the calculated shifts of A iso ŽMg.. The results show that the Pauli repulsion between weakly overlapping electron densities can be accurately described by means of the functional derivative of the non-additive kinetic energy which is approximated using the gradient-dependent explicit energy functional. These results are in line with our previous analyses w8,7x of the applicability of the approximate kinetic energy functionals of the Lee, Lee, and Parr w17x form which indicated that for such functionals the PW91 enhancement factor is the best
approximation of the exact Tsnadd functional for small overlaps between electron densities of subsystems. Acknowledgements The author is grateful to Prof. J. Weber helpful for discussions, to Prof. D.R. Salahub for providing a copy of the deMon program, and to Prof. V. Malkin and Dr. O. Malkina for the subroutines for calculating the hyperfine tensor. Financial support from the Federal Office for Education and Science, acting as Swiss COST office, is gratefully acknowledged. This work is also a part of Project 20-49037.96 of the Swiss National Science Foundation. References w1x V. Barone, in: D. Chong ŽEd.., Recent Advances in Density Functional Theory, vol. I, World Scientific, Singapore, 1995, pp. 287–334. w2x W. Kohn, L.J. Sham, Phys. Rev. 140 A Ž1965. 1133. w3x V.G. Malkin, O.L. Malkina, L.A. Eriksson, D.R. Salahub, in: P. Politzer, J.M. Seminario ŽEds.., Theoretical and Computational Chemistry, vol. I, Density Functional Calculations, Elsevier, 1995. w4x L.A. Eriksson, in: M. Springborg ŽEd.., Density Functional Methods: Applications in Chemistry and Materials Science, Wiley, 1997. w5x P. Cortona, Phys. Rev. B 44 Ž1991. 8454. w6x T.A. Wesołowski, A. Warshel, J. Phys. Chem. 98 Ž1993. 5183. w7x T.A. Wesołowski, H. Chermette, J. Weber, J. Chem. Phys. 105 Ž1996. 9182. w8x T.A. Wesołowski, J. Chem. Phys. 106 Ž1997. 8516. w9x J.P. Perdew, Y. Wang, in: P. Ziesche, H. Eschrig ŽEds.., Electronic Structure of Solids’91, Academie, Berlin, 1991, p. 11. w10x T.A. Wesołowski, Y. Ellinger, J. Weber, J. Chem. Phys. 108 Ž1998. 6078. w11x N. Vulliermet, T.A. Wesołowski, J. Weber, Collabor. Czech. Acad. Sci. 63 Ž1998. 1447. w12x T.A. Wesołowski, N. Vulliermet, J. Weber, J. Mol. Struct. ŽTHEOCHEM. 458 Ž1999. 151. w13x L.A. Eriksson, J. Chem. Phys. 103 Ž1995. 1050. w14x T.A. Wesołowski, J. Weber, Chem. Phys. Lett. 248 Ž1996. 71. w15x J.P. Perdew, Y. Wang, Phys. Rev. B 33 Ž1986. 8800. w16x J.P. Perdew, Phys. Rev. B 33 Ž1986. 8822. w17x H. Lee, C. Lee, R.G. Parr, Phys. Rev. A 44 Ž1991. 768. w18x A. St-Amant, Ph.D. Thesis, Universite´ de Montreal, ´ 1992. w19x L.B. Knight Jr., C.B. Cleveland, R.F. Frey, E.R. Davidson, J. Chem. Phys. 100 Ž1994. 7867.
Chemical Physics Letters 311 Ž1999. 87–92 www.elsevier.nlrlocatercplett
Application of the DFT-based embedding scheme using an explicit functional of the kinetic energy to determine the spin density of Mgq embedded in Ne and Ar matrices Tomasz Adam Wesołowski
)
UniÕersite´ de GeneÕe, de Chimie Physique, 30 quai Ernest-Ansermet, CH-1211 GeneÕe 4, Switzerland ` Departement ´ Received 12 January 1999; in final form 24 June 1999
Abstract The formalism of the Kohn–Sham equations with constrained electron density is extended to the spin-polarized case. The isotropic hyperfine coupling constants Ž A isoŽMg.. of Mgq embedded in a Ne or Ar matrix represented using a cluster are calculated and compared to that of free Mgq. For the largest basis set used, the calculated values Ž222.9 and 210.4 gauss for Ar and Ne, respectively. agree with experimental measurements Ž222.4 and 211.6.. The shifts of A isoŽMg. relative to the values for free Mgq are rather basis-set-independent. q 1999 Elsevier Science B.V. All rights reserved.
1. Introduction The determination of the hyperfine structure of an isolated radical by means of quantum-chemistry methods usually involves computer-time-expensive calculations. Within the framework of conventional ab initio approaches, reliable results can be obtained provided the applied theoretical approach assures a good description of the electronic correlation. Moreover, since the isotropic hyperfine coupling constant is proportional to the spin density at the nuclear cusp, accurate determination of the hyperfine structure implies more severe criteria for the basis set completeness than the corresponding calculations of the energy and other global observables. Therefore, the application of conventional ab initio methods is limited to relatively small systems w1x. ) Fax: q41-22-702-6518; e-mail: [email protected]
Density functional theory ŽDFT. in its Kohn– Sham ŽKS. formulation w2x has been successfully applied for deriving the hyperfine structure parameters Žfor review, see Refs. w1,3,4x.. Compared to conventional post-SCF ab initio methods, the KS calculations make it possible to include electron correlation effects on the spin density at lower computational cost and to converge with the basis set more rapidly. As a consequence, much larger systems can be studied. In the case of a radical forming a complex with other molecules, both the KS formalism as well as the conventional ab initio methods involve the quantum-mechanical treatment of the whole system comprising the radical and the molecules Žor atoms. representing the environment. In this work, an alternative strategy for deriving the spin density of an embedded molecule is proposed. Instead of applying KS theory to the whole system, the spinrelectron densities of subsystems are treated separately. The subsystems correspond to
0009-2614r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 Ž 9 9 . 0 0 7 4 5 - 9
88
T.A. Wesołowskir Chemical Physics Letters 311 (1999) 87–92
the radical and to the embedding molecules, respectively. For a given electron density of the embedding molecules, the spin density of the radical is determined using the spin-polarized version of the KS equations with constrained electron density ŽKSCED. w5,6x. Similarly to the spin non-polarized case, the KSCED method makes it possible to derive the spin density of the embedded radical by means of the KS-like equations in which the effective potential includes one part which describes the isolated radical and the other one arising from the interactions with the embedding molecules. The interaction part of the one-electron KSCED effective potential includes the non-additive kinetic energy component which is approximated using an explicit functional of the kinetic energy. This is an essential feature of KSCED. Our previous studies on molecular complexes formed by closed-shell molecules showed that the KSCED energies and KSCED electron densities depend critically on the approximate functional used to derive the non-additive kinetic energy w7,8x. The best agreement between the supermolecular KS energies and the KSCED ones was obtained with the gradient-dependent kinetic energy functional ŽPW91K. which uses the same gradient dependency as the exchange functional of Perdew and Wang ŽPW91. w9x. Indeed, the KSCED calculations using the PW91K kinetic energy functional combined with the PW91 exchange-correlation functional led to the interaction energies of other weak complexes such as the ones formed by benzene and small diatomics ŽN2 , O 2 , and CO. w10x agreeing within 0.1 kcalrmol with experimental measurements. Subsequently, the KSCED equations using the PW91K kinetic energy were applied to study the properties of CO physisorbed on metal oxide surfaces ŽMgOŽ100. and ZnOŽ1010.. w11,12x. In these studies, the KSCED equations were used as a first-principle embedding scheme in which the surface was represented by means of a cluster with the electron density calculated with the absence of CO and kept frozen for all considered geometries of the adsorbed molecule. The shifts of the CO stretching frequencies agreed within a few cmy1 with the experimental measurements. The geometry of CO and the stretching frequency were not affected if the electron density was allowed to relax. In this work, the KSCED effective potential which applies the PW91K kinetic energy functional is ap-
plied for deriving the spin density of an embedded radical. For the studied complexes formed by the Mgq cation and rare-gas atoms, the effect of the surrounding atoms on the radical’s spin density arise mainly from the Pauli repulsion between the overlapping electron densities. The supermolecule KS calculations by Eriksson w13x showed that the hyperfine structure of Mgq involved in complexes with rare-gas atoms can be well described at this level of theory.
2. Methods 2.1. The formalism The derivation of the KSCED equations in the spin-non-polarized case can be found elsewhere w5,6x. The spin-polarized version is outlined below. The spin density of the whole complex is represented as a sum:
r s s r 1s q r 2s ,
s s a or b ,
Ž 1.
where the densities r 1 and r 2 integrate to the total number of electrons N1 and N2 , respectively. Similarly to the KS theory, the electron density r 1s is represented by means of one-electron functions: N1 s < s <2 r 1s s Ý nŽ1.i cŽ1.i ,
s nŽ1.i s 0 or 1 .
Ž 2.
i
For a given r 2s, r 1s which minimizes the total energy, the functional Ew r a , r b x is derived from the following one-electron equations: s ,eff s y 12 = 2 q V KSCED Ž r. Ž r 1a , r 1b , r 2a , r 2b . cŽ1.i s s s ´Ž1.i cŽ1.i Ž r. ,
is1, N1 ,
Ž 3.
s,eff with the one-electron effective potential V KSCED taking the form: s ,eff V KSCED s
ZA
Ý y < ryR A
A
<
q
H
r 1 Ž rX . q r 2 Ž rX . d rX < rX y r <
q Vxcs Ž r 1a q r 2a , r 1b q r 2b . q
d Tsnadd r 1a , r 1b ; r 2a , r 2b dr 1s
,
Ž 4.
T.A. Wesołowskir Chemical Physics Letters 311 (1999) 87–92
where: Ø r i s r ia q r ib , i s 1,2; Ø Tsnadd w r 1a , r 1b ; r 2a , r 2b x s Ts w r 1a q r 2a , r 1b q r 2b x y Ts w r 1a , r 1b x y Ts w r 2a , r 2b x where Ts w r s x denotes the kinetic energy of non-interacting electrons in a reference system as defined in the KS theory w2x; Ø Vxcs is conventional exchange-correlation potential Vxcs s Ž d Exc . r Ž dr s . ; Ø the index A runs through the nuclei of subsystems 1 and 2. Eqs. Ž3. and Ž4. provide the spin-polarized version of equation Ž14. of Ref. w5x and equation Ž20. of Ref. w6x which where derived for overlapping closeshell systems. s,eff . The KSCED effective potential Ž V KSCED given by Eq. Ž4. can be formally partitioned into two components s ,eff s ,eff ,0 s ,eff ,int V KSCED s V KSCED q V KSCED ,
Ž 5.
where s,eff,0 Ø VKSCED does not depend on r 2 and it can be identified with the conventional KS effective potential of the isolated subsystem 1: s ,eff ,0 VKSCED s VKsS,eff Ž r 1 . s Ý y
A
H
q
ZA < r yRA <
r 1 Ž rX . d rX q Vxcs Ž r 1a , r 1b . ; < rX y r <
Ž 6.
s,eff,int Ø VKSCED collects all the terms arising from the interactions between subsystems:
s ,eff ,int VKSCED s
ZB
Ý y < ryR B
B
<
q
H
r 2 Ž rX . d rX < rX y r <
q Vxcs Ž r 1a q r 2a , r 1b q r 2b . yVxcs Ž r 1a , r 1b q
d Tsnadd r 1a , r 1b ; r 2a , r 2b dr 1s
89
For a given spin density Ž r 1s ., derived from the KSCED calculations the isotropic hyperfine coupling constant Ž A iso ŽMg.. is calculated as: 4p A iso Ž Mg . s gb g Mg b Mg² Zs :y1 3 = r 1a Ž r Mg . y r 1b Ž r Mg . , Ž 8. where Ø g and g Mg are electron and nuclear magnetons; Ø b and b Mg are electron and nuclear gyromagnetic ratios. 2.2. The embedding electron density: r 2 Since the electron density of the whole complex is represented using the electron densities of subsystems, it is possible to apply different levels of theory Žor computational effort. to derive each of them, depending on the a priori knowledge of the system under investigation. This offers the possibility of introducing additional approximations or simplifications. The embedding electron density r 2 is used for expressing analytically the embedding effective potential. Following Eq. Ž3., any change of r 2 involves a change of r 1. In this work, two ways for obtaining r 2 were applied. In the first, r 2 was taken from the KS calculations for the isolated embedding subsystem. In the second, it was obtained from the ‘freezeand-thaw’ cycle w14x in which r 1 and r 2 exchange their roles in subsequent KSCED calculations until they converge at values minimizing the total KSCED energy. The first way implies that the relaxation of the electron density of the embedding subsystems resulting from the interactions is neglected. In this work, it was assumed that r 2a s r 2b , i.e., the embedding electron density is not spin-polarized. 2.3. Localization of the electron density of subsystems
;
Ž 7.
Ø the indexes A and B run through the nuclei of subsystems 1 and 2, respectively. The terms in Eq. Ž7. represent the nuclear attraction, electron–electron repulsion, exchange-correlation, and the Pauli repulsion between the overlapping electron densities, respectively.
As has been discussed elsewhere w8x, for each subsystem the electron density Ž r i . can be localized near its nuclei or it can be allowed to spread over the whole complex. The former case offers substantial computational savings owing to the fact that a smaller number of basis set functions is needed to expand the electron density of the selected subsystem than is needed for the whole complex. In such a case,
90
T.A. Wesołowskir Chemical Physics Letters 311 (1999) 87–92
solving Eqs. Ž3. and Ž4. requires a computational effort comparable to that needed to obtain the spin density of the isolated Ži.e., not interacting with other molecules. subsystem. Since the basis sets used for expanding the fragments’ electron densities are not complete, the use of the localized basis sets is not well suited for describing such effects as the charge transfer between the subsystems or the formation of covalent bonds between them. To take into account such effects, the basis sets must be complete or at least include all functions describing the complex correctly by means of the supermolecular calculations. Such KSCED calculations involve, however, a computational effort exceeding that of the KS calculations for the whole complex. It is therefore not recommended for practical use. It is applied here to illustrate the effect that the incompleteness of the basis sets used to expand the fragment’s electron density has on the KSCED results. 2.4. Approximate functionals The PW86 w15x exchange and the P86 w16x correlation functionals were used to approximate the exchange-correlation functional in both the KS and in the KSCED calculations. Previous supermolecular KS studies of the same systems by Eriksson w13x showed that this approximation of the exchange-correlation functional is sufficiently accurate to describe the spin density. The non-additive kinetic energy part of the effective potential wŽ d Tsnadd w r 1a , r 1b ; r 2a , r 2b x.rŽ dr 1s .x was derived analytically w7x using the gradient-dependent kinetic energy functional. The used approximate kinetic energy functional was constructed following the Lee, Lee and Parr scheme w17x in which the same gradient-dependency as the one of the PW91 exchange functional w9x was used. 2.5. Numerical implementation The calculations were performed using a modified version of the deMon w18x program into which the KSCED formalism was implemented. A fine grid was used with 128 radial points in all calculations. The calculations were made using DZVP-type functions basis functions applied previously by Eriksson
in studies of the same system w13x. DZVP basis sets for neon Žcontraction pattern: Ž621r41r1) .. and TZVP for argon Žcontraction pattern: Ž73111r 6111r1) . were used in all calculations whereas three different basis sets were used for magnesium: Ø basis A ŽDZVP with contraction pattern 6321r 411)r1) ; Ø basis B ŽDZVP with decontracted valence orbitals and the contraction pattern 411111111r21111r1; Ø basis C Žfully decontracted DZVP basis sets: 111111111111r111111r1. The auxiliary basis functions for fitting the charge density and Vxc were derived from the original functions used by Eriksson w13x for the same systems Ž5,2;5,2. for Ne, and Ž5,4;5,4. for Mg and Ar. The modified functions, which can be labeled using the same convention as: Ž5,2;10,4. for Ne, and Ž5,4;10,8. for Mg and Ar, retained the original charge density part but the exchange-correlation part was supplemented by the functions fitting Ž d Tsnadd w r 1a , r 1b ; r 2a , r 2b x.rŽ dr 1s .. The functions used for fitting the kinetic energy part were obtained by a simple procedure in which all the coefficients of the original exchange-correlation functions were multiplied by two. The multiplication factor is justified by the r 2r3 dominant power dependence of the kinetic energy part of the effective potential. For each complex, the magnesium isotropic hyperfine coupling constants were calculated using a cluster comprising eight atoms to represent the matrix. The cluster atoms occupy the corners of the cube and the Mg–X ŽX s Ne or X s Ar. distance corresponds to the minimum at the supermolecular ˚ R MgAr KS potential energy surface Ž R MgNe s 2.97 A, ˚ s 3.39 A. w13x. 3. Results Table 1 shows the isotropic hyperfine coupling constants for Mgq embedded in the cluster of eight rare-gas atoms ŽNe or Ar.. The electron density Ž r 2 . taken from the KS calculations for the isolated X 8 ŽX s Ne or Ar. cluster was used to derive the ems,eff,int . bedding effective potential Ž V KSCED . The absolute values of the isotropic hyperfine coupling constant for free Mgq and embedded Mgq depend strongly on the basis set used for magnesium. The effect of the matrix on A iso ŽMg., as measured
T.A. Wesołowskir Chemical Physics Letters 311 (1999) 87–92
91
Table 1 Isotropic hyperfine coupling constants Ž A iso ŽMg. in gauss. of free Mgq ŽKS. and embedded Mgq ŽKSCED: non-relaxed electron density of the embedding atoms and the localized expansion of the electron density of subsystems.. The basis sets A, B, and C, correspond to different decontractions of the DZVP basis for magnesium Basis set
A iso ŽMgq .
A iso ŽMg8Neq .
A iso ŽMg8Arq .
A iso ŽMg8Arq . –A iso ŽMg8Neq .
A B C Exp w19x.
y197.3 y216.1 y209.1
y209.4 y229.9 y222.9 y222.4
y197.2 y216.8 y210.4 y211.6
12.2 13.1 12.5 10.8
by the difference between A iso ŽMg.Žembedded. and A iso ŽMg.Žfree., is not affected by decontracting the basis set. As a consequence, the difference between the isotropic hyperfine coupling constants of Mgq embedded in either a Ne or Ar matrix – Ne Ž Ž D AAr Mg. s A iso ŽMg8Arq–A iso ŽMg8Neq. does iso not vary with the change of the basis sets Ž12.2–13.1 gauss. and it agrees to within about 2 gauss with the experimental value Ž10.8 gauss w19x.. For the largest basis set Žfully decontracted DZVP., the absolute values of A iso ŽMg. are in excellent agreement the experimental measurements. 3.1. The effect of the relaxation of r 2 on the Mg q spin density The calculations described in the previous section do not take into account the modification of the electron density of the embedding molecules on the
calculated isotropic hyperfine coupling constant. To include such deformations, the spin density of Mgq was calculated using the r 2 obtained by means of the ‘freeze-and-thaw’ cycle of KSCED w14x. Such a cycle converges quickly for the Ne and Ar complexes studied. The values of A iso ŽMg. obtained in the third iteration differ by less than 0.1 gauss from the ones derived from the fully converged ‘freezeand-thaw’ cycle. The values of A iso ŽMg. derived using relaxed r 2 are shown in Table 2. It can be seen that the r 2 relaxation involves a very small Žusually - 1 gauss. change of the calculated values of A iso ŽMg.. 3.2. The effect of the charge transfer between the subsystems on the Mg q spin density Table 2 shows the isotropic hyperfine coupling constants derived using the supermolecular expansion of both r 1 and r 2 . It can be seen that allowing
Table 2 Isotropic hyperfine coupling constants Ž A iso in gauss. of embedded Mgq derived from the KSCED calculations using Ža. relaxed embedding electron density Ž r 2 ., and Žb. supermolecular expansion of the spin density of the subsystems. The basis sets A, B, and C correspond to different decontractions of the DZVP basis for magnesium. Values in parentheses denote the A iso changes with respect to the KSCED results derived using non-relaxed r 2 and localized expansion of r s1 and r s2 Basis set
A iso ŽMg8Neq .
Relaxation of r 2 : A y209.0 Žq0.4. B y229.5 Žq0.4. C y222.7 Žq0.2. Supermolecular expansion of r 1 and r 2 : A y207.2 Žq2.2. B y226.1 Žq3.8. C y219.2 Žq3.7. Exp w19x.
y222.4
A iso ŽMg8Arq .
A iso ŽMg8Arq . –A iso ŽMg8Neq .
y196.9 Žq0.3. y215.6 Žq1.2. y209.8 Žq0.6.
12.1 Žq0.1. 13.7 Žq0.8. 12.9 Žq0.4.
y197.0 Žq0.2. y214.9 Žq1.9. y207.7 Žq2.7.
10.2 Žy2.0. 11.2 Žy1.9. 11.5 Žy1.0.
y211.6
10.8
92
T.A. Wesołowskir Chemical Physics Letters 311 (1999) 87–92
the spin density to spread over the whole complex reduces the magnitude of A iso ŽMg. in all cases. The effect is larger for the neon complex than for that of – Ne Ž argon leading to the reduction of D AAr Mg.. iso
4. Conclusions The difference between the magnesium isotropic hyperfine coupling constants measured in Ne or Ar – Ne Ž matrices Ž D AAr Mg.. is well reproduced by iso means of the KSCED effective potential which is expressed as a function of the electron density of the atoms surrounding the radical. The KSCED calculations which use a non-relaxed electron density of embedding atoms and localized spin density of the radical lead to the values of the isotropic hyperfine coupling constants which agree to within 1 gauss with experiment for the largest basis set used. Such calculations involve a substantial reduction of computational effort compared to the supermolecular KS calculations. The absolute values of A iso ŽMg. depend strongly on the size of the basis orbital basis set used for magnesium. The relative value obtained as the difference between the A iso ŽMg. calculated for the complex and for the isolated Mgq which measures the effect of the matrix is rather basis set insensitive. – Ne Ž Similarly, the difference D AAr Mg. does not sigiso nificantly vary with the change of the magnesium basis set. The effects of the relaxation of the embedding electron density and the charge transfer between the interacting subsystems were also determined within the framework of the KSCED calculations. For the studied systems, they are not significant amounting to less than 20% of the calculated shifts of A iso ŽMg.. The results show that the Pauli repulsion between weakly overlapping electron densities can be accurately described by means of the functional derivative of the non-additive kinetic energy which is approximated using the gradient-dependent explicit energy functional. These results are in line with our previous analyses w8,7x of the applicability of the approximate kinetic energy functionals of the Lee, Lee, and Parr w17x form which indicated that for such functionals the PW91 enhancement factor is the best
approximation of the exact Tsnadd functional for small overlaps between electron densities of subsystems. Acknowledgements The author is grateful to Prof. J. Weber helpful for discussions, to Prof. D.R. Salahub for providing a copy of the deMon program, and to Prof. V. Malkin and Dr. O. Malkina for the subroutines for calculating the hyperfine tensor. Financial support from the Federal Office for Education and Science, acting as Swiss COST office, is gratefully acknowledged. This work is also a part of Project 20-49037.96 of the Swiss National Science Foundation. References w1x V. Barone, in: D. Chong ŽEd.., Recent Advances in Density Functional Theory, vol. I, World Scientific, Singapore, 1995, pp. 287–334. w2x W. Kohn, L.J. Sham, Phys. Rev. 140 A Ž1965. 1133. w3x V.G. Malkin, O.L. Malkina, L.A. Eriksson, D.R. Salahub, in: P. Politzer, J.M. Seminario ŽEds.., Theoretical and Computational Chemistry, vol. I, Density Functional Calculations, Elsevier, 1995. w4x L.A. Eriksson, in: M. Springborg ŽEd.., Density Functional Methods: Applications in Chemistry and Materials Science, Wiley, 1997. w5x P. Cortona, Phys. Rev. B 44 Ž1991. 8454. w6x T.A. Wesołowski, A. Warshel, J. Phys. Chem. 98 Ž1993. 5183. w7x T.A. Wesołowski, H. Chermette, J. Weber, J. Chem. Phys. 105 Ž1996. 9182. w8x T.A. Wesołowski, J. Chem. Phys. 106 Ž1997. 8516. w9x J.P. Perdew, Y. Wang, in: P. Ziesche, H. Eschrig ŽEds.., Electronic Structure of Solids’91, Academie, Berlin, 1991, p. 11. w10x T.A. Wesołowski, Y. Ellinger, J. Weber, J. Chem. Phys. 108 Ž1998. 6078. w11x N. Vulliermet, T.A. Wesołowski, J. Weber, Collabor. Czech. Acad. Sci. 63 Ž1998. 1447. w12x T.A. Wesołowski, N. Vulliermet, J. Weber, J. Mol. Struct. ŽTHEOCHEM. 458 Ž1999. 151. w13x L.A. Eriksson, J. Chem. Phys. 103 Ž1995. 1050. w14x T.A. Wesołowski, J. Weber, Chem. Phys. Lett. 248 Ž1996. 71. w15x J.P. Perdew, Y. Wang, Phys. Rev. B 33 Ž1986. 8800. w16x J.P. Perdew, Phys. Rev. B 33 Ž1986. 8822. w17x H. Lee, C. Lee, R.G. Parr, Phys. Rev. A 44 Ž1991. 768. w18x A. St-Amant, Ph.D. Thesis, Universite´ de Montreal, ´ 1992. w19x L.B. Knight Jr., C.B. Cleveland, R.F. Frey, E.R. Davidson, J. Chem. Phys. 100 Ž1994. 7867.