J. Phys. Chem. 1992, 96, 7294-7301
1294
Simple Analytical Representation of Atomic Electron Charge Denslttes, Electrostatic Potentials, and Local Exchange Potentials L.Fernandez Pacios Departamento Quimica y Bioquimica, ETSI Montes, Universidad Politecnica, 28040 Madrid, Spain (Received: February 27, 1992)
A procedure presented before (Pacios, L. F. J . Phys. Chem. 1991.95, 10653) to represent atomic electron charge densities p(r) and electrostatic potentials is here generalized to include also local exchange potentials and extended to cover main-group elements up to krypton. Based on physically suitable simple functions to describe p(r), a method to reprcsent atomic potentials explicitly treating the classical (electrostatic) and quantum (exchange) parts is developed. The corresponding parameters for main-group atoms H-Kr are given, and the reliability of the whole model is discussed.
Introduction It is a fact now well established that the spherically averaged charge density p(r) for atomic ground states is a monotonically decaying function with a unique maximum located at the nucleus, p(r=O) = po, and the rate of the exponential decrease changing at certain radial intervals.’ These intervals happen to enclase the positions of minima in the radial distribution function D(r) = 4 d p ( r ) , fact used by Politzer and Parr2 to define a boundary radius, r h , between core and valence regions for first-row atoms. Very recently? these and other criteria (as for example those based on Vzp(r))have been analyzed to define atomic shell boundaries. In the first part of this work: an approximate procedure to represent atomic electron densities and electrostatic potentials based on the behavior of p(r) described above was presented for Li-Ne atoms. Since there is only a minimum in D(r) (apart from r = 0 and r = -) for these atoms, it seems appropriate to take a sum of two exponential terms, one for the core part and the other for the valence, as a minimum choice for an approximate representation of &). The aim is to develop a procedure simple enough to compute electron densities and related properties (electrostatic potentials, for example) in large molecults in an easy and fast way while maintaining the basic information from accurate ab initio atomic results encoded into the analytical functionS chosen. At present the utility of electrostatic molecular potentials as powerful tools in studies concerning biomolecular systems is und ~ u b t e d .For ~ small molecules it is relatively easy to obtain them directly from ab initio quantum calculations up to a high degree of accuracy. For large biomolecules, however, this approach is clearly prohibitive. It must be taken into account that the electrostatic potential usually needs to be computed numerous times for different spatial orientations of the molecule in a determined environment.6 The ultimate goal of the procedure presented in this work is the construction of atomic instead of electronic functions as building blocks for the further definition of intermolecular potentials based on atomic pair interactions. Apart from this evident computational simplification, another appealing feature is that one has an explicit analytical form for every atomic contribution to the interaction being considered. Moreover, as shown in this work, the electrostatic and exchange electronic effects can be explicitly described so that the model atomic potential obtained comprises the representation of both classical (electrostatic) and quantum (exchange) parts. This conceptually desirable feature gives the procedure a clear advantage over other semiempirical approximations. The general formulation of the model is given in next section. Although previously shown for Li-Ne atoms: the procedure is here generalized and extended in a revised form to main-group elements up to Kr. This section gives the definition and construction of the atomic functions used to represent atomic electron densities and electrostatic potentials. The analysis of results and performance of the model in atoms is presented in the third section.
The fourth part of this work is entirely devoted to the construction of an empirical exchange potential derived within this model. A brief summary and outline of the more important conclusions are finally presented. Atomic units are used through the paper. CoastructiOa of Atomic Functions for p ( r )
The starting point for the development of the procedure here presented is the atomic wave functions determined from the optimized potential model (OPM) by Talman et al.’ In this treatment a system of selfamistent equations is derived, in which an effective local central potential U(r)is obtained as the solution of a linear integral equation. If exchange terms in (H) are initially neglected, the solution of the problem is U(r)= U&, where this last term is just the Hartree potential arising from the charge distribution. The actual integral equation derived is one for U(r)-UH(r)and can be interpreted as an equation for an effective local central exchange potential Uex(r). In the self-consistent scheme, U(r)is varied to minimize (H) so that the final potential is the optimum in a variational sense. The whole exposition of the OPM method and discussion of results may be found in ref 6.
There are two reasons to choose the OPM procedure as the starting point for our model. First, it provides explicit separate radial electrostatic and exchange potentials in numerical form. Second, despite its relative simplicity and easiness of use, it shows an accuracy at the HartreeFock (HF) level, as was demonstrated in the previous paper! The strategy for the development of our method is the definition of a certain analytical model p(r) upon which all the atomic terms should be constructed. The corresponding parameters in this model density will be determined by setting up certain conditions imposed with respect to the reference (OPM) atomic results. Before fixing any particular analytical form, let us state the general framework where our model density will be defined. Recalling the atomic virial theorem VT = -2T, where VT and T are respectively potential and kinetic energies, and since the total energy may be written as E = T + VT,one has E = ‘/2VT. If correlation effects are at this initial stage ignored, the potential energy can be expressed then as
V’ =
VN,
+ V, + Vx
(1)
where V,, is the nuclear-electrons attraction energy, V, the interelectronic Coulomb repulsive interaction, and V,the exchange contribution. The other initial property important for our purpases, the atomic electrostatic potential, is defined as V(r) =
-r -
Ii - i l
As it is well-known, V(r)is related to the electron density through the Poisson’s equation V2V(i) = 4~rp(i),which under spherical symmetry can be written in the form
0022-3654/92/2096-1294S03.0Q/O 0 1992 American Chemical Society
The Journal of Physical Chemistry, Vol. 96, No. 18, 1992
Electron Charge Densities and Electrostatic Potentials
azv(r)
r
2 av(r)
a9 +-=r ar
. I
+ Bir) exp(-Bir)
(3)
4rp(r)
If the electron density is somehow known, the electrostatic potential V(r) may be easily obtained by means of a certain "screening f ~ n c t i o n "customarily ~-~ defined as 4(r) = rV(r)/Z. Thus, V(r) = Zt$(r)/r, and eq 3 is modified to
Substituting the model density into eq 6 and recalling A M Pexp(-Br) dr =
n1 -
B"+' '
B>0
( n integer positive), the nuclear attraction energy is
(4)
The relevant fact here is that all the contributions to the atomic total energy written as E = I/*( VN, V, V,) can be expressed in terms of V(r)and p ( r ) (the exchange potential V, will be the subject of the next section). The nuclear-electrons attraction energy is
+ +
+
The expression for the sum VNe 2V, as given by eq 10 is slightly more involved to obtain. Substituting (1 1) and (12), collecting terms, and solving the various integrals making use of (1 3), one gets the final expression M M AiAj(2Bi + E,) VN, 2V, = -32d-C-C (15) i = ~ j =B?(Bj ~ Bj)3
+
or, under spherical symmetry VN,
= -4uZJmrp(r) 0 dr
(6)
+
So far we have been considering energy relations. Now properties directly related to the density are analyzed, thus, the normalization condition for the model p ( r ) is
On the other side, the Coulomb interelectronic potential energy is The density at the nucleus is simply
(7)
M
Po
If the eq 2 is multiplied by p ( 9 and integrated, one obtains
This way, V, can be written in the forms
I i - tl
I
zp(r3 dir
V(r) p ( i ) d7 (9)
which, after introducing the definition ( 5 ) for VN, and imposing spherical symmetry, gives
VN, + 2V, = -4*Jm9V(r)p(r) dr
If a certain r, is taken as boundary between two atomic regions as done with the minimum in the radial distribution function (rk) for defining core and valence domains? it is possible to calculate two separate charges, say core N, = 4r.for4p(r) dr and valence N, = 4r.frmm9p(r)dr. If the indefinite integral exp(Br) exp(Br) dr = -[ l + (1 -&)'I (18) B3 is taken into account, the expressions needed to compute these charges are
M
exp(-Bir)
(11)
This analytical choice for p ( r ) is directly based on the exponential decaying behavior discussed above. For Li-Ne atoms two terms were used in ref 4;in this paper the length of this expansion is set to M = 3 for Na-Ar and M = 4 for K-Kr main-group elements (transition metals not included). We turn now to the determination of the parameters in (1 1). To accomplish this task, a number of conditions need to be imposed so that one had certain reference data (OPM results) reproduced. After analytically fixing p(r), it is direct to set all the atomic properties in tenns of (1 1). Thus,solving for 4(r) with the model p(r), eq 4 gives the "screening function" r
The electrostatic potential is then
M
Ai
N, = 4 r C - ( 2 i=l E?
- exp(-Bir,)[(l
M Ai N, = 4u-C- exp(-Bir,)[(l i=lB?
(10)
Up to this point, the two potential contributions to the total energy have been written in terms of p(r) and V(r),both mutually related through eq 3. Now we introduce here for the electron density the functional form used in the previous paper? that is, a sum of exponentials: i
(17)
i= 1
19
did?'= 2V, =
p(r) =
= CAi
+ Bir,J2 + 11) + Bir,,.J2 + 11
(19)
Adding up these two formulas, it is immediate to confirm NE+ N, = Z, with Z given by eq 16. This radial value r, must be a minimum in the radial distribution function D(r), 90 that finding the derivative ( w 4 / w r = ? m= 0 with D(r) = 4 d p ( r ) and the p ( r ) model (eq 11) leads to the relation M
CAi(2 - Birm)exp(-Birm) = 0
i= 1
(20)
An additional test of reliability for the electron density is also provided by the radial atomic expectation values (F),p =
-2,-1,1,2:
We have thus sufficient analytical relations to determine the parameters which reproduce the reference OPM atomic results. Actually, there are more conditions than needed, so that depending on the length M chosen in (1 l), some of them will serve for additional tests of performance of the model. As a consequence, a considerable degree of arbitrariness remains at the choice of the concrete procedure to fix the set Ai, Bi.
72%
The Journal of Physical Chemistry, Vol. 96, No. 18, 1992
Since the main objective of this model is the development of atomic functions to be further implemented in calculations on intermolecular interactions, special care must be taken of the reproduction of the long-range region in the atomic electron density. Thus, rather than accurately reproducing a vast list of atomic properties, we have searched a model density showing a satisfactory resemblance to the overall (especially the outer) shape of the OPM density, while keeping an acceptable precision in the computed atomic properties. It must be taken into account that the use of extremely simple functions is the essential restriction which imposes an equilibrium among all the possible requirements. After an exhaustive exploratory work, the next basic conditions have been finally maintained for all the lengths M considered: reproduction of the OPM V,, and V . energies as given by eqs 14 and 15, normalization condition (16), and reproduction of the OPM density at the nucleus (17). For Li-Ne atoms M = 2, and the four parameters in (1 1) are completely determined from these four conditions; this case was the subject of the previous paper: where an algorithm devised to find numerically the parameters was presented (this procedure is here applied also to helium). For Na-Ar atoms, M = 3, and six parametem are then needed. In this case, the above conditions are used to fix the three linear Ai parameters and one of the Bi. The other two are varied by successive lesser increments until: (i) the form of the resulting D(r) function, and (ii) the accuracy of the f m l atomic properties computed for every whole set {A&}, are satisfactory. To this end, a computer program has been written to systematically vary the "free" Bi parameters, compute the atomic properties in terms of the expressions above presented, and plot the resulting radial distribution function compared to the reference OPM. In all cases, the initial guess for the arbitrary Bi parameters is made from a previous nonlinear optimization fit procedure which uses the m r r r~utine.~ For K-Kr atoms we set M = 4 and eight parameters are therefore needed for p(r). In addition to the above four general conditions, the reproduction of (9) is also imposed, on one side to reduce the number of B, parameters to vary, on the other side to ensure an overall correct shape in atoms like these, where the radial distribution function presents many oscillations. In all cases it is worth mentioning that slight changes in the parameters of the final {A&} set cause large deviations in the results of the atomic properties. Thus, the model p(r) given by expression 11 with the final parameters listed in Table I represents an optimum compromise between simplicity and reliability for an approximate electron density which will serve as basis for the construction of atomic potentials. Arguments in support of this affirmation are presented in next section.
Discdon of R d t s The fmal set of Ai, Bi parameters in the model density (1 1) for the maim-group elements H-Kr obtained from the procedure presented in the last section is given in Table I. Li-Ne parameters given in ref 4 are here repeated for completion purposes, and helium is also included with M = 2. Since the exact solution for the H atom is = d/2 exp(-r), its density can be written in the form of our model density as p(r) = A exp(-Br), with A = d , B = 2. It is then immediate to verify Z = 8~A.6'= 1, VNe= 4~zA.6' = -1, E '/z( V) = '/2( V N ~ ) = -'I2,and radial expectation values as given in the general expression 21, ( f 2=)4?rAB1 = 2, ( f ' ) = 4 r A B 2 = 1, ( r ) = 2 4 r A 8 4 = 3/2, and (9)= 967rAB5 = 3, all exact values. The electrostatic potential in the general form (12) with M = 1 and A = d = 0.318 310, B = 2 is then for hydrogen
Table I1 shows atomic data reproduced by the model density: nuclear charges obtained from the normalization condition (la), density at the nucleus (17), and VNe,V, energies calculated with (14) and (15). Referencevalues have been obtained from pndous OPM calculation for every atom. In this method, the Coulomb interelectronic repulsion energy V, is explicitly given, whereas
Pacios
i
NITROGEN Solid OPM density Dashed Model density
Figure 1. OPM radial distribution function D(r) (solid curve) and D(r) computed with the model p(r), eq 11 (dashed curve), for nitrogen. 18.0
PHOSPHORUS Solid . OPM density Dashed: Model density
r (bohr)
Figure 2. OPM radial distribution function D(r) (solid curve) and D(r) computed with the model p(r), cq 1 1 (dashed curve), for phosphorus. 50.0 45.0
,
I
1
1 ARSENIC
Solid OPM density Dashed Model density
i
Lq 5.0
:
---__ r (bohr)
Figare 3. OPM radial distribution function D(r) (solid curve) and D(r) computed with the model p(r), eq 1 1 (dashed curve), for anenic.
nuclear attraction energies VN,have been calculated as the difference between the total potential energy and the sum (V, Vd, V,being the OPM exchange energy also oxplidtly detnmined. Densities at the nucleus p(0) have been obtained extrapolating the numerical OPM electron densities by means of a rational function extrapolation (Bulirsch-Star algorithm'0). For hydrogen, exact values are given. The shape of radial distribution functions D(r) = 4 d p ( r ) computed with the model density (1 1) and compared to their comsponding OPM functions arc plotted in F i i 1-3 for some
+
Electron Charge Densities and Electrostatic Potentials TABLE I: Psrameters Ai and Bi in the Model Electron Density (Eq 11) for the Atoms H-Kr H Ai 0.318310 Bi 2.000000 1.982685 1.612575 He Ai Bi 4.840638 2.961300 13.76849 0.045482 Li Ai Bj 5.818505 0.972372 35.24546 0.141982 Be A, 7.888538 1.175743 E, 71.52347 0.410814 B A, 9.983889 1.478643 Bi C Ai 126.5168 1.005506 12.10368 1.817852 Bi N A, 204.0372 2.112261 Bi 14.24224 2.165886 0 A, 307.8554 3.959890 16.39772 2.516477 Bi F Ai 441.7051 6.818640 Bj 18.56940 2.867753 Ne Ai 609.2808 10.99983 E, 20.75692 3.219070 Na Ai 811.6359 22.64350 0.066924 23.52895 4.198842 0.987565 Bi 1055.512 38.90142 0.242820 Mg Ai B, 26.22141 5.124824 1.231856 1351.191 51.83015 0.324150 AI Ai E, 28.37559 5.589321 1.261974 1701.654 64.19435 0.535419 Si A, 30.32586 5.983461 1.397054 B, 2099.882 86.60264 0.990446 p Ai 32.62216 6.675897 1.592433 Bi 2550.706 118.5469 1.715 107 s Ai 35.08001 7.509562 1.792874 Bi 3054.353 163.9730 2.838155 CI Ai 2.003272 37.70089 8SO6661 Bi 3612.514 225.5703 4.364711 Ar Ai 40.46129 9.608664 2.204522 Bi 4304.480 228.8470 6.493966 0.078027 K Ai 42.03668 9.487663 2.705441 0.914152 Bi 5028.654 281.4227 10.89296 0.315758 Ca Ai 44.74341 10.16653 3.477912 1.140437 Bi Ga Ai 19728.65 630.7752 53.15526 0.325868 64.42690 12.66935 4.248192 1.277828 Bi 21738.24 683.8334 75.96176 0.635832 Ge Ai 13.46359 4.770980 1.426197 66.36155 Bi 23893.28 717.3065 110.8526 1.246079 AS Ai 1.610779 68.25318 14.26108 5.429911 Bi 26171.07 762.4529 153.9620 2.062352 SC Ai 70.29363 15.19501 6.065248 1.772866 Bi 28571.57 823.1301 205.8293 3.077716 Br Ai 72.43167 16.33326 6.671889 1.919073 Bi 31104.75 880.1117 279.8113 4.674739 Kr A, 7.439718 2.085725 74.83648 17.52711 Bi
selected atoms. As is apparent from these plots, the model p ( r ) curve reproduces well both the innermost and outermost parts of the reference OPM curves, a fact found in all the atoms studied in this work. The simplicity of the analytical model density imposes a severe restriction to the possibility of reproducing the intermediate oscillations, and thus the model curves show much smoother transitions from inner (core) to outer (valence) regions. It must be remembered, however, that the precise reproduction of electron density shapes is not the main purpose of this model. Directly related to the existence of different regions in radial distribution functions are the properties presented in Table 111. rml is the first position of minimum in model D(r) functions, determined solving eq 20 with the Newton-Raphson algorithm. The partition of charge into core, No and valence, N,, contributions is set taking as boundary between these regions the outermost minimum in D(r), rmin,as prescribed by P o l i t ~ e r .The ~ ~ ~OPM rmhvalue is almost coincident with that tabulated by Boyd" with H F limit densities as it has been demonstrated beforeS4Partial charges given in second and third columns of Table I11 are obtained by equaling this radius to rmin (19) and compared to the corresponding Nc, N, OPM values in last columns, determined by numerically integrating OPM electron densities. For first-row
The Journal of Physical Chemistry, Vol. 96, No. 18, 1992 7297
TABLE Ik OPM Atomic Reference Data Reproduced by the Model D-iW
atom H He Li Be B
C N
0
F Ne
Na Mg AI
Si P
S C1 Ar
K
Ca Ga Ge As Se Br Kr
V,
Z
P(0)
VNe
1 2 3 4 5 6
0.3183lob 3.595260 13.81397 35.38744 7 1.93428 127.5223 206.1495 31 1.8153 448.5237 620.2806 834.3463 1094.6567 1403.3458 1766.3833 2187.4746 2670.9683 3221.1643 3842.4485 4539.8994 5321.2856 20412.902 22498.676 24722.687 27089.543 29603.609 32269.348
-1 .ob -6.74824 -17.14531 -33.63223 -56.89054 -88.05590 -128.09129 -177.96443 -238.64408 -3 11.09968 -389.67725 -478.98182 -578.42563 -689.23534 -81 1.84012 -946.68589 -1094.2239 -1254.9097 -1422.7340 -1602.9579 -4607.5423 -4969.6306 -5347.8209 -5742.3841 -6 153.5976 -6581.7347
7 8 9 10 11 12 13 14 15 16 17 18 19 20
31 32 33 34 35 36
2.051 19 4.06186 7.15476 1 1.57849 17.74622 25.96754 36.55921 49.83988 66.12892 80.00402 95.77867 112.78035 131.86486 153.14580 176.75687 202.83901 231.53682 257.20605 284.78784 834.83546 896.71333 961.25087 1028.5552 1098.7359 1 171.8950
a V,, is the nuclear electrons attraction energy, V, the repulsive interelectronic Coulomb energy, and p(O), the total density at the nuclcus. All values are in atomic units. bExactvalues.
TABLE IIk Position of First Minimum in the RadiDl Distribution Function D ( r ) ,rml(bohr), and Partition of Charge into Core, N, and Valence, N, Parts, Taking as Boundary Limit the Outermost Minimum in the OPM D ( r ) Function for Ns-Kr Atoms (See Text) model d r ) OPMa
atom Na Mg A1
Si P
S CI Ar K
Ca Ga Ge
As Se Br Kr
rml
0.266 0.231 0.211 0.199 0.184 0.170 0.157 0.146 0.137 0.123 0.089 0.087 0.084 0.082 0.079 0.076
NC 10.236 10.247 10.279 10.411 10.424 10.396 10.339 10.280 10.615 10.610 11.397 1 1.367 11.404 11.420 11.407 11.513
N" 0.764 1.753 2.721 3.589 4.576 5.604 6.661 7.720 8.385 9.390 19.603 20.633 21.596 22.580 23.593 24.487
rml
0.269 0.240 0.216 0.195 0.179 0.164 0.152 0.141 0.132 0.124 0.074 0.071 0.068 0.066 0.064 0.062
Nc
N"
10.237 10.322 10.326 10.326 10.327 10.330 10.330 10.330 10.337 10.345 10.942 10.972 11.013 11.033 11.039 11.073
0.763 1.678 2.674 3.674 4.673 5.670 6.670 7.670 8.663 9.655 20.058 21.028 21.987 22.967 23.961 24.927
OPM calculations, this work.
Li-Ne atoms there is only one minimum, and the core-valence partition is unique. rml,N,, N , values for these atoms were presented in Tables I1 and V in ref 4 and are not repeated here. The analysis of these properties demonstrates some interesting features. Although no precise reproduction of charges should be expected, there is a reasonable agreement between model and OPM charges, especially good if one thinks of the fact that middle regions in D(r) are precisely the worst represented by the a p proximate p(r) functions. For Li-Ne atoms the differences in valence charges are invariably lower than 0.08 units: for Na-Ar lower than 0.10, while for K-Kr amount up to around 0.40 units. Two considerations are pertinent for the latter case: on one side, we have now radial distribution functions showing more complex behavior and thus the approximation in using the model function (1 1) with M = 4 is cruder. On the other, the valence charges amount to higher values, greater than 20.0 for Ga-Kr atoms. Differences in rmlvalues are lower than 0.07 bohr in Li-Ne and
7298 The Journal of Physical Chemistry, Vol. 96, No. 18. 1992
Pacios
TABLE Iv: R.ctirl Exmetation V d m (Atomic Udb)from Model Density C.lerhtiom for Atom H-Kr
Hb He Li
Be B
C N 0 F
Ne Na Mg
AI
Si P
S
c1
Ar K
Ca Ga Ge As
se Br
Kr
2.0 11.99 30.32 57.66 93.52 138.3 192.3 255.7 328.8 411.8 502.1 603.7 718.1 844.8 979.7 1124 1278 1442 1621 1803 4634 4960 5297 5643 5998 6355
1.o 3.374 5.715 8.408 11.38 14.68 18.30 22.25 26.52 31.11 35.43 39.92 44.49 49.23 54.12 59.17 64.37 69.72 74.88 80.15 148.6 155.3 162.1 168.9 175.8 182.8
1.5 1.853 4.742 6.288 7.022 7.387 7.61 1 7.766 7.882 7.971 11.00 12.37 13.80 14.53 15.04 15.45 15.76 16.03 19.80 21.77 23.45 24.30 24.96 25.48 25.89 26.29
3.0 2.361 16.40 19.41 17.75 15.42 13.47 11.91 10.66 9.645 26.75 29.16 33.43 32.89 31.16 29.44 27.65 26.12 5 1.28 56.61 41.02 42.24 42.11 41.46 40.57 39.59
11.98 30.17 57.49 93.40 138.3 192.5 256.1 329.3 412.5 506.4 610.8 727.4 852.4 987.8 1134 1290 1456 1639 1828 4593 4910 5237 5576 5926 6287
3.374 5.713 8.409 11.38 14.68 18.30 22.25 26.52 31.12 35.44 39.93 44.5 1 49.25 54.15 59.20 64.40 69.75 74.94 80.21 148.7 155.4 162.1 169.0 175.9 182.9
1.855 4.693 6.111 6.810 7.196 7.451 7.637 7.780 7.894 10.85 12.26 13.72 14.55 15.11 15.52 15.83 16.08 19.46 21.26 23.41 24.28 24.93 25.45 25.89 26.26
2.370 15.17 17.12 15.82 14.06 12.53 11.28 10.24 9.381 27.19 29.64 33.51 32.81 31.18 29.39 27.67 26.07 5 1.28 56.61 41.02 42.24 42.11 41.46 40.57 39.59
"OPM calculations, this work. *Exact values.
TABLE V clgp Condition for Model Density Cdculatioas
atom H He Li Be B
C N
0 F Ne Na Mg A1
Si P
S CI Ar
K Ca Ga Ge As
se
Br Kr
CUSP
-22
error (%)
-2.000 -3.998 -5.803 -7.862 -9.935 -12.023 -14.119 -16.221 -18.331 -20.446 -23.003 -25.446 -27.528 -29.432 -3 1.58 1 -33.835 -36.183 -38.607 -40.339 42.828 -62.670 -64.544 -66.402 -68.373 -70.407 -72.678
-2 -4 -6 -8 -10 -12 -14 -16 -18 -20 -22 -24 -26 -28 -30 -32 -34 -36 -38 -40 -62 -64 -66 -68 -70 -72
0.0 0.05 3.3 1.7 0.65 0.19 0.85 1.4 1.8 2.2 4.6 6.1 5.9 5.1 5.3 5.7 6.4 7.2 6.2 7.1 1.1 0.85 0.61 0.55 0.58 0.94
lower than 0.01 bohr in Na-Ca, whereas as expected, larger errors are again encountered in Ga-Kr: 0.01 5 bohr in values typically around 0.070 bob. The overall agreement, however, is reasonably good if one considers that no condition concerning spatial partition has bem impoaed on the construction of the approximate functions. Radial expectationvalues computed with (21) are given in Table IV compared to the corresponding OPM numerically calculated integrals. Despite the approximate nature of the electron density to which these radial properties correspond, the overall agreement is remarkable. ( f 2 ) errors are less than 0.5% for Li-Ne and around 1.1% for Na-Kr. ( i l ) deserves a special comment. Since VN, is defined by eq 6, it is immediate to check than ( f ' )= -VNJZ and thus the condition of reproducing this potential energy implies the reproduction of (r-I). It must be emphasized however that OPM VN, values, as stated before, are not computed directly
and that ( i l ) integrals are numerically calculated from the radial grid for which the OPM density is computed. Thus, the extent to which the above relation is satisfied can be viewed as a further test of validity for the model. Whereas ( f Z and ) ( f l ) are properties which depend heavily on the behavior of the electron density in the innermost radial density, ( r ) and (9) are most influenced by medium- and long-range radial values. For Li-Ne atoms the four parameters in the model density are determined by the four conditions which (in particular those concerning p(0) and VN,), favor the reproduction of the inner regions, and thus these atoms show much better (iZ) and ( f ' ) values than ( r ) and ( 3 ) .Na-Ar atoms have already six parameters, and the correspondingly higher flexibility allows the improvement of ( r ) , (9) values. For fourth-row atoms K-Kr the additional condition of reproducing (9) implicitly produces a very good agreement in ( r ) . It must be emphasized again that the restriction of the extreme simplicity in p(r) precludes an exceedingly good precision simultaneously in energy-related properties (which favor behavior in the inner atomic region) and long-range density related prop erties. Considering the number and the stringent nature of the tests so far analyzed, the overall performance of the model density seems quite acceptable. Further information on the quality of the model representation of p(r) is provided by the cusp condition. As a corollary to Kato's theorem, Steiner12established the relation between the atomic density at the nucleus, p(O), and its derivative:
(aP(r)/ar),=o= -2a40) (22) If the model density (1 1) is introduced, the cusp condition gives M
M
-22 = C(AiBi)/CAi i= I
i- I
(23)
It must be emphasized that the expression (22) is a completely general condition to be fulfilled by the exact wave function obtained for a spinless n-electron exact Hamiltonian from the nonrelativistic Schrodinger equation for any atomic state.12 Table V presents this condition for our model density. The column labeled Cusp gives the value of the quotient in (23) for all the atoms considered in this work. The right column gives the relative error with respect to the 2 2 value. In any case this error exceeds 7 . 5 2 , and 13 of the 25 atoms included in this table (H is not counted) present errors lower than 2.0%. If the rather approximate
The Journal of Physical Chemistry, Vol. 96, No. 18, I992 7299
Electron Charge Densities and Electrostatic Potentials
Solid OPM Electrostatic Potential Dashed. Model Electrostatlc Potential
Solid : OPY Electrostatic Potential Dashed: Model Electrostatic Potential
NITROGEN
ARSENIC
-.
'\
10 - 7 1
00
20
40 60 r (bohri
80
i
PHOSPHORUS
-'
10 0 r
Flgm 4. OPM elcctmstatic potential V(r)(solid curve) and model V(r), eq 12 (dashed curve), for nitrogen.
Solid : OPM Electrostatic Potential Dashed: Model Electrostatic Potential
10
i
i
(h-l-!
Figwe 6. OPM electrostatic potential V(r) (solid curve) and model V(r), eq 12 (dashed curve), for arsenic.
this approach the exchange potential is related to the electron density through the approximated expression
where atomic units are used and where a is a parameter for which a number of procedures exist'5 in order to fix its value. Computational methods based on the Xa approximation have been succtssfully applied to a wide range of atomic and molecular problems but such exchange potentials exhibit the wrong behavior Uex(r) 0 at large r, since p ( r ) 0 at long distances. As a continuation to our model density (1 1) we propose a functional relation of the type (24), say Vx(r)= -C, [ ~ ( r ) ] ~but /~, in order to keep the correct behavior at large r, the term - exp(-Or)l/r (25) is added. This behaves like -l/r when r m and tends to -6 when r 0. Our model exchange potential is then empirically defined as 1 - exp(-@r) V,(r) = -Cx[p(r)I1l3(26) r The exchange energy for this potential is
-
Flgme 5. OPM electmetatic potential V(r) (solid curve) and model V(r), eq 12 (dashed curve), for phosphorus.
nature of the model is taken into consideration, this result can be regarded as quite acceptable. FQum 4-6 show electrostatic potentials V(r)plots for the same atoms of the F i i 1-3. OPM numerical potentials (solid curves) are plotted together with model V(r)potentials (dashed curves) computed with eq 12. Since about 10 o r d m of magnitude in the eledraatatic potentials are covered for the range of radial distances present in OPM numerical grids, a logarithmic scale is chosen for displaying V(r). As seen, the agreement is fairly good, especially for N and h. For all the atoms studied, both potentials are virtually indistinguishable in the inner atomic region and differ very little for medium and long distances. The phosphorus case is representative of the worst agreement found in all the V(r) analyzed, and even so the differenca are noticeable only when the electrostatic potential has fallen down to values around 0.001 au.
Model Exchtge Poteatid Asetatad at t h e b q h h g of the second section, the main reason to &me the OPM atomic wave function as reference for building our model is that it provides an explicit local exchange potential UJr) while keeping the overall accuracy comparable to the H F limit level. The essence of the OPM procedure is precisely the ConStruCtion of a total effective potential which at long distances is completely dominated by the exchange term.' Moreover, one of the most important advantages of this formulation over other approximate treatments of exchange effects is that it shows at large r, the correct behavior Uex(r) -l/r (in atomic units). In local density schuna3" thereare several ways to approximate exchange potentials and energies: one of the best known is the Xa method by Slater," based on free-electron gas models. In
-
-
-
-
Ex = (Vx(r)) =
Y 2 J h
Vx(r) fl
(27)
which under spherical symmetry and substituting (24) gives
E, = - 2 ~ C ~ J ~ $ p ( r )dr~ / ~2xJmrp(r)[ 1 - exp(-@r)] dr (28) The first integral with p(r) given by (1 1) has been numerically calculated by using a modified Simpson d e l 6 with a practical infinity of 50.0 bohr. This routine gives numerical results coincident up to the seventh decimal place with that provided by the numerical integration method implemented in the MATHEMATICA package.l' Let us denote this numerical integral by Z(p4I3). The second term in (28) can be easily obtained in analytical form; in doing so, one finds for the total exchange energy I
E, = -2*CXZ(p4/') - 2 x e A j l 1-1
Biz
-
1
-1 (Bj
+ /5)2
(29)
Since at this point the only contribution to the total atomic energy not determined yet in our model is precisely the exchange one, we impose to (29) the condition of reproducing OPM E, values. This fmes one of the two parameters in the model exchange potential (26), say C,. To determine the parameter we use the fact of that at r = 0 this potential has the finite value V,(O) = -Cx[p(O)]'l3 - (3 (30) which may be compared to the corresponding V,(O)obtained for
7300 The Journal of Physical Chemistry, Vol. 96, No. 18, 1992 00 -1 0
-2
NITROGEN
1 I
-30
TABLE VI: Atomic Exchange Energies, E,, Vduea at the Npfkw for Model, V,(O), rsd OPM, Vxom(0), Exchange Potentials, rsd C, .ad B Panmeten in tbe Modd Eq 26 for He-Kr Atoms (All Vdues in Atomic Units) atom E." VAO) V,OPM(O)" c. 8
c
1
j /
i
l
He
r
', '
I
i
Li
Be
1
h
-60 -7 0
i
B C N 0 F
__----
__...
Pacios
_ _ - - +
model
Ne 1
-9 0 0 001
,
,
,
,
,
,
01
0 01
.
I
,
1
,
,
.
.
I
Na Mg A1
.
10
1
r (bnhr)
Si
P S C1 Ar
00
I , , , ,
,
~
I , , , , , ,
,
,
'
I
,
'
'
,
,
,
'
, ' " ' , 3
' "
;..-I--I /,
f?
-20
-40
.-
-60
-
,--.
PHOSPHORUS
I '
1 f
L
v
',
-80
-10.0
/
1
Ca
1
Ga Ge
I'
I -8.0 -
.--'model '
"""'
0 01
'
i '
"""'
0.43641 0.45479 0.45137 0.47685 0.49680 0.51163 0.52665 0.54034 0.55384 0.55113 0.55144 0.56185 0.57397 0.58323 0.59161 0.59976 0.60847 0.61190 0.61020 0.69869 0.70458 0.70903 0.71274 0.71637 0.71810
1.21 1.77 2.23 2.60 2.92 3.22 3.49 3.75 4.00 4.44 4.82 5.14 5.43 5.72 6.00 6.27 6.52 6.81 7.13 8.95 9.13 9.32 9.52 9.72 9.95
i
I-llf
g; '
-1.68708 -2.30346 -3.13083 -3.97639 -4.77876 -5.57930 -6.38054 -7.18371 -7.98996 -8.57399 -9.37861 -10.2 176 -1 1.OS07 -1 1.8828 -12.7158 -13.5504 -14.3874 -15.0577 -15.8762 -24.9086 -25.7454 -26.5818 -27.4188 -28.2574 -29.0990
1
- 12.0 -14.O1""' 0.001
AS Se Br Kr
-1.87855 -2.86123 -3.71189 -4.58316 -5.42057 -6.24238 -7.06125 -7.88612 -8.72333 -9.62843 -10.5032 -11.4304 -12.3683 -13.2910 -14.2083 -15.1276 -16.0503 -16.9421 -17.7831 -28.0449 -29.0206 -29.9752 -30.9257 -31.8808 -32.8122
OPM calculations, this work.
f'
I
-
K
i
-1.02560 -1.78046 -2.66552 -3.74055 -5.00222 -6.45897 -8.12087 -9.99845 -12.10250 -14.00877 -15.98458 -18.05874 -20.24206 -22.53902 -24.95519 -27.49630 -30.16828 -32.65706 -35.18971 -73.47218 -77.39842 -8 1.38623 -85.44550 -89.58411 -93.80849
'
'
"""'
'
'
"""'
10
1
-12.0
ARSENIC
-
r (CAhr)
Figure 8. OPM local exchange potential V,(r) (solid curve) and model Va(r),eq 26 (dashed curve), for phosphorus.
the OPM numerical exchange potential by the same extrapolating technique used before.I0 The fl parameter is then determined under the condition of minimum error between the OPM V,(O) value and (30),while keeping the reproduction of the exchange energy. Table VI shows atomic exchange energies E,, exchange model potentials at the nucleus, V,(O), together with their OPM corresponding values, VxoPM(0),and C,, fl parameters for He-Kr atoms. In Figures 7-9 are plotted both model and OPM exchange potentials for the same atoms selected before for plotting purposes; the curve -l/r is included in every case to illustrate the large-r behavior. As is clearly apparent, the model V, potential reproduces reasonably well the overall shape of the reference w e , especially the transition between different domains. Wiggles and peaks in OPM curves are not reproduced: this is obviously expected from the simple formulation of our model exchange potential. It must be taken into account that we have imposed a priori functional forms under computational convenience criteria while maintaining the physical description of the essential nature of every atomic effect. Thus, it is of no interest for our purposes to reproduce exactly the kinks exhibited by exchange potentials, features on the other side only possibly meaningful within the particular OPM procedure. Summary pad Conclmiom The well-known monotonical decreasing behavior of the electron charge density in neutral atoms changing its rate of decrease at certain radial intervals allows its physically meaningful description in the form of sum of exponential terms. Taking this assumption as the basis for a simple analytical representation of atomic electron densities, a model has been presented which defines the
--16.0
i
-20.0
-24.0
-28.0
P J -
// /I
oPM
;-j/f
_ _ _ _ -_- _- - model
-320~'"'"'
0 001
i
'
'
"""'
0 01
I
I '
'
I
I
. . . . . . . . . . . . . . . . . . . . . . . .
01
7
(bonr)
IC
F i i 9. OPM local exchange potential V J r ) (solid curve) and model V&), eq 26 (dashed curve), for arsenic.
classical electrostatic potential and a local exchange potential as functions of that model p(r) density. The parameters in the functional form of p(r) are determined by imposing the conditions needed to reproduce selected results of previously chosen atomic calculations. All the properties of interest can be easily expressed in simple analytical form in terms of the model density and therefore in terms of those parameters. In this work, the OPM procedure by Talman et a1.7 has been chose as the reference atomic treatment due to two appealing features. First, it provides a variational effective local atomic potential explicitly constructed as sum of the two contributions, classical electrostatic and local exchange. Second, its accuracy is almost comparable to that H F limit. In despite of imposing the severe restriction of taking very few terms in the analytical expansion of p(r) (2 in H t N e , 3 in Na-Ar, and 4 in K-Kr), the determination of other atomic properties within the model demonstrate an acceptable reliability. But the most important conclusion is that electron densities, electrostatic potentials and exchange local potentials are very satisfactorily described by the model here presented. This lends support to its validity as ground for the further development of atomic pair
7301
J. Phys. Chem. 1992,96,7301-7307
potentials to be included in intermolecular interactions studies, the ultimate goal toward which this model is directed. The usefulness of this approximate p(r) and their derived p e tentials in local density type mol& calculations will be analyzed in the near future. In this context, the inclusion of the model into efficient techniques for the decomposition of intermolecular interaction energies is one of the immediate problems to be treated. Related to this application, the estimation of intra- and intermolecular correlation effects by means of any suitable local density existing approximation with our analytic model p(r) will be also explored. This could permit a fast and easy computation of approximate correlation energies without basis set dependency problems. Acknowledgment. Financial support from the Spanish Direccion General de Investigacion Cientfica y Tecnica (DGICYT), under Grant No. PS89-0025, is gratefully acknowledged.
References and Notes (1) Weinstein, H.; Politzer, P.; Srebenik, S.Theor. Chim. Acta 1975,38, 159.
(2) Politzer, P.; Parr, R. G. J. Chem. Phys. 1976, 64, 4634. (3) (a) Schmider, H.; Sagar, R. P.; Smith, V. H. J . Chem. Phys. 1991,94, 8627. (b) Kohout, M.; Savin, A.; Preuss, H. J. Chem. Phys. 1991,95,1928. (4) Pacios, L. F. J . Phys. Chem. 1991, 95, 10653. (5) (a) Politzer, P., Truhlar, D. G., Eds.; Chemical Applications of Atomic and Molecular Electrostatic Potentials; Plenum: New York, 1981, (b) Buckingham, A. D.; Fowler, P. W.; Huston, J. M. Chem. Rev. 1988,88,963. (6) Sokalski,W. A.; Sneddon, S. F. J. Mol. Graphics 1991, 9, 74. (7) (a) Talman, J. D.; Shadwick, W. F. Phys. Rev. A 1976, 14, 36. (b) Talman, J. D.Compur. Phys. Commun. 1989, 54, 85. (8) Politzer, P. J. Chem. Phys. 1980, 72, 3027. (9) Chandler, J. P. QCPE 1976, 13, 307. This routine is based on the SIMPLEX optimization procedure: Nelder, J. A.; Mead, R. Compur. J. 1965, 7, 308. (10) 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. (11) Boyd, R. J. J. Chem. Phys. 1977, 66, 356. (12) Steiner, E. J . Chem. Phys. 1963, 39, 2365. (13) Sham, L. J.; Kohn, W. Phys. Reu. 1966, 145, 561. (14) Slater, J. C. The Self-Consistent Field for Molecules and Solids; McGraw Hill: New York, 1974; Vol. 4. (15) Nagy, A. Inr. J. Quantum Chem. 1987.31, 269. (16) Press, W. H.; et al. in ref 10, Chapter 4. (17) Wolfram, S.J. Mathemcrtica. A System for Doing Mathematics by Computer, 2nd ed.; Addison-Wesley: Redwood City, CA, 1991.
Anisotropic Molecular Polarizabilities, Dipole Moments, and Quadrupole Moments of (CH2)2X, (CH9)2X,and C4H4X(X = 0, S, Se). Comparison of Experimental Results and ab Initio Calculations Michael H. Coonan, Ian E. Craven, Mark R. Hesling, Geoffrey L. D. Ritchie,*,' and Mark A. Spackman Department of Chemistry, University of New England, Armidale, New South Wales 2351, Australia (Received: March 11, 1992)
Measurements of the temperature dependence of the Cotton-Mouton effects, together with the Rayleigh depolarization ratios and mean polarizabilities,of species in the series (CH2)2X,(CH3)2X, and C4H4X (X = 0, S) are reported and analyzed to yield optical-frequency (632.8 nm) principal polarizabilities of the relevant molecules. In addition, the previously developed 6-31G (+sd+sp) basis set is used to obtain SCF and MP2-level mean and anisotropic polarizabilities and also dipole and quadrupole moments for all nine molecules in the series (CH2)2X,(CHJ2X, and C4HJ (X = 0,S,Se). The results elucidate the electronic charge distributions in these molecules.
Introduction The determination, by experimental and theoretical methods, of accurate molecular polarizabilities and multipole moments continues to be of interest, not only because of the intrinsic value of these properties but also because their ab initio calculation is an important test of the methods of computational quantum chemistry. However, relatively few studies, especially interactive experimental and theoretical studies, have been reported of extended series of molecules, from which useful structural generalizations can be drawn. In this investigation, measurements of the temperature dependence of the Cotton-Mouton effect (magnetic-field-induced birefringen~e)~-~ were combined with Rayleigh depolarization ratios and mean polarizabilities to derive the optical-frequency prinoipal polarizabilities of molecules in the series (CH2)2X,(CH3)2X,and C4H4Xwith X = 0 and S. The results were complemented and extended by SCF and MP2-level calculations, with the previously developed 6-3 l G (+sd+sp) basis set? of the mean and anisotropicoptical-frequency polarizabilities, and also the electric dipole and quadrupole moments, of all nine species in the series (CH2)2X,(CH3)zX,and C4H4Xwith X = 0, S,and Se (Le., oxirane, thiirane selenirane; dimethyl ether, sulfide, and selenide; furan, thiophene, and selenophene). Analysis of the data elucidates the electronic charge distributions in these molecules. 0022-365419212096-7301$03.00/0
Theory
The low-density molar Cotton-Mouton constant, ,C, of a diamagnetic species is, in SI units6 ,C = 2nV,&3(n2
+ 2)2]-1[(nll- nl)K2]g=o
(1)
= ( N A ~ 0 2 / 2 7 0 ~ o )+l A(W-' ~ [ 4 x X X- X) + olyy(xJy-
x) + ~ z z ~ x -XI11 z z
(2)
where eq 1 is a measure of the refractive index difference, nilnL, induced in the gas by the magnetic induction B, and eq 2 is the theoretical relationship3between the observed birefringence and fundamental molecular properties. In eq 2, Aq (=qafl,afl 1 /3qee,flfl) is the anisotropy in the magnetic hyperpolarizability; am,aw,a,,and xu, xyy,xzare, for the molecules of interest here, the principal elements of the optical-frequency polarizability and the magnetizability, respectively; and x (=1/3xea)is the mean magnetizability. Obviously, eq 2 has the form ,C = P + Q T l (3) so that the Cotton-Mouton constant of an anisotropic molecule should exhibit a linear dependence on the reciprocal of the temperature. If the effect is measured over a range of temperature such that a plot of ,C against T'can be reliably extrapolated to T I = 0, the intercept P, given by 0 1992 American Chemical Society
