Valence electronic densities and electrostatic potentials in ab initio effective core potential procedures

Chemical Physics 147 ( 1990) 335-34 1 North-Holland

Valence electronic densities and electrostatic potentials in ab initio effective core potential procedures L. Fernandez Patios Departamento Quimica y Bioquimica, ETSI Mantes, Vniversidad Politecnica. 28040 Madrid, Spain

Received 26 March 1990

Over the past decade a variety of effective core potential (ECP) formulations have appeared in the literature so that a considerable number of valence basis sets and ECP operators are at present available. In this work it is demonstrated that despite their disparity, ECPs share common physical features related to the behavior of the electronic density previously found by Politaer on atoms. A connection in terms of electrostatic potentials for valenceonly electronic densities is also presented for valence energies computed by means of different ECP models.

1. Introduction

The reliability of ab initio effective core potentials (ECP ) in electronic structure calculations for atomic

and molecular systems is at present a well-established issue. A large body of evidence demonstrates that ECP procedures reproduce all-electron results to a high degree of accuracy [ 1,2 1. They are used not only by computational advantage reasons but, even more important for systems containing heavy atoms, because when formulated within a proper relativistic scheme, it is possible to develop relativistic effective potential (REP) operators which constitute a simple and efftcient way to introduce relativistic effects in molecular calculations [ 2,3]. Diverse extensive tabulations of ECP and REP op erators expressed in the form of analytical expansions, as well as optimized valence basis sets have ap peared in the literature over the past years (see for example refs. [ 4-61 and references therein). A considerable number of ECP formulations have been thus proposed, showing a similar performance in atomic and molecular calculations. This disparity of models can be contemplated as a consequence of the existing freedom to define pseudo-orbitals [ 7 1, which allows the generation of very distinct effective potentials. However, it must be pointed out that no connection with any physical meaning is usually considered apart from the reproduction of one or another type 0301-0104/90/$03.50

of all-electron (AE ) results, normally orbital energies and radial expectation values for valence electrons in atoms. Of course, a number of comparisons with experimental data in calculations including ECPs have also appeared in the literature, especially for systems containing transition elements [ 1,2]. However, such tests do not represent any direct validation of the effective potential used but, much more important, of the type of calculation carried out, i.e. kind of molecular wavefunction employed, degree of correlation treated, quality of the basis set, inclusion of relativistic effects, and so on. Since every effective potential is constructed from atomic wavefunctions, following a completely theoretical procedure based on an arbitrary partition of the electron space [ 7 1, there are no physical reasons to generate one or another type of ECP operators. It is the main purpose of this work to analyze the connection of ECP formulations with the core/valence separation issue as stated by Politzer et al. some years ago, from a physical point of view concerning electronic distributions (see chapters 1,2 of ref. [ 8 ] for a detailed discussion on this subject). This type of analysis seems pertinent for several reasons. Firstly, it provides a means to extract the physical significance common to any ECP procedure, as related to the core/valence partition in the atomic space. Secondly, it gives a precise idea of the degree of accuracy attained for properties not usually considered in ECP

0 1990 - Elsevier Science Publishers B.V. (North-Holland)

336

L. Fernandez Patios /Ab initio efective core potential procedures

calculations, as it is the case for electronic charge distributions. Thirdly, it establishes some grounds for finding electrostatic potential-valence energy relationships which can be further employed to obtain properties directly from electronic densities and/or electrostatic potentials in large molecular systems [ 81 containing heavy elements. If REP operators are considered, these relationships may serve to introduce relativistic effects in an approximate manner with no additional formal or computational difficulty. There are at least two reasons to explore these relationships: First, it is well known that quantum calculations give frequently more accurate one-electron properties (as electronic density, for example) than total energies. Second, and more important from a practical point of view, the electrostatic potential and the electronic density are real physical properties which can be determined experimentally, so that the theoretical information concerning them is directly veritiable in many situations as, for example, those related to molecular reactivity. 2. Core radii and valence charges Some years ago, Weinstein et al. [ 9 ] demonstrated that spherically symmetric all-electron Hartree-Fock charge densities for ground-state atoms can be represented in terms of monotonically decreasing functions of r. Whereas near the nuclei this radial dependence is sharply decreasing, it becomes less pronounced at a certain radial interval, although no single point can be associated with this change. There is however at least one minimum in the radial charge density R(r) = 4nrzp( r), located within the interval where the change on the radial dependence occurs. The position r, of this minimum is found to vary systematically along rows of the periodic table. Politzer and Parr [ lo] used this point to define the boundary limit between two distinct regions of the atom, and associated a physical meaning to this r value as the radius of the surface separating core and valence atomic regions. For atoms with principal quantum numbers greater than 2, r, is the outermost minimum in the total radial density function R(r). The physical significance of r, is highlighted by the fact that a procedure based on the Thomas-Fermi model, proposed by Politzer and Parr [ 10 ] to com-

pute the valence region energy, correctly predicts the experimental values, as commented below. If other distances from the nucleus are chosen to define r,, large discrepancies between valence energies and the corresponding ionization potentials are found, so that a clear uniqueness may be associated with the radius r,. After this introduction it seems pertinent to ask whether these core radii show any relationship with the radial distances implicitly taken as boundaries to separate core and valence regions in pseudopotentials of effective core potential formalisms. First of all, it must be remembered that the inherent nonuniqueness of pseudowavefunctions [ 71 allows certain arbitrariness in constructing ECPs. This fact is the cause of the considerable number of different ECP procedures so far available [ 1,2]. But with the only exception of model potentials as those by Huzinaga et al. [ $61, the majority of pseudo-orbitals have in common a smooth behavior and nodeless form. As a consequence, the core part of valence orbitals (pseudo-orbitals) keeps no resemblance with the AE counterpart. However, the pseudo-orbital must reproduce the behavior of the AE wavefunction in the outer atomic region. Most ECP procedures present a certain radial distance from which pseudo-orbitals and AE orbitals coincide. Let us call this distance rfit. It is obvious that the value of r,, will depend on the particular scheme followed to construct pseudo-orbitals. If we search any physical significance to be associated with this kind of core/valence partition, it must be expected that rfi, and r, values exhibit some relationship. As an example of ECP to carry out this analysis, we choose the averaged relativistic effective potential (AREP) [ 111 based on the shape-consistent technique proposed by Christiansen et al. [ 121 to define pseudo-orbitals. Details of the whole procedure are given in ref. [ 111 so that only the construction of pseudo-orbitals will be here recalled in order to get an idea of the meaning of r,, in that formalism. The shape-consistent technique was proposed with the aim to correct various errors encountered in earlier ECP formulations, whose origin was shown to lie in an improper partitioning of valence and core spaces [ 121. If the pseudo-orbital must represent the valence electron density accurately, one can define the radial dependence of the pseudo-orbital xv(r) as

L. Fernandez Patios/Ab initioeffectivecore potentialprocedures

x”(r)=v”/v(r)+F”(r)

(1)

7

where vV( r) is the AE radial valence orbital and F,(r) is a function chosen to cancel the radial oscillations. in the core region eliminating thus nodes and going smoothly to zero in the valence region. Fv( r) is determined by matching a seven-term polynomial to wV(r) at a certain point which is precisely I+~,.F,(r) is then taken as the difference between the polynomial and w”(r) for r values lower than rfit and it is set equal to zero elsewhere [ 111. As a consequence of this technical procedure to define pseudo-orbitals, rfit represents the innermost point at which a smooth nodeless function can be obtained to construct xv(r) . The valence orbital vV( r ) is thus partitioned into a pseudo-orbital and the corresponding core segment, -F,(r), and the boundary limit between these two parts is r,,. Although it is evident that this radial distance is a particular parameter within the procedure to construct a particular type of ECP, it exists a remarkably accurate linear dependence between the core radii r, by Politzer (tabulated by Boyd [ 13 ] ) and the pseudo-orbitals rfitradii. This relationship is shown in fig. 1 for atoms of the two first rows of the periodic table. The correlation coefficients for first and second row are 0.9974 and 0.9980 respectively. If one considers that the physical and

3.0 -

2.5 3

.3

2.0 -

5 1.5 -

1.0 -

0.5 -

'%.b

potentials at

nuclei

3.5 -

c

formal models underlying r, and rfi, quantities are rather different, this linear relationship is striking and suggests a common information encoded in both definitions of the core regions in atoms. It is interesting to continue this analysis focusing on the core/valence separation of electronic charges. Table 1 gives the valence charges for first- and second-row atoms integrated over the outer atomic region defined from r,. These charges have been computed by using a number of representative ECP models with different basis sets and effective potential operators implied. Only the model potential method (column labelled MP) gives wavefunctions with the correct AE nodal structure while tlie others yield pseudo-orbitals showing no resemblance with AE orbitals in the core region. Despite the disparity inherent to all the valence-only ECP calculations, valence charges in table 1 are remarkably consistent. AE values computed by Politzer and co-workers from atomic Clementi-Roetti wavefunctions include the contributions from core orbitals to the total charge in r,-m and they are thus slightly larger. However, there is again an excellent linear correlation between AE HF and AREP values, being these a representative choice among ECP models. In fact, separate fits for charges in first- and second-row atoms give correlation coefficients 0.999997 for both cases.

3. Valence energies and electrust&

4.0

337

0.;

1.b

1.;

2.b

2.;

3.b

3.L

r, (a.u.) Fig. 1. Plot of radial boundaries between core and valence parts of pseudo-orbitals in the AREP procedure, r,, [ I 11, versus outermost minima in radial charge densities, r, [ 13 1,for first- and second-row atoms.

Assuming the R(r)minimum, r,, defines a core region in such a form that the electrons outside can be treated as moving in the effective resulting field, Politzer and co-workers derived within the framework of the Thomas-Fermi theory a formula to determine the valence region energy [ 8,10,16 1. The energies of outer electrons thus computed are remarkably close to the respective sums of ionization potentials only if the r, value is used, whereas large discrepancies are found for a number of other radial distances [ 16 1. To make this comparison, Politzer took the sum of the first Z, ionization potentials, Z,: being the valence charge (last column in table 1). Since Z, is not in general an integer number, only the appropriate fraction of the last ionization potential, determined by linear interpolation, was included in the sum. This

L. Fernandez Patios /Ab initio effitive corepotential procedures

338

Table 1 Valence charge integrated in r,-ao for first- and second-row atoms ‘) Atom

Li Be B C N 0 F Ne Na Ma Al Si P S Cl Ar

Valence-only ECP calculations MP b’

AREF =)

PWG *’

SBK ”

a.944 1.912 2.860 3.804 4.757 5.711 6.669 7.644 0.782 1.650 2.626 3.590 4.555 5.521 6.470 7.426

0.943 1.914 2.867 3.817 4.772 5.728 6.689 7.666 0.779 1.644 2.621 3.587 4.553 5.525 6.485 7.449

0.942 1.913 2.866 3.820 4.768 5.723 6.682

0.946 1.918 2.868 3.817 4.772 5.732 6.693 7.672 0.781 1.646 2.623 3.593 4.566 5.541 6.501 7.467

0.696 1.597 2.606 3.572 4.563 5.526 6.506

“r,valuesaregiveninref. 113). b, Model potential by Huzinaaa et al. [ 61. ‘) AREP calculations. [ 111. d)ECPbyPettersson,Wah&penandGropen (41. ‘)ECPbyStevens,BaschandKrauss I) HF functions by Clementi and Roetti. Values taken from ref. [ 131.

comparison was used by Politzer to highlight the physical meaning of rc. A similar analysis for different ECP procedures is shown in table 2, where valence energies are compared to the sum of the first NVexperimental ionization potentials, NVbeing now the number (integer) of valence electrons. Since ECP calculations refer to the atom resulting after formally removing core electrons (pseudoatom), these values are total energies for this particular problem and correspond to the basis sets given in every model (see footnotes to table 2). Although it should be possible to obtain lower energies by using larger basis sets, values in table 2 are given to demonstrate the good agreement with experimental ionization potentials for valence electrons. Of course, this comparison is not intended for giving a precise idea of the quality of ECP calculations. (SCF in all cases) but for showing that the core/ valence separation inherent to any ECP procedure implicitly defines a physically meaningful valence region in a parallel manner to that given in the Politzer study. There is a complementary point of view on the core/valence separation which concerns electrostatic

[IS].

WHe’

AEHFs’

0.782 1.649 2.627 3.592 4.565 5.541 6.487 7.458

0.956 1.946 2.908 3.869 4.837 5.815 6.800 7.793 0.785 1.675 2.671 3.675 4.674 5.668 6.666 7.669

“ECPbyWadtandHay

[15].

potentials. It is well known that energies of atoms and molecules and electrostatic potentials at their nuclei are related [ 8 1. Since the electrostatic potential is a real physical property obtainable from the electronic density function, considerable interest has been devoted to explore the connections between these two quantities. Let us recall very briefly some points relevant to our purposes. The electrostatic potential at the point r produced by the electrons and nuclei of a system having an electronic density function p(r) is given in atomic units by

where 2, is the charge on nucleus A located at RA. The electrostatic potential V(r) is also rigorously related to the electronic density by Poisson’s equation vW(r)=47r&r). If the system is a single atom, the electrostatic potential at the nucleus is then if,=-

Ip(r) r

dr

.

L. Fernandes Patios /Ab initio effdive corepotentid procedures

339

Table 2 Total valence energy (hat-tree) for tint and second-row atoms Atom

Valence-only

Li Be B C N 0 F Ne Na MS Al Si P s Cl AT

ECP calculations

MP ”

AREP b,

SBK c,

-0.1955 -0.9587 -2.5369 -5.3183 -9.6519 - 15.6804 -23.8994 - 34.6553 -0.1814 -0.7826 - 1.8743 -3.6705 -6.3454 -9.9366 - 14.7158 -20.8482

-0.1964 -0.9619 -2.5337 - 5.3044 - 9.6222 - 15.6412 -23.8362 - 34.5699 -0.1823 -0.7862 - 1.8770 -3.6703 -6.3387 -9.9261 - 14.6952 -20.7931

-0.1960 -0.9602 -2.5341 - 5.3027 -9.6177 - 15.6224 -23.7931 - 34.4954 -0.1816 -0.7836 - 1.8739 - 3.6645 -6.3264 - 9.9029 - 14.6495 -20.7373

WHd’

Exp. ,Y IP =)

-0.1806 -0.7810 - I .8695 - 3.6757 -6.3152 -9.8738 - 14.6808 -20.6732

0.1981 1.0118 2.6232 5.4397 9.8099 15.9157 24.21 IO 35.0433 0.1889 0.8335 1.9573 3.7899 6.4966 10.1629 15.0256 21.2314

‘) Model potential by Huzinaga et al. [ 61. ‘) AREP calculations [ 1 I 1. C~ECPbyStevens,BaschandKrauss[14]. “)ECPbyWadtandHay[lS]. e, Sum of the first N, experimental ionization potentials from ref. [ 17 1.

For an atom with N electrons and nuclear charge 2, the nonrelativistic total energy is related to the electrostatic potential at its nucleus, PO, by the exact expression [ 18 ]

B(Z r)=rUZ

r)lZ,

(8)

eq. (7) can be written as (9)

Although this relation is rigorous, the evaluation of its integral is so difficult that it has a scarce practical application in exact treatments. There are however approximations to obtain simpler forms of this equation. Let us first assume spherical symmetry for p(r). Expression (4) is then al v, = -4x

I 0

rp( Z, r) dr ,

(6)

If we define ~6 = (+/at),=,, a relationship between the screening factor and V. is found:

0!!f z,r=O=Zvo4. ar

(10)

Therefore

avo az=Qb+z$g,

(11)

which, substituted into ( 5) gives an approximate expression for the energy in terms of the function pb

(12)

and the Poisson equation has the form 0

V2V(z, t) =($3+

;(3z=4V(Z,

r) .

(7)

If a function bp(Z, r) (usually known as screening factor) is defined as

This equation, derived by Politzer and Parr [ 18 1, is much easier to apply than eq. (S), although it requires a particular form for the screening factor, or equivalently for the electrostatic potential. This point

340

L. Fmandtz

Pacm /Ah uwo C~~~IIYPcow porenrralpowdurrs

is again related to the core/valence separation. As stated above, there is considerable evidence indicating that core and valence parts of electronic densities have different behavior. As a consequence, their relationships with electrostatic potentials are also distinct. Thus, while a Thomas-Fermi type equation I’( Z, r) =/?[p( Z, r) 1”’ provides a good representation for the V-p relation in the core region [ 16 1, in the valence region it is necessary to consider an expression of the type I’(Z!,r)=a(Z)p(%,r).

,

-D$Yr)

where D and u arc parameters to be determined. expression implies an electrostatic potential L’(Z,.r)=

$

exp(

and an electronic P(%,,r)=

-D%",r).

(Ida) This

(14b)

density

$Fexp(_DZ:r),

so that a(Z) = (DZ:)‘/4n in ( 13). These formulae represent improved forms of expressions introduced by Politzer [ I6] to obtain valence energies for atomic outer regions defined from r,. With eq. ( l4c) it is easy to verify ,I)

4x

r2p( Z,. r)

row

atoms

(13)

A particularly simple expression showing the correct behavior for valence electronic densities can be encountered from a screening factor defined as (p(Z,. r)=exp(

Second

dr= Z,

.

1 0

The definition ( 14a ) leads 10 dp; = - DZ: and therefore the electrostatic potential at the nucleus is V0= -DZ:+'.From this & expression, the energy ( 12) reduces 10 the simple equation

-.l

sl,

0;

1

.o In

1.5

20

23

2,

FIN. 2. Logarithmtc plot of toul vrlcncc encrgxs computed oath AREP opcra~ors.E [ I I 1, versus valence chraa. Z. for lirsl- and second-row atoms.

way 10 that employed by Politzer for AE electronic densities. It should be therefore expected that when applied to energies computed by means of any ECP procedure, expression ( 15) holds for a given row of the periodic table, i.e. parameters D.u must be associated to a given core for atoms with valence charge Z,. If this is the case, the plot of In( -E) versus In Z, must be a straight line. Fig. 2 gives this plot for AREP energies [ I I 1. The correlation coefficient for tint-row atoms is 0.9993 and for second row 0.999 I. Since total energies for different ECP models are quite similar, as values in table 2 demonstrate. this same good linear correlation can be expected for any effective potential formulation. In fact, the extension of this analysis to all the ECP models considered above, leads 10 linear correlations with coefficients invariably better than 0.999.

4. coachsims E=-

- D z:+2. a+2

(15)

It is inmediate to demonstrate that Hellmann-Feynman theorem, iK/aZ= V, is fulfilled if Z= Z,. The derivation of formula ( 15 ) has been made under the assumption that electronic densities in the valence region exhibit a behavior which can be described by an expression of the type ( 14c). in a similar

The analysis presented in this work demonstrates that although effective core potential formulations are primarily devised as mere computational tools, they enclose physical information which can be traced in terms of the behavior ofelectronic densities. The entity resulting aAer the core/valence separation in atomic wavefunctions is made to construct ECP op

L. Fernandez Patios /Ab initio effective core potential procedures

erators (pseudoatom), not only represents an adequate simpler system to carry out molecular treatments, but keeps the essential information concerning electronic density properties. If one considers that any ECP procedure implies several modifications of atomic treatments, leading to functions (pseudo-orbitals) which show in general a different behavior with respect to all-electron orbitals, that is not a priori evident. Moreover, since there is a variety of ECP formulations, it is a meaningful finding that the very distinct models analyzed share a common electronic density behavior. Despite the absence of core electrons, this behavior is completely analogous to that encountered by Politzer for all-electron atomic wavefunctions. If one thinks that a number of molecular properties are directly related to the electronic density, or equivalently to the electrostatic potential, the search for reliable procedures to determine these properties without necessarily performing complex quantum calculations, presents an evident interest. The experimental studies concerning molecular reactivity constitute a field where such theoretical treatments should represent valuable information.

341

References [ 1] M. Krauss and W.J. Stevens, Ann. Rev. Phys. Chem. 35 (1984) 357.

[ 2) P.A. Christiansen, W.C. Ermler and KS. Pitzer, Ann. Rev. Phys. Chem. 36 (1985) 407.

[ 3 ] K. Balasubramanian and KS. Pitzer, Advan. Chem. Phys. 67 (1987) 287. [4] L.G.M. Pettersson, U. Wahlgren and 0. Gropen, J. Chem. Phys. 86 ( 1987) 2 176. [ 5) V. Bonifacic and S. Huzinaga, J. Chem. Phys. 60 (1974) 2779. [ 6 ] S. Huzinaga, L. Seijo, Z. Barandiaran and M. Klobukowski, J. Chem. Phys. 86 (1987) 2132. [ 71 L. Szasz, Pseudopotential Theory of Atoms and Molecules (Wiley, New York, 1985). [ 81 P. Politzer and D.G. Truhlar, eds., Chemical Applications of Atomic and Molecular Electrostatic Potentials (Plenum Press, New York, I98 1). [9] H. Weinstein, P. Politzer and S. Srebenik, Theoret. Chim. Acta 38 (1975) 159. [lo] P. Politzer and R.G. Parr, J. Chem. Phys. 64 (1976) 4634. [ 111 L. Femandez Patios and P.A. Christiansen, J. Chem. Phys. 82 (1985) 2664. [ 121 P.A. Christiansen, Y.S. Lee and K.S. Pitzer, J. Chem. Phys. 71 (1979) 4445. [ 131 R.J. Boyd, J. Chem. Phys. 66 (1977) 356. [ 141 W.J. Stevens, H. Basch and M. Krauss, J. Chem. Phys. 81 (1984) 6026. [ 151 W.R. Wadt and P.J. Hay, J. Chem. Phys. 82 (1985) 284. [ 161 P. Politzer, J. Chem. Phys. 72 (1980) 3027. [ 171 C.E. Moore, Natl. Stand. Ref. Data Ser. Natl. Bur. Stand. 34 (1970). [ 181 P. Politzer and R.G. Parr, J. Chem. Phys. 61 (1974) 4258.