J . Phys. Chem. 1991, 95, 10653-10658
10653
A Simplified Representation of Atomic Electron Densities and Electrostatic Potentials L. Fernandez Pacios Departamento Quimica y Bioquimica, E TSI Montes, Uniuersidad Politecnica, 28040 Madrid, Spain (Received: May 31, 1991)
A procedure to represent atomic electron charge densities p ( r ) and electrostatic potentials V(r) by means of simple functions is presented. The aim is to construct fast to compute analytical functions for the systematic calculation of p(r) and V(r) maps in systems containing large molecules. Taking as starting point the optimized potential model for obtaining accurate atomic local potentials, the procedure proposed allows an overall representation of these properties in a simplified yet reliable manner. The parameters obtained for the first-row atoms Li-Ne are given and the quality of such an analytical representation is discussed.
Introduction There is an obvious interest in obtaining accurate electron charge densities for studying the structure and properties of biomolecular systems. The considerable amount of data directly related to the electron density is not the only reason for this interest.’ In addition, an increasing number of applications of density functional formalisms are currently attracting great attentiom2 Moreover, the elegant theory by Bader on the old ideas of “atoms in molecules”, incorporating topological features to the study of molecular systems is built upon the electron density and its L a p l a ~ i a n . ~ Other important property directly obtainable from the electron density, the molecular electrostatic potential, is actually one of the most useful concepts in many fields of theoretical chemistry? Extensive applications ranging from calculation methods to determine molecular energies to studies of enzymatic activity have been devised for electrostatic potentials. Since it is a real physical property, studies centered on it provide a more physically meaningful description than methods based on atomic charges, necessarily defined in a more or less arbitrary manner. Most of the above-mentioned applications require, however, electron density or electrostatic potential maps, i.e., evaluation of these properties for a mesh formed by a large number of points. Therefore, an efficient method for the computation of these maps happens to be necessary. Accurate electron densities may be easily obtained from ab initio wave functions for small molecules but, since the computation of electrostatic potentials or electron densities is a problem which depends on the number of electrons, this task turns formidable for a large molecule with hundreds or even thousands of electrons. In order to systematically compute these maps for biomolecules, a simple yet reliable method is clearly needed. An optimum equilibrium between reliability and simplicity is of course not easy. To devise such methods an appropriate starting point should be establishing atomic rather than electronic functions to represent electron densities and electrostatic potentials. Along these lines is inscribed this work. The main goal is to develop a method to reliably encode the information of accurate all-electron atomic wave functions into analytically simple functions. Apart from the simplification represented by a computational procedure in which the building blocks are atomic functions instead of electronic ones, a supplementary advantage over conventional molecular semiempirical formulations is that one has an explicit analytical form for every atomic contribution to the total property being considered. Moreover, as will be shown (1) Fliszar, S . Charge Distributions and Chemical Effects; Springer: New York, 1983. (2) Kryachko, E. S.; Ludefia, E. V. Energy Density Functional Theory of Many-Electron Systems; Kluwer: Dordrecht, 1990. (3) Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Oxford University Press: Oxford, UK, 1990. (4) Politzer, P., Truhlar, D. G., Eds., Chemical Applications of Atomic and Molecular Electrostatic Potentials; Plenum: New York, 1981. Buckingham, A. D.; Fowler, P. W.; Hutson, J. M. Chem. Rev. 1988, 88, 963.
0022-3654/91/2095-10653$02.50/0
in a forthcoming paper, it is also possible to construct atomic functions for the representation of exchange and correlation effects within the same model. This paper is organized as follows. In next section, the atomic procedure used as starting point (optimized potential model by Talman et aL5) is briefly considered along with the initial considerations on the different behavior exhibited by the atomic electron density. The method for constructing the above-mentioned functions is presented in the third section with application to first-row elements and, finally, a summary and main conclusions are included.
Optimized Potential Model Calculations and Density Partitioning The most widespread procedure to determine atomic wave functions is the Hartree-Fock (HF) method in which the electrons are considered as moving in an average effective potential originated by the nucleus and the electrons. The wave function is obtained by solving a set of integrodifferential equations involving the nonlocal exchange potential, therefore HF wave functions can not be determined from a local potential formulation. If one desires to find an atomic effective local potential for the electrons, other procedures have to be considered. One of the oldest and best known is the Thomas-Fermi statistical model which may be used only for semiquantitative purposes due to its rather approximate nature. Also developed long time ago is the Slater Xa model’ in which a semiempirical relation is established between the exchange potential and the electron charge density. A more recent approach and one almost coincident with HF results is the optimized potential model (OPM) by Talman et al.s In this procedure an effective local potential, variationally optimized, is obtained through next scheme. Consider a central potential U(r),solve until self-consistency the expressions arising from the Schrodinger equation for the single-particle orbitals, calculate ( H ) ,and then vary U(r)to minimize (H). The exposition of the OPM method as well as a discussion of results can be found in refs 5 and 6. The main reason for choosing this model as starting point for our work is that it provides an effective local potential along with an exchange potential while keeping the degree of accuracy at the HF level.6 Atomic orbitals, electron density, and both potentials are given in numerical form for a radial logarithmic mesh uniformly distributed in a certain variable related to r so that points are closely spaced near the nucleus and uniformly spaced for large r val~es.~ To derive a simplified representation for electron densities, it is useful to start at the virial theorem V, = -2T, where V, is the total potential energy and T the kinetic energy. Since E = T
+
VI
9
E = Y2V1 ( 5 ) Talman, J. D.; Shadwick, W. F. Phys. Rev. A 1976, 14, 36. Talman, J. D. Comput. Phys. Commun. 1989,54, 85. (6) Aashamar, K.; Luke, T. M.; Talman, J. D. Phys. Reu. A 1979, 19,6.
0 1991 American Chemical Society
10654 The Journal of Physical Chemistry, Vol. 95, No. 26, 1991
Fernandez Pacios
+ V , ) Potential Energies, and Total Energy ( E ) for First-Row
TABLE I: Kinetic Energy (T), Nudea~Electronic( VNJ and Interelectronic ( V , Atoms Li-Ne"
T atom Li Be
OPMb 7.431 56 14.5706 24.5248 37.6538 54.2882 74.7584 99.3947 128.528
B C N
0 F Ne
" All energies in hartrees.
Vec + OPMb 2.281 40 4.489 24 7.837 94 12.7440 19.5086 28.4383 39.8414 54.0264
-VNe
H Fc 7.432 75 14.5730 24.5290 37.6883 54.4005 74.8093 99.4091 128.547
OPMb 17.1453 33.6322 56.8905 88.0559 128.091 177.964 238.644 311.100
H F' 17.1467 33.6348 56.8972 88.1360 128.352 178.076 238.667 311.132
Optimized potential model calculations, this work.
atom Li Be B
C N
0 F Ne
NC
N"
rmin
Nc
1.763 1.096 0.800 0.622 0.503 0.419 0.356 0.307
2.044 2.053 2.09 1 2.130 2.161 2.185 2.199 2.206
0.956 1.947 2.909 3.870 4.839 5.815 6.801 7.794
1.764 1.096 0.800 0.621 0.502 0.419 0.356 0.307
2.044 2.054 2.092 2.131 2.163 2.185 2.200 2.207
Nv 0.956 1.946 2.908 3.869 4.837 5.815 6.800 7.793
4.0
Ormi,,(bohr) is the position of minimum in the radial distribution function 4r&(r). Optimized potential model calculations, this work. Hartree-Fock analytic values from ref 11.
If correlation effects are ignored, the potential can be expressed as Thus written, V,, is the nuclear-electrons attraction potential energy, V, the interelectronic repulsive interaction, and V, the exchange contribution to the electron-electron potential energy. Table I shows T,VN,, V , + V,, and E energies for first-row elements. OPM results are compared with the analytic, essentially exact, Hartree-Fock values by Clementi-Roetti.8 Although slight differences are found in some values, especially for N and 0, data given in this table illustrate the reliability of this method. One must keep in mind that the energies chosen for this comparison are near-HF limit values. But our attention will be mainly focused on electron densities and electrostatic potentials; these properties as determined by means of OPM calculations are now considered. It is known that the electron density p(r) is a function which decreases monotonically from its maximum at the nucleus p(r=O) = po, changing its rate of exponential decrease at certain radial interval^.^ This behavior defines as many different exponential regions as there are n quantum numbers. For each atom the transition from one region to other takes place over a radial interval, so that it cannot be associated to any single point. However, Politzer and Parrlo have shown that, within this transition interval, there is a well-defined point, corresponding to the minimum in the radial distribution function D(r) = 4 d p ( r ) , r = rmin.These authors have proposed t h a t this point defines a boundary surface separating core and valence regions in atoms. In case of atoms with principal quantum numbers greater than 2, this definition is applied to the outermost minimum in D(r). Boyd" has determined rminvalues for all the atoms Li-Xe, using Clementi-Roetti H F wave functions: as well as charges, invariably (7) Slater, J. C. The Selfconsistent Field for Molecules and Solids: Quantum Theory of Molecules and Solids; McGraw Hill: New York, 1974; Vol. 4. Slater, J. C. The Calculation of Molecular Orbitals; Wiley: New
York, 1979. (8) Clementi, E.; Roetti, C. A t . Data Nucl. Data Tables 1974, 14, 177. (9) Weinstein, H.; Politzer, P.; Srebenik, S. Theor. Chim. Acta 1975, 38, 159. (10) Politzer, P.; Parr, R. G. J . Chem. Phys. 1976, 64, 4634. (11) Boyd, R. J. J . Chem. Phys. 1977, 66, 356.
OPMb 7.432 35 14.5724 24.5278 37.6581 54.2946 74.7677 99.4079 128.546
H F' 7.432 73. 14.5730 24.5291 37.6886 54.4009 74.8094 99.4093 128.547 ~~
h OPM Electron Density
H F'
rmin
-E
H Fc 2.281 18 4.488 76 7.839 24 12.7591 19.5507 28.4576 39.8490 54.0379
Hartree-Fock analytic calculations from ref 8.
TABLE 11: Electron Density Partitioning into Core (O-rmin)Charge, N, and Valence (r--m) Charge, N,, for First-Row Atoms" OPMb
vx
- - - - - - - : Be
1.o
2.0
r,,,h = 1.096
3.0
r (bohr)
F i 1. Logarithm of the electron density p(r) from OPM calculations for beryllium and neon. Full squares indicate the position of r,,, distance of minimum in the radial distribution 4 d p ( r ) . TABLE III: Values of Electron Densities, p,,, and Electrostatic Potentials, V , (atomic units), at the Nucleus for First-Row Atoms Li-Ne
VO
PO
atom Li Be B C N 0
F Ne
OPM" 13.8140 35.3874 71.9343 127.522 206.150 311.815 448.524 620.281
HFb 13.8342 35.4277 71.9846 127.555 206.135 311.975 448.610 620.146
OPM" -5.71536 -8.408 19 -11.3785 -14.6766 -18.2999 -22.2472 -26.5183 -31.1130
"Optimized potential model calculations, this work. analytic calculations from ref 8.
HFb -5.71555 -8.40870 -11.3794 -14.6893 -18.3360 -22.2595 -26.5186 -31.1132 Hartree-Fock
noninteger numbers, integrated in the different regions defined by the minima in D(r). It has been recently shownI2that the point r,, has also a physical meaning in effective core potential (ECP) f0rma1isms.l~ This finding gives a more remarkable significance to rminif one considers that ECP formulations are developed without any a priori defined radial distance to separate core and valence regions. Table I1 gives HF and OPM rminvalues for Li-Ne atoms and the corresponding charges integrated in O-rmin(core charge) and rmi,-m (valence charge). Since OPM electron densities are obtained in numerical form, the position of the minimum in D ( r ) has been determined by means of a four-point Lagrange interpolating po1yn0mial.l~ The charges have been calculated by integrating numerically the function 4 d p ( r ) in the proper in(12) Fernandez Pacios, L. Chem. Phys. 1990, 147, 335. (13) Fernandez Pacios, L.; Christiansen, P. A. J . Chem. Phys. 1985, 82, 2664. For a general exposition see: Szasz, L. Pseudopotential Theory of Atoms and Molecules; Wiley: New York, 1985. (14) Press, W. H.; Flannery, B. P.; Teukolsky, S.A.; Vetterling, W. T. Numerical Recipes. The Art of Scientific Computing, Cambridge University Press: Cambridge, UK, 1989; Chapter 3.
Electron Densities and Electrostatic Potentials
The Journal of Physical Chemistry, Vol. 95, No. 26, 1991
tervals. Since the independent variable is not equally spaced in the radial mesh, parabolas have been fitted to data points, three at a time, but using the parabola to represent the integrand only in the first panel, except for the last three points, for which the parabola is used for both panels. This technique is known to be equivalent to the usual Simpson rule. As is clearly apparent from the data in Table 11, there is almost no difference between OPM and H F integrated charges and rminvalues. Figure 1 shows the behavior mentioned before for the electron density. The logarithm of p(r) versus r for beryllium and neon is plotted in order to illustrate the change of the exponential decrease rate. The positions of rmi,are also indicated. A final test on the reliability of OPM calculations oriented toward our purposes is given in Table 111. Values of electron densities and electrostatic potentials at the nucleus are shown for Li-Ne atoms, compared with the corresponding H F results. The atomic electrostatic potential V(r) = - r
ss?
PV’) dr“
17 - 7’1
(3)
10655
If r- is taken as the boundary radius between core and valence regions, the partitioning of density (6) into the corresponding charges, gives for the inner (core) part
N, = 47r[Acsr* 0 ? exp(-B,r) dr
+ A, Srhr2 exp(-B,r) 0
dr]
After substituting I? exp(Br) dr = (eBr/B3)[(l - Br)2 + 11 and a bit of manipulation, the core charge is 41rA, N , = -[2 - [(l + Bcrmin)2+ 11 exp(-B,rmi,)]
+
8 :
In exactly the same manner one gets for the valence charge
Nv = 4 ~ 1 ? : p(r) dr the expression
takes at the nucleus the value V, = -
d7
(4)
+
Therefore, it can be calculated as
Adding eqs 10 and 11 it is straightforward to confirm that N, N , = Z with Z given by eq 8. The minimum condition in D(r) = 47dp(r)
where the sum runs over occupied atomic orbitals with occupation numbers ni. HartreeFock Vo values in Table I11 have been determined in this form. In OPM calculations one gets the electrostatic potential and the electron density in numerical form at every point of the radial mesh utilized. Voand po in this case have been obtained by extrapolating to r = 0. Since extrapolation is a generally more risky operation than interpolation, two different methods have been employed: the usual Lagrange polynomial and the rational function extrapolation (Bulirsch-Stoer algorithm).I4 In both cases four points have been chosen and comparable estimated errors are found, of the order of 10-4-10-3 in po and 10-5-10-4 in Vofor Li-Ne atoms. Final extrapolated values presenting the least estimated error are then selected for po and
for the model density (6), leads to the relation Ac(2 - & m d
exP(-Bcrmin) + A d 2 - Bvrmin) exP(-Bvrmin) = 0 (12) which allows one to obtain r,, from the set of parameters A,, B,, A,, B,. Let us now turn to the determination of the electrostatic potential arising from the electron density (6). Both properties are related by means of the Poisson’s equation OzV(7)= 4rp(7),which under spherical symmetry may be written in the form
VO. Development of Atomic Functions for p ( r ) and V ( r ) The evidence presented above lends support to the choice of the OPM method as a ground to develop approximate functions for representing total electron densities p(r) and electrostatic potentials V(r). At this first stage, the discussion is limited to first-row atoms Li-Ne. The before mentioned behavior of p(r), illustrated in Figure 1, allows one to establish two radial regions in p(r): core (O-rmin) and valence (rmin-=). Let us therefore assume that the total p(r) can be represented by means of a sum of two exponential functions p(r) = A, exp(-B,r) + A, exp(-Bvr) (6) where subscripts in parameters A and B refer to the core and valence prts. The normalization condition imposed to this density is, introducing spherical symmetry Z=
s
a*V(r) + -2-av(r) a? r ar
- 47WM
In order to find V(r) from the known p(r) let us define a ‘‘SLreening functionn4J04(r) = rV(r)/Z. This way, V(r)= Z4(r)/r and the Poisson’s equation takes the form d24(r) - - y47rr P ( r ) dr2
which, after substituting ( 6 ) and solving for 4(r), gives r
.
The electrostatic potential is therefore r
pd7=
4 r [ A , X m r 2exp(-B,r) dr
+ A , I m r 2exp(-B,r)
dr] (7)
Recalling that SrI“ exp(-Br) dr = n!/Bn+I,B > 0 and n integer positive, the condition (7) is 87rA, 8rAv
z=-
Bc3
+- Bv3
[ F]
=
v,/z
r=O
which, after using (14), is found to be r
The density at the nucleus is simply Po = A , + 4
The values of V(r) at the nucleus, Vo,is found from the relationlo
(9)
.
Fernandez Pacios
10656 The Journal of Physical Chemistry, Vol. 95, No. 26, 1991
At this point we have already obtained the functional forms for p ( r ) and V(r). It is then possible to set conditions related to atomic energies for determining the parameters needed. If the virial theorem (1) is expressed in terms of potential energy contributions (2), one gets E =
1/2(VNe
+ Vee + V x )
(17)
where correlation effects have been ignored. The nuclear-electrons attraction energy
TABLE I V Parameters in the Model Electron Density (6) for First-Row Atoms Li-Ne
atom Li Be B
C N
0
F Ne
BC 5.818505 7.888 538 9.983 889 12.103679 14.242 241 16.397 720 18.569401 20.756 924
B"
A" 0.045 482 0.141 982 0.410 8 13 1.005 507 2.112260 3.959 890 6.818641 10.999 834
0.972 372 1.175 742 1.478 643 1.817853 2.165885 2.516477 2.867 754 3.219 066
TABLE V Core (N,) and Valence (N,) Charges, Radii rmin (bohr), and Electrostatic Potentials at the Nucleus, V o (atomic units), for the Model Ekctron Density
for the density (6) is easily obtained: VN, = -4nZ
4 13.768 488 35.245458 71.523467 126.51682 204.03720 307.855 42 441.705 10 609.28075
1 - r p(r) d r = -4nZ
atom Li Be B C N 0 F Ne
0
Note looking at (16) and (19) that the relation VN,= ZVo is in fact fulfilled. Let us consider now the second term in the energy expression (17), the interelectronic potential energy
If the electrostatic potential, as defined by eq 3, is multiplied by p ( 7 ) and integrated over 7, the result
NC 2.059 2.097 2.154 2.203 2.239 2.263 2.277 2.281
N" 0.941 1.903 2.846 3.797 4.761 5.737 6.723 7.719
rmin
VO
1.765 1.199 0.869 0.670 0.536 0.441 0.371 0.317
-5.715 10 -8.408 06 -11.3781 -14.6760 -18.2988 -22.2456 -26.5161 -31.1100
By setting (23) equal to (24), it is possible to get one equation which allows us to obtain B, from one given value for B, x3
+ a,x + a2 = 0
(25)
where ai
= ( P - P O Y ~ ) / ( P O Y-~ Q)
allows one to write V,, in the formIs
s
d7-
1V(r)
p(7)
d7
which, under spherical symmetry and using eq 18, gives VN,
+ 2V,,
= -4n
L=
r2V(r) p ( r ) d r
(21)
Introducing the expressions (6) for p ( r ) and (15) for V(r), solving the integral, and arranging terms, the result is
To find the four parameters from those conditions the next algorithm has been devised: (i) Set an initial B, trial value from a previous least-squares fit to the valence (rmin-m) p ( r ) part using the numerical OPM density. (ii) Given the OPM values for po, rminrand VN,, and the B, in (i), compute a,, a2 in (25) and solve this equation for x (B,). Due to the setting of this equation, the solution B, thus obtained fulfills the normalization condition and reproduces the VN, energy, with A, given by either (23) or (24) and A, = po - A,. (iii) For every B, one has a complete set of parameters in (ii). Making use of eq 22, compute the VN, 2V, energy and compare with the corresponding OPM value. (iv) Change E, by successive lesser increments and repeat steps (ii) and (iii) for every B, until the OPM VN, 2Ve, sum is reproduced up to a desired level of accuracy (chosen as lo4 au). For all the atoms studied, the final B, from which final values for the other parameters are fixed happens to lie within a small interval around the least-squares B, in step (i). For example, separate least-squares fits to the core and valence parts of the OPM electron density in fluorine, gives A, = 427.788, B, = 16.2130 (correlation coefficient 0.9975) and A, = 4.732 49, B, = 2.670 59 (correlation coefficient 0.9994). Final parameters from the above algorithm are A, = 441.705, B, = 18.5694, A, = 6.818 64, B, = 2.86775. This general resemblance can be interpreted as an additional support for the functional form (6) as a reasonable one to represent total electron densities in atoms. Final parameters for Li-Ne atoms are listed in Table IV. The model density (6) with these parameters is thus normalized, presents the correct value at the nucleus, and reproduces the proper OPM nuclear-electrons attraction energy. As far as exchange effects have been ignored, it also reproduces (along with its derived electrostatic potential) the interelectronic repulsion energy. Note, however, that no calculation of either N,, N , charges or rmlnradii has been imposed as condition for determining the parameters. It is therefore expected that there will be some deviation from the corresponding OPM values. Core and valence charges cal-
+
+
where BT = B, + B,. Summing up, we have sufficient information to determine the four parameters in (6) imposing the next conditions: 1. Normalization condition, eq 8. 2. Reproduction of the electron density at the nucleus, eq 9. 3. Reproduction of the VN, energy, eq 19. 4. Reproduction of the V,, energy, eq 22. Writing B;' = x, BC1 = y, -VN,/kz = Q, Z/8n = P,and using (9), the parameter A, which satisfies conditions 1 and 2 can be expressed in terms of the exponents B,, B, as A, = ( P - POv3)/(X3 - Y 3 )
(23)
Equivalently, from condition 3, A, =
( Q - P O ~ ~ ) / (-Xy~2 )
(15) Politzer, P.J . Chem. Phys. 1980, 72, 3027.
(24)
The Journal of Physical Chemistry, Vol. 95, No. 26, 1991
Electron Densities and Electrostatic Potentials TABLE VI: Cusp Condition, Eqs 26 and 27, for the Model Electron Density atom -22 cusp re1 error, % Li Be B C N 0
F Ne
-6 -8 -10 -12 -14 -16 -18 -20
-5.8025 -7.8616 -9.9353 -12.0226 -14.1185 -1 6.2214 -18.3307 -20.445 9
3.29 1.73 0.65 0.19 0.85 1.38 1.84 2.23
I 6.0
‘O’O
1
O%a
p ( r ) OPM
1
1 .o
I
2.0
3.0
r (bohr) Figure 3. Radial distribution function D(r) for fluorine obtained from OPM atomic calculations (solid line) and using the model electron density p(r) given by eq 6 (dashed line).
D(r) = 4nr2p(r) Solid:
D(r) = 4n?p(r) Solid:
2.0 .
CARBON
tI I
FLUORINE
Doshed: p(r) MODEL
I
tn
ttII\I
10657
p(r) OPM
Dashed: p ( r ) MODEL
TABLE VII: Electrostatic Potential V ( r ) (atomic units) at Some Radial Distances (bohr) Selected from the Logarithmic Mesh in OPM Numerical Calculations for C and F Atoms“ carbon V(r) fluorine V(r) r OPM model OPM model
Figure 2. Radial distribution function D(r) for carbon obtained from OPM atomic calculations (solid line) and using the model electron density p ( r ) given by eq 6 (dashed line).
culated with the model density via formulas (IO) and (1 I), respectively, rmi,values computed from (12) with the NewtonRaphson algorithm, and electrostatic potentials at the nucleus obtained from eq 16 are shown in Table V. These values must be compared to the corresponding OPM data in Table 11. As is apparent from those figures, the model p(r) produces a displacement of charge toward the core region of about 0.07 for all the atoms except Li. The position of minimum in D(r), rminris invariably shifted outwards, although it is closer to OPM values as the number of electrons increase. Electrostatic potentials Vo, on the other side, are in close agreement with OPM Voin Table 11. It must be emphasized that, although the reproduction of V,, energies is one of the conditions imposed to find the parameters in (6), OPM Vovalues correspond to the extrapolation r 0 of numerical electrostatic potentials, whereas Vodata in Table V have been calculated via formula (16). This good agreement can thus be viewed as a further test of validity for the model. Further information on the quality of the approximate p(r) function is provided by the analysis of the cusp condition derived by Steiner as a corollary to Kato’s theorem.I6 This condition establishes a relation between the values at the nucleus of p(r) and its derivative:
-
[ y]
= -2Zp(O)
r=O
Introducing the model density (6), one readily finds
- ACBC + AVBV = -22 Ac + A” Table VI gives this expression calculated with the final parameters under the column labeled “cusp”, along with its relative error. Considering the extreme simplicity of the functions used, these errors can be regarded as quite acceptable. The shape of the radial distribution function computed with the model p(r) is illustrated in Figures 2 and 3 for carbon and (16) Steiner, E . J. Chem. Phys. 1961, 39, 2365. Kato, T.Commun. Pure Appl. Math. 1957, 10, 151.
0.001 033 5793.4760 5793.4764 8685.7122 8685.713 1 0.002 186 2729.693 0 2729.693 2 4090.039 7 4090.0402 0.004626 1282.4678 1282.4679 1919.2104 1919.2108 0.009780 598.8626 598.8629 893.837 7 893.838 4 0.018 232 314.488 5 314.489 1 467.3760 467.3772 0.038404 141.871 5 141.873 3 208.823 5 208.8236 0.071 162 70.533 81 70.53680 102.5624 102.5475 0.130976 33.357 82 33.355 92 47.91 1 32 47,839 87 14.99050 0.238 221 14.96044 21.088 89 20.95624 0.424776 6.495 779 6.430 192 8.21 1667 8.148 847 0.734463 2.538464 2.514976 2.397501 2.471 872 1.281 752 1.OOOOOO 1.304001 0.931 950 1.004989 1.467033 0.420915 0.464842 0.201 152 0.228884 0.107 180 2.076761 0.134 132 0.031 814 0.036064 2.630 819 0.032 687 0.045 418 0.006 569 0.006 972 3.045431 0.013803 0.020516 0.002 101 0.002063 3.977443 0.002 120 0.003 535 0.000 176 0.000 136 5.034636 O.OOO274 0.000494 0.000011 0.000009
‘Model” values are computed with eq 15.
fluorine, respectively. The model curve presents a smoother behavior than the OPM one, especially in the transition from core to valence regions. The inner region and the valence outermost part are satisfactorily reproduced. The agreement between both functions improves, however, as 2 increases and thus, for fluorine the model curve shows a slightly better resemblance to its OPM counterpart. The analysis of electrostatic potentials reveals a very good agreement between OPM numerical data and values computed with the model expression (15). Table VI1 shows this comparison for carbon and fluorine as representative examples of a behavior found to be identical for the other atoms in the row. Some radial distances from the logarithmic mesh used in OPM calculations, covering a range of about 9 orders of magnitude in V(r),have been selected. As this table clearly illustrates, the model potential presents a satisfactory reliability, especially in the innermost atomic region where the ability of (15) for reproducing the reference OPM V(r) is remarkable. As far as the aim of this work is the construction of simple functions to represent atomic electron densities and electrostatic potentials, no consideration of exchange effects has been yet made. However, the main reason for choosing the OPM method as starting point for our model is the possibility of constructing in addition a separate exchange potential. The extension of this model to cover exchange terms as well as the proper analysis of such a potential is now being investigated and will be the subject of a forthcoming paper.
J. Phys. Chem. 1991, 95, 10658-10662
10658
Summary and Conclusions The known behavior shown by atomic electron densities in the form of exponential decrease with change of slope at a certain radial interval allows a simple representation as the sum of two exponential functions. The information of a previously obtained accurate density can thus be encoded into a functional form especially suited for fast computations of maps in systems containing large molecules. A semiempirical procedure to determine the parameters needed for such a representation has been devised in this work. Imposing the reproduction of nuclear attraction energies, interelectronic repulsion energies without exchange, and electron density values a t the nucleus along with the normalization condition, this procedure yields an optimum set of parameters for every atom. From this model electron density a simple analytic expression for the electrostatic potential V(r) has been derived through the Poisson equation. As a direct consequence, the model reproduces also the V(r) value a t the nucleus. Although no explicit reproduction of the electrostatic potential taken as reference (OPM)
is imposed, this procedure yields functions showing a remarkable agreement between both potentials, especially in the innermost radial region. Since the atomic calculations chosen as starting point correspond to a local potential formulation with quality of near-Hartree-Fock limit accuracy, the densities and electrostatic potentials under consideration are highly reliable. From that formulation it is also possible to find an empirical potential compatible with the model presented in this paper to represent the exchange contribution. This task along with the extension to atoms with higher principal quantum numbers is currently being developed at this laboratory. The true proof of validity for these simple representations of electron densities and electrostatic potentials will be, of course, their performance in calculations on large molecular systems. This point is now under study and will be presented in the near future. Acknowledgment. Financial support from the Direccion General de Investigacion Cientifica y Tecnica (DGICYT), under grant no. PS89-0025, is gratefully acknowledged.
Ab Initio Molecular Orbital Study of the GeH,O Potential Energy Surface Suk Ping So Chemistry Department, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong (Received: June 5, 1991)
The geometries of seven structures on the singlet GeH20 potential energy surface and triplet hydroxygermylene (HGeOH) have been optimized using the HF/3-21G(*) and MP2/3-21G(*) methods, and their vibrational frequencies have been computed. Electron correlation errors are corrected up to the MP4SDTQ level. It has been shown that the ground state of hydroxygermylene is a singlet and that, in contrast to the carbon analogues, germanone (H2GeO)lies higher energetically than hydroxygermylene. trans-HGeOH has been found to be slightly more stable than cis-HGeOH only after electron correlation is taken into consideration. With correction for x r q o i n t vibrational energies, MP4SDTQ activation energies for the reactions trans-HGeOH cis-HGeOH, H2Ge0 trans-HGeOH, and H2Ge0 H2 + GeO have been calculated to be 11.4, 34.2, and 66.6 kcal mol-', respectively.
-
-
Introduction The importance of electron-correlation effects in the study of reaction potential surfaces is well established.' Correlation is known to make a large contribution to covalent bond energies, and its neglect may lead to errors of over 50 kcal mol-' in the calculated values.2 In addition, activation energies, usually involving transition-state structures with partially broken bonds, are also very sensitive to correlation effect^.^ Extensive a b initio calculations on the structures of hydroxycarbene, namely cis-HCOH and trans-HCOH, have been reported4-10 including the use of C18 and MP'O techniques. In (1) Bauschacher, C. W., Jr.; Haber, K.; Schaefer, H. F. 111; Bender, C. F. J. A m . Chem. SOC.1977,99, 3610. (2) (a) Davis, J. H.; Goddard, W. A. 111. J . A m . Chem. SOC.1977, 99, 71 11. (b) Davis, J. H.; Goddard, W. A. 111; Harding, L. B. J . A m . Chem. SOC.1977,99,2919. (c) Harding, L. B.; Goddard, W. A. 111. J . Am. Chem .Sac. 1977. 99. _.. -. , . ., 4520. - -. (3) Pople, J. A.; Krishnan, R.; Schlegel, H. B.; Binkley, J. S . Int. J . Quant. Chem. 1978, 14, 545. (4) Pople, J. A.; Binkley, J. S . ; Seeger, R. Inr. J . Quunt. Chem. Symp. 1976 - - . -, -10-, -1 . (5) Krishnan, R.; Pople, J. A. Int. J . Quanr. Chem. 1978, 14, 91. (6) DeFrees, D. J.; Levi, B. A.; Pollack, S . K.; Hehre, W. J.; Binkley, J. S.; Pople, J. A. J . A m . Chem. SOC.1979, 101, 4085. (7) (a) Jaffe, R. L.; Hayes, D. M.; Morokuma, K. J. Chem. Phys. 1974, 60, 5108. (b) Jaffe, R. L.; Morokuma, K. J . Chem. Phys. 1976, 64, 4881.
-
addition, transition states for the following reactions have also been examined theoretically:8J*'2
-
trans-HCOH HzCO
-+
H2CO
cis-HCOH
(IC)
trans-HCOH
(2C)
H2+CO
(3C)
Recently, Withnall and AndrewsI3 carried out an infrared investigation of the photochemical reaction of germane and ozone in solid argon. Using the techniques of isotopic substitution and filtered photolysis, they identified, among other species, hydroxygermylene (HGeOH) and germanone (H,GeO) and reported a complete assignment of their infrared spectra. However, the experimental structures of these species are not yet known. Theoretically, Trinquier et al.'4*'5optimized the geometries of (8) Goddard, J. D.; Schaefer, H. F. 111. J . Chem. Phys. 1979, 70,51 17. (9) Lucchcse, R. R.; Schaefer, H. F. 111. J . A m . Chem. SOC.1978, 100,
298. -__.
(10) Harding, L. B.; Schlegel, H. B.; Krishnan, R.; Pople, J. A. J. Phys. Chem. 1980, 84, 3394. (11) Miller, W. H. J . A m . Chem. SOC.1979, 101, 6810. (12) Scuseria, G.E.; Schaefer, H. F. 111. J. Chem. Phys. 1989, 90, 3629. (13) Withnall, R.; Andrews, L. J. Phys. Chem. 1990, 94, 2351. (14) Trinquier, G.; Pelissier. M.; Saint-Roch, B.; Lavayssiere. H. J . Organomer. Chem. 1981, 214, 169.
0022-3654/91/2095-10658%02.50/0 0 1991 American Chemical Society
