Homonuclear diatomic systems: differential equation for mirror plane electron density, and fragment model

Journal of MolecularStructure (Theochem), 288 (1993)235-243 0166-1280/93/$06.00 0 1993- ElsevierScience Publishers B.V. All rights reserved

235

Homonuclear diatotiic systems: differential equation for mirror plane electron density, and fragment model J.M. L6peza, J.A. Alonsoa, N.H. March*yb aDepartamento de Fisica Tedrica, Universidad de Valladolid, E-4701 1 Valladolid, Spain bTheoretical Chemistry Department, University of Oxford, 5 South Parks Road, Oxford OX1 3UB, UK (Received 25 March 1993; accepted 10 April 1993) Abstract

The mirror plane electron density is an important characteristicof molecularbonding in homonucleardiatomic systems.Therefore, by means of the concept of the Pauli potential an ordinary differentialequation is constructed for the mirror plane density amplitude. The effective potential energy entering this Schriidinger equation is determined exactly for Hz and approximately, by a low-order density gradient form of density functional theory, for the supermolecularion (I&. Finally, a fragment model for representing the mirror plane density is proposed.

1. Introduction

In an earlier study by Alonso et al. [I] the mirror plane electron density has been shown to characterize the ground-state potential energy curve in the hydrogen molecule ion Hz. In this case, it was also possible by exploiting the separability of the ground-state wavefunction in confocal elliptic coordinates to set up an exact ordinary (rather than partial) differential equation for the square root of the mirror plane density, i.e. the electron density amplitude. Such separability is not, in general, feasible for a multielectron homonuclear diatomic system, which is the main focus of the present work. Therefore, starting from general well-established ideas involving the so-called Pauli potential [2-41, the Schriidinger (partial differential) equation for the ground-state density amplitude [ p(r; R)]‘12 in a homonuclear diatomic at internuclear separation R is, at first formally, reduced to an equation for * Corresponding

author.

the mirror plane density amplitude. The way this can be done is set out in Section 2 below and in the Appendix. The differential equation for the ground state density amplitude is characterized by an effective potential energy, and this is constructed exactly for Hi in Section 2.2, thereby relating to our earlier study on this system [l]. An approximate form for this potential energy is displayed in Section 3 for the supermolecular ion (K&)2 by means of a low-order density gradient form of density functional theory. Finally, a representation of the mirror plane density by a “fragment model” is proposed in Section 4 and tested in the (K&)2 supermolecular ion. 2. Differential equation for the mirror plane density amplitude

Suppose we have a non-interacting system of electrons moving in a potential energy V(r). Then the introduction of the Pauli potential V,,(r) bypasses the need for constructing the one-electron wavefunctions [2] generated by V(r) and allows

236

J.M. Ldpez et al./J. Mol. Struct. (Theo&m)

one to calculate the electron density amplitude [p(r)] ‘I2 G A(r) directly from the Schrodingerlike equation

We next assert that As(s), by analogy with Eq. (l), satisfies a Schrddinger equation

d2A&) + a2A&) 3X2

V2A(r) +$

[E- V(r) - VJr)]A(r)

= 0

6Tw

(2)

6p(r)

where Ts[p] is the (as yet unknown) single-particle kinetic energy functional and Tw[ p] is the Weizsacker [6] inhomogeneity kinetic energy Tw[p] = &/ydr

ay2

2m + F

ko - Us;

R)IAo(s)

(1)

Vr,(r) can be formally written, following Holas and March [5] as 6Ts --vp(r) = Sp( r)

288 (1993) 235-243

(3)

=0

(5) Given Eq. (5), the central issue in determining the mirror plane density amplitude Ao(s) is to find the effective potential energy Y(s; R). This quantity V(s; R) can be written down by comparing directly Eqs. (5) and (1); however, because this is somewhat formal and the details are a little complex, the argument is given in the Appendix. Below, we immediately take two examples by which we can demonstrate graphically that V(s; R) in Eq. (5) has a rather simple spatial dependence (s-dependence). 2.2. Exact form of l”(s; R) for the hydrogen molecule ion Hz+

2.1. Formal reduction of Eq. (1) to an ordinary dlflerential equation for the mirror plane density amplitude in a homonuclear diatomic system

Let us first work in Cartesian coordinates with the z axis along the bond, and the origin at the midpoint z = 0. Then the plane perpendicularly bisecting the line joining the two identical nuclei - the mirror plane - has evidently electron density p(x, y, 0; R) at z = 0, say ps(x, y; R). But now, because of axial symmetry in the ground state, x and y always enter through the combination s=(x2+y 2) ‘I2 . Thus the mirror plane electron density amplitude [po(x, y; R)]‘/2 will be written as simply A&), the R dependence no longer being indicated explicitly for notational compactness. Clearly, in the second term on the left-hand side of the Schrodinger equation (Eq. (l)), one can put z = 0 everywhere, and one then has

We next make contact with our earlier study of Hz in which, as noted above, the ground-state wavefunction was separated in confocal elliptic coordinates 5 and n. On the mirror plane n = 0 and in the notation of Pauling and Wilson [7] X(S) = Ao(4

where Ao(s) is the mirror plane density amplitude at a perpendicular distance s from the bond axis. But now one has the exact differential equation for X(t) for H,f , namely $

R)]&(s)

= 0

[(p-

I)%]

+ (-xt2+2D5+P)x=o

(7)

where X measures the ground-state energy E(R), D = R/a0 and p= p(X) is given by (see Pauling and Wilson [7])

j++&X2+4X3 8505

~),_)+(q,_._(~),-, +$ [E- Vo(s;R) - I’&;

(6)

- 0.000013X4 - 0.0000028X’ + . . .

Now X(l) = X(2r/R) on the mirror plane, where r = [s2 + (R 2/4)]1’2, i.e. E2 = 1 + 42/R 2. Trans-

forming Eq. (7) to variable s and using Eq. (6),

J.M. Ikpez et al./J. Mol. Struct. (Theochem)

288 (1993) 235-243

231

one is led to the result

“3

jsi+$A:+ (z+$)k,

+[-A(1 +g) +2(1+$+p] (9)

xAo=O

where Ah and Ai stand for dAO/ds and d2Ao/d.s2 respectively. Using s = (x2 +y2)‘12, and d’A&)/dx*+ a2Ao(s)/ay2 = (?32Ao/i2s2) + (l/s)dAo/&, one can rearrange Eq. (9) to the form of Eq. (5), to obtain -51

A:+s

0

Ab

2

s

4

(a.u.)

Fig. 1. Plot of effective potential energy ^Y(s; R) for H: from Eq. (11) for equilibrium separation 4. This is determined to within the eigenvalue co(&).

2S2f$ =o

(10)

This allows one to write an exact explicit form for the effective potential energy Y(s; R) appearing in Eq. (5) for this example of Hg. The result is

(11) For the equilibrium separation &, this effective potential energy -LT(s;4) is plotted in Fig. 1. It is seen to have the shape of a “Coulomb potential with cut-off’. That is, it tends to a constant as s + 0, and it approaches its asymptotic limit for large S, like l/s. What this argument shows is that Eq. (5) can be formulated for Hz, with “Ir(s; R) given explicitly

in Eq. (11). It is to be noted that Eq. (5) must be solved “self-consistently”, both Ao(s) and X appearing in T(s; R). 3. The example of the supermolecular ion (K&J2 One wishes, of course, to use Eq. (5) in a multielectron homonuclear diatomic system. This is done below, within a low-order gradient expansion form of density functional theory (DFT), going back to the work of Weizsacker [6] referred to in connection with the Pauli potential in Eq. (2). The Euler equation of the Weizsacker problem takes the form of an expression for the chemical potential p, which, of course, is the same at any point in the electronic charge cloud of the molecular system under consideration. This equation for p reads (see, for example, Lopez et al. [S])

where 3

3h2 Ck=iG

0 ii

213

(12)

238

J.M. Ldpez et al./J. Mol. Struct. (Theochem) 288 (1993) 235-243 0.0025

1

A

I

R=31.38

a 0 0005 E E 2

0 0000

s

(au.)

Fig. 2. Mirror plane density h(s) for (Kj& from the selfconsistent solution of Eq. (12) corresponding to a low-order gradient expansion form of density functional theory. Curves are for different values of the separation R between the centres of the two K& fragments: 26.38, 31.38 and 33.38 a.u.

Equation (12) has been solved self-consistently for the (Na&)2 supermolecular ion in Ref. 8 (the interested reader should consult that reference for details of the model). The calculations have been extended here to the case of (K$,)2. The mirror plane density p(r, z = 0; R) G p&) for the (K&)2 supermolecular ion is plotted in Fig. 2 for three different values of the distance between the centres of the two K& super-ions R.WeizsZcker’s original proposal corresponds to X = 1 in Eq. (12) (see Eq. (3)); Eq. (12) was also solved in Ref. 8 for X = 0.5, a value of X which has been found to be appropriate for describing most properties of simple metal clusters [9]. The results in Fig. 2 correspond to x = 0.5. It is relevant to recall at this stage that the shape of the curve giving the total energy of the Hz molecular ion as a function of the separation between the two H nuclei is very different from the curve giving the total energy of the supermolecular ion (K&)2 (or (Nai)2) as a function of the separation between the two super-ions K&. Hl is a typical bound molecule, whose binding energy curve has a minimum at the equilibrium inter-

nuclear separation. The dissociation products for very large separation are H + H+. In contrast, however, although the fused super-ion (&0)2’ is stable with respect to the dissociated products K,+,+K &, there is a “fusion” barrier for the reaction K& + K 2’0+ (K~o)~+; in other words the curve giving the energy of the system as a function of the separation between the two K& superions has a maximum (instead of a minimum) at intermediate separations (see Ref. 10 for the case (NaJ2’). The physics of (charged) cluster fragmentation has been discussed elsewhere [8,10] and is not treated here. We merely mention that the smallest interfragment distance used in Fig. 2, (R = 26.38a.u.) corresponds to a situation in which the two K$ fragments are very close to each other, although the united cluster limit (K4,,)2+ has not yet taken place. However, the largest of the three (R = 33.38 a.u.) roughly corresponds to the separation at which the maximum of the barrier occurs. As utilized by March and Murray [3] in early work on DFT, Eq. (12) for the specific case X = 1 becomes a Schriidinger equation for the density amplitude [p(r,R)]'j2 G A(r),which is precisely of the form of Eq. (l), but now with the Pauli 0010 7

0 008

0 006

Fig. 3. sps(s) plotted vs. s from Fig. 2 for the same values of interfragment separation R. The position of the maximum varies with R and can be taken as a “characteristic length” entering the mirror plane density.

J.M. Lbpez et al./J. Mol. Strut.

239

(Theochem) 288 (1993) 235-243

However, knowing ps(s) in Fig. 2, we have employed Eq. (5) directly, in the form

potential given explicitly by VJr; R) = $ck[p(r; R)]2’3

(13)

One can follow the route set out in the Appendix, but now with known V,,(r), and reduce Eq. (12) to an ordinary differential equation for the mirror plane density amplitude A&). This affords, therefore, an alternative route to determining the spatial dependence of p&) shown for (K&)2 at three values of R in Fig. 2. To exhibit a characteristic length in ,&), the quantity s&s) is also plotted in Fig. 3. The movement of the position of the maxima with R is the feature that is noteworthy here. 3.1. Form of effective potential Y(s; R) for (K&)2 Since the full three dimensional Euler equation (Eq. (12)) was solved by the method set out in Ref. 8 it is obviously possible to determine Y(s; R) by a procedure parallel to that set out in the Appendix.

(14) and the left-hand side is plotted in Fig. 4. The general shape of T(s; R) is seen to be quite similar to that in Fig. 1 for Hi. We must note that what is lost in using a “reduced” Schriidinger equation to calculate the mirror plane density amplitude Ao(s) is the absolute magnitude of p,-&) fixed, say, by the bond midpoint density pm. This is, of course, determined by direct use of Eq. (12) in the present example, plus the normalization condition Jp(r; R)dr = N,I, where NC1is the total number of (valence) electrons. This in turn fixes the chemical potential p(R) for a given N,t (see also Fig. 6 below). This matter of normalization can be dealt with, at least in principle, by representing the mirror plane density pa(s) in terms of a fragment model, as set out in Section 4.

4. Spherical “fragment” model constructed to generate the mirror plane density Once one limits attention to the mirror plane density pa(s) alone, its shape can be determined via the density amplitude Ao(s) = [pa(s)]‘/* from Eq. (5) for chosen internuclear separation R; but the absolute value is not thereby known, as has already been noted. Therefore, adopting the long standing chemical idea of “atoms in molecules”, let us write for the two-centre density p(r; R) the approximate form IP(r;R)l”“= -1.0

,II/IIJI,I,IIII/,/II,JII’I”“,“““’II/IIII””

0.00

10.00

20.00

s

30.00

40.00

50

(O.U.)

Fig. 4. Effective potential Y(s; R) to within a constant (compare with Fig. 1 for Hl) required to generate the mirror plane density amplitude A,(s) = [pa(s)]‘/*, with h(s) as in Fig. 2.

J’(R){

[Pr(rl)l”“+[Pf(r2)1’/.) (15)

where pr(r) denotes a spherical “fragment” density, and rl and r2 (as usual) represent distances measured from nuclei 1 and 2. Finally N(R) in Eq. (15) is a normalization factor which depends on n. The choice n = 2 resembles in DFT the

240

J.M. Lbpez et al./J. Mol. Struct. (Theochem) 288 (1993) 235-243

LCAO method based on wavefunctions. Lee et al. [l l] for the H2 molecule have argued for n = 3, which choice we can describe intuitively by noting that it represents the superposition of one-centre Dirac-Slater exchange potentials. The advantage of using Eq. (15) on the mirror plane (r, = r2 = r, say) is now evident. One obtains [p(s,z = 0; R)]%

[p&)]““=

2N(R)[pf(r)]“” (16)

To illustrate Eq. (18), the Weizsiicker plots of the mirror plane density for (K&) in Fig. 2 have been utilized to construct the spherical fragment density pr(r) for r > R/2. The reduced quantity pf(r)/pf(R/2) is plotted in Fig. 5 for the three different values of R shown in Fig. 2, the variable 5 = 2r/R now being employed (see also below). It is important to note that the long-distance behaviour of pf(r) is dominated by the exponential decay

where R2

r=

(19)

112

( ) s2+-

(17)

4

Raising Eq. (16) to the power n then yields PO(S)= [2Jv-(R)]“Pf

(18)

Hence, for a given R, the spatial dependence of the spherical fragment density pf(r) for r > R/2 is precisely determined by p&), independently of the “law of combination” (i.e. the choice of n) of one-centre fragments in Eq. (15).

12 1

0.2

00

R=

,I~(,1111/11~11/11,,11~1”~“/‘~“n~”l””””’I””””

080

100

1 20

1 40

1.60

1.80

2

Fig. 5. Reduced fragment density pf(r)/pf(R/2) extracted using Eq. (18) with p&) taken from Fig. 2. Note that the mirror plane density h(s) only determines the “reduced” fragment density plotted for r > R/2.

where Z,(R) is the molecular ionization potential at separation R. Thus pf(r) at large r contains within it this molecular property Zm:it is not solely characterized by the isolated atoms; rather it reflects “atoms in molecules”, as already emphasized.

4.1. Approximate 3caling”property

offragment

density

Of course, if there was no relation between the fragment densities pf{[s2 + (R2/4)]1/2} in Eq. (18) for different values of R, then there is no gain in economy over the direct use of the mirror plane density Z&S). But the asymptotic form of Eq. (19) at large r now motivates an examination of an approximate “scaling” property of the fragment density. Thus, we rewrite Eq. (19) in terms of the reduced coordinate c = (ri + r2)/R of the confocal elliptic system, but now with < = 2r/ R on the mirror plane rl = r2. The exponential factor becomes then exp [-R(2Z,(R))‘/2<]. But the chemical potential p(R) is related to Z,(R). If we adopt the approximate Mulliken form ]p] w (Z, + EA)/2 and neglect the (small) affinity EA, then one can rewrite Eq. (19) as pf(r) z5 B[exp (-c$)lRm

This asymptotic argument motivates the plot, with f(E) the “scaled fragment” replacing the

241

J.M. Lbpez et al./J. Mol. Struct. (Theochem) 288 (1993) 235-243

asymptotic form exp (-at)

-0.18

in Eq. (20), of

Limit at R=m

(21) 5 d -0.20 :

i

0 .o -0.22

y

6 e -0.23

1

Distance

(a.u)

Fig. 6. Chemical potential p(R) vs. distance R calculated from Eq. (12) for both (K&), and the isolated ion K&. This information is utilized in constructing Fig. 7.

where f(c) has also been “normalized” such that f(t = 1) = 1. Using the chemical potential p(R) for (K&)2 plotted in Fig. 6, f(c) has been constructed using Eq. (21) from the plots of Pti”&(c) in Fig. 5. The spread of curves of f(r) in Fig. 7 for the three values of R is now much reduced: this is the “approximate scaling” we were seeking. The curve labelled R = 26.38 a.u. is, in fact, a case of strong interaction of (K&0) ions, whereas the other cases are more akin to conventional equilibrium distances in quantum chemistry. Using the average of R = 31.38-33.38 a.u. in Fig. 7, we obtain the mirror plane densities shown in Fig. 8. The “fragment” extracted from the larger R values gives a useful approximation to the mirror plane density except for the “strong interaction” case, i.e. for R = 26.38 a.u. This is akin to approaching the united atom limit in conventional

R=26.38

-10

4

4216 0

1

2

3

4 0.8

1 .o

1.2

1.4

16

18

2

F

Fig. 7. Plot of dimensionless fragment density f(E) defined through Eq. (21) for (K$), using the reduced fragment density in Fig. 5, plus chemical potential p(R) from Fig. 6.

Fig. 8. Approximate “recovery” of mirror plane electron density for two largest values of R for (K&)r from average of two large R curves for f(t) in Fig. 7.

242

J.M. Lbpez et al.lJ. Mol. Struct. (Theochem) 288 (1993) 235-243

quantum chemical systems: here any fragment model can be expected to break down. Returning to the point about normalization of the electron density: we must note that whereas an unnormalized “fragment” density pr(r) can be extracted for r > R/2 from the mirror plane density p,-,(s), “atomic” information must be added to construct a normalized pr(r) in to r = 0. Once such an “extrapolation” is made, the “combination law” shown in Eq. (15) can be used to construct an approximate twocentre density which can be normalized for any chosen n. 5. Summary An ordinary differential equation has been set up for the electron density amplitude A&) = ‘j2,where pa(s) is the mirror plane density [PO(f)1 as a function of the perpendicular distance from the bond axis of the homonuclear diatomic system under consideration. The differential equation derived for Ao(s) is characterized by an effective potential energy Y(s : R) and this is constructed exactly for Hl. An approximate form is displayed for Y(s; R) in the supermolecular ion (Klo)z. Finally, the representation of the mirror plane density pa(s) by a spherical “fragment” model is proposed. 6. Acknowledgements This work has been supported by DGICYT (Grant PB 89-0352). One of us (N.H.M.) wishes to acknowledge that his contribution to the present work was carried out during his tenure of the BBV Foundation Chair at the University of Valladolid during the Autumn of 1992.

7. References 1 J.A. Alonso, J.M. Lopez 79 (1993) 1143.

and N.H. March, Mol. Phys.,

N.H. March, Phys. Lett. A, 113, (1985) 66,476. N.H. March and A.M. Murray, Proc. R. Sot. London, Ser. A, 256 (1960) 400. C. Herring and M. Chopra, Phys. Rev., A37 (1988) 31. A. Holas and N.H. March, Phys. Rev., A44 (1991) 5521. C.F. Weizslcker, Z. Phys., 96 (1935) 431. L. Pauling and E.B. Wilson, Introduction to Quantum Mechanics, Dover, New York, 1985. J.M. Lopez, J.A. Alonso, F. Garcias and M. Barranco, Ann. Phys., 1 (1992) 270. 9 Ll. Serra, F. Garcias, M. Barranco, J. Navarro, L.C. Balbas and A. Matlanes, J. Phys., Condens. Matter, 1 (1989) 10391. 10 M. Barranco, J.A. Alonso, F. Garcias and J.M. Lopez, in M. Brenner, T. Liinnroth and F.B. Malik (Eds.), Clustering Phenomena in Atoms and Nuclei, Springer, Berlin, 1992, p. 305. 11 C. Lee, S. Zdravkovic and R.G. Parr, J. Chem. Phys., 92 (1990) 2114.

8.

Appendix: Formal derivation of effective potential energy V(s; R) in Eq. (5) determining the mirror plane density amplitude Ao(4 =

[Pow11’2

To proceed from Eq. (4), one simply adds and subtracts (a2Ao/ax2) + (a2Ao/~y2) to obtain

$2+$+{[(~)z~o-$J] +[($)z=o-%q +(z)z~o} +

$

-

vpo(s; M-40(4

[E- Vo(s;R)

=0

(AlI where we remind the reader that s = (x2 + y2)‘12. One next multiplies and divides the term in braces { } in Eq. (Al) by Ao, in order to combine it with the potential energy terms V. and VP0 already present in Eq. (Al). By comparing the result then obtained with Eq. (5) one obtains the

J.M. Ldpez et al/J. Mol. Struct. (Theochem) 288 (1993) 235-243

rather formal expression to determine the effective potential energy Y(s; R) $J - Y(s; R) = E - Vo(s; R) - V&; R) +~{a,,,=,-Ao,+a,,,=o - AoYY + &,z=o}

(A2)

243

where the notation A,,,_o = (c?~,~/~x~),=,, A oxx = d2Ao/ih2, etc. has been employed. What Eq. (A2) demonstrates is the undoubted existence of the effective potential energy V(s; R) in Eq. (5). In spite of the somewhat complicated nature of Eq. (A2), it is shown by means of the examples in Sections 2 and 3 of the main text that, in fact, the spatial dependence (s-dependence) of “Ir(s; R) is rather simple in these cases.