Partition theory: A very simple illustration
aa r X i v : . [ c ond - m a t . o t h e r] J un Partition theory: A very simple illustration
Morrel H. Cohen
Department of Physics and Astronomy, Rutgers University,126 Frelinghuysen Rd., Piscataway, NJ 08854, USA andDepartment of Chemistry, Princeton University, Washington Rd., Princeton, NJ 08544, USA
Adam Wasserman
Department of Chemistry and Chemical Biology,Harvard University, 12 Oxford St., Cambridge MA 02138, USA
Kieron Burke
Department of Chemistry, University of California at Irvine,1102 Natural Sciences 2, Irvine, CA 92697, USA
We illustrate the main features of a recently proposed method based on ensemble density func-tional theory to divide rigorously a complex molecular system into its parts [M.H. Cohen and A.Wasserman, J. Phys. Chem. A , 2229 (2007)]. The illustrative system is an analog of thehydrogen molecule for which analytic expressions for the densities of the parts (hydrogen “atoms”)are found along with the “reactivity potential” that enters the theory. While previous formulationsof Chemical Reactivity Theory lead to zero, or undefined, values for the chemical hardness of theisolated parts, we demonstrate they can acquire a finite and positive hardness within the presentformulation.
1. INTRODUCTION
In a series of recent papers [1–3], two of us have devel-oped a rigorous method for dividing a complex systeminto its parts based on density-functional theory [4–8].The underlying theory, partition-theory (PT), was usedto construct a formulation of chemical reactivity theory(CRT) [3] which, for the first time, is consistent with theunderlying density-functional theory [8, 9] and is richerin structure than the preexisting CRT [10–13].In PT [1–3], a sharp definition of the individual partsinto which the whole system is partitioned is achievedfirst by selecting the nuclei of each putative part andmaintaining these in the positions in which they occur inthe whole and then requiring that the sum of the electrondensities of the parts, each of which is treated as thoughisolated, add up exactly to the electron density of thewhole (the density constraint). The electron densities ofthe parts are then to be determined by minimizing thesum of the density functionals of the individual parts withrespect to the densities of the parts subject to the densityconstraint. The density functional used, that of ref.[8](PPLB), allows for the existence of noninteger numbers ofelectrons on each part, necessary e.g. for the definitionsof electronegativity [12] and hardness [13], key indicesof chemical reactivity [3], and for incorporating covalentbonding between inequivalent parts.The minimization proceeds via a Legendre transforma-tion, which introduces a reactivity potential v R ( r ) as theLagrange multiplier of the density constraint. Thus, theformalism can become computationally complex. Firstthe electron density of the whole system must be deter-mined. Then, the densities of the parts must be deter-mined simultaneously with v R , all of which is requiredto set the stage for the determination of mutual reactivi- ties between parts, though certain self-reactivities can bedetermined for each species alone without reference to alarger system [3].Accordingly, in the present paper, we develope the par-tition theory in detail for an extremely simple systemto exhibit its main features explicitly. The illustrativesystem is an analog of the hydrogen molecule in whichthe electrons move in one dimension along the molecularaxis without interacting, and the nuclear Coulomb poten-tials are replaced by attractive delta-function potentials.As a consequence of these extreme simplifications, manyquantities of interest can be determined analytically in atransparent manner, including the electron density of themolecule, of its parts (the “atoms”), and the reactivitypotential at all internuclear separations.In Section 2, the model is defined and the moleculardensity obtained. In Section 3, the parts are defined,shown to have one electron each, and a polar representa-tion for their wave functions found which facilitates theminimization. In Section 4, the minimization is carriedout, resulting in an Euler equation for the polar angle β ( x ) of that representation. β ( x ) is found in Section 5and used to determine the reactivity potential v R in Sec-tion 6. The principle of electronegativity equalizationformulated in refs.[2] and [3] is shown to hold in Section7. Also in Section 7, the hardness [3] of the isolated Hatom is calculated, shown to be nonzero, and correlatedwith the strength with which its electron is bound. Thus,despite the fact that the model is a caricature of the realsystem, meaningful features of the partition theory areindeed illustrated by it, as discussed in the concludingSection, 8.
2. 1D-H ; INDEPENDENT ELECTRONSMOVING IN ATTRACTIVE δ -FUNCTIONPOTENTIALS IN ONE DIMENSION Our task is to partition an analog of the H moleculein which two electrons move independently in δ -functionnuclear potentials in one dimension into parts, analogsof H atoms. Each H atom has, by symmetry, only oneelectron, so the need for the PPLB density functionalis avoided. Indeed no explicit use of density-functionaltheory is required for either the molecule or the atoms.The ground-state wave function ψ and energy E of anisolated H atom are (atomic units are used throughout): ψ ( x ) = √ Ze − Z | x | , (2.1) E = − Z / . (2.2)In Eq.(2.1), ( − Z ) is the strength of the δ -function po-tential. To draw the analogy closer to real hydrogenicatoms, one could equate Z to the nuclear charge.The ground-state energy E ( N = 1) of one electronmoving independently in the two δ -function potentialscentered at x = ± a is E ( N = 1) = − κ /
2, where κ satisfies κ = 2 Z/ (1 + tanh κa ) . (2.3)The corresponding wavefunction is: ψ M ( x ) = Be κ ( a −| x | ) , | x | > a = B cosh κx cosh κa , | x | < a (cid:27) , (2.4)where B = κ / (cid:20) κa cosh κa + tanh κa (cid:21) − / ; (2.5)Note that κ → Z as a → κ → Z as a → ∞ (separated atom limit).The two-electron molecular electron density is givenby: n M ( x ) = 2 | ψ M ( x ) | , (2.6)and the total energy of the molecule is E M ( N = 2) = 2 E M ( N = 1) = − κ , (2.7)where N is the number of electrons in the molecule. Thechemical potential of the molecule is therefore µ M = E (2) − E (1) = E (1) = − κ / . (2.8)
3. PARITY DECOMPOSITION
We now partition the molecule into two parts α =1 ,
2, each having a real one-electron wave function ψ α ,localized around − a and + a respectively, so that n M ( x )is given by n M ( x ) = n ( x ) + n ( x ) , (3.1) where n α ( x ) is the electron density of each part α = 1 , ψ α ( x ) = p n α ( x ) . (3.2)They are mirror images of each other, ψ ( x ) = ψ ( − x ) , (3.3)and both are normalized.We now decompose the ψ α into their symmetric, ψ s ( − x ) = ψ s ( x ), and antisymmetric, ψ a ( − x ) = − ψ a ( x ),parts by a rotation within the function space they span, ψ = 1 √ ψ s + ψ a ) , ψ = 1 √ ψ s − ψ a ) ; (3.4) ψ s = 1 √ ψ + ψ ) , ψ a = 1 √ ψ − ψ ) . (3.5)The rotation leaves “lengths” within the space invariantso that n M = ψ s + ψ a . (3.6)We next introduce β = β ( x ), a polar angle in the functionspace, ψ s = √ n M cos β , ψ a = √ n M sin β , (3.7)so that ψ , = p n M / β ± sin β ) (3.8)Because the ψ α are non-negative, | β | cannot exceed π/ β must be an odd function of x , to ensure ψ α is also odd. This also guarantees normalization of ψ α .
4. THE EULER EQUATION FOR β ( x ) To apply PT [2, 3], begin with the original Hamiltonian H = 12 X i =1 , p i − Z X i =1 , [ δ ( x i − a ) + δ ( x i + a )] . (4.1)Then divide the system into overlapping regions, eachwith a given number of electrons. In this case, we chooseone electron on the left, and the other on the right. Thuswe have two 1-electron problems: H α = p v α , v , = − Zδ ( x ∓ a ) . (4.2)The PT problem is to minimize E = ( ψ , H ψ ) + ( ψ , H ψ ) , (4.3)subject to normalization of the wavefunctions, but alsoto the constraint that the total density equal the origi-nal molecular density, Eq.(3.1). (Without the latter con-straint, we’d obviously find ψ , = ψ ( x = ∓ a )). In thepolar representation of Sec.3, both density and normal-ization constraints are automatically satisfied, so the par-tition problem becomes simply minimizing E as a func-tional of β . That functional is E = Z dx (cid:26) (cid:20) n ′ M n M − n ′′ M + n M ( β ′ ) (cid:21) + 12 n M [( v + v ) + ( v − v ) sin 2 β ] (cid:27) . (4.4)Varying it yields δ E = Z dx { n M β ′ δβ ′ + ( v − v ) n M cos 2 βδβ } . (4.5)Integrating by parts, as usual, leads to δ E = 2 nβ ′ δβ | x =+ ∞ x = −∞ + Z dx (cid:26) ddx (cid:18) n M dβdx (cid:19) + ( v − v ) n M cos 2 β } δβ . (4.6)For E to be stationary with respect to arbitrary varia-tions δβ of β , both terms contributing to δ E in Eq.(4.6)must vanish. The Euler equation which results from thevanishing of the second term in Eq.(4.6) is − ddx (cid:18) n M dβdx (cid:19) + ( v − v ) n M cos 2 β = 0 , (4.7) ddx (cid:18) n M dβdx (cid:19) + Z ( δ ( x − a ) − δ ( x + a )) ×× n M cos 2 β = 0 . (4.8)The vanishing of the first term in Eq.(4.6) sets the bound-ary condition at infinity on the Euler equation (4.8).There are two possibilities, the vanishing of β ′ at infinityor the fixing of β there so that δβ must vanish. As weshall see in Section 5, imposing the latter results in anunacceptable divergence in β ′ at infinity. We thereforeimpose the boundary condition β ′ ( x ) = 0 , | x | = ∞ . (4.9)
5. SOLVING FOR β ( x ) Eq.(4.8) becomes ddx (cid:18) n M dβdx (cid:19) = 0 , | x | 6 = a , (5.1)subject to the boundary conditions Eq.(4.9) and β ( a − ) = β ( a + ) ≡ β a β ′ ( a − ) − β ′ ( a + ) = Z cos 2 β a (cid:27) x = a , (5.2) β ( − a + ) = β ( − a − ) = − β a β ′ ( − a − ) − β ′ ( − a + ) = Z cos 2 β a (cid:27) x = − a . (5.3) The general solution of (5.1) is dβ ( x ) dx = c n M ( x ) , (5.4) β ( x ) = Z x dx ′ c n M ( x ) + c . (5.5)where c and c are constants. As implied above in Sec-tion 4, if c does not vanish β ′ diverges exponentially atinfinity, according to Eq.(5.4), because n M goes exponen-tially to zero, so, in accordance with Eq.(4.9), c vanishesfor | x | > a , and β ( x ) is constant there, β ( x ) = β a , x > a = − β a , x < − a (cid:27) . (5.6)For | x | < a we can rewrite Eq.(5.6) as β ( x ) = Z x − a dx ′ c n M ( x ′ ) − β a , (5.7)which implies that β a = 12 Z a − a dx c n M ( x ) . (5.8)From (5.8) we can relate c to β a via Eq.(2.6), c = 2 κB β a cosh κa tanh κa . (5.9)Inserting (5.9) for c into Eq.(5.4) and the result into theBC (5.2) or (5.3) produces an equation for β a , β a = Z κ sinh 2 κa cos 2 β a . (5.10)Inserting Eqs.(5.9) and (2.6) into Eq.(5.7) yields the re-markably simple result β ( x ) = tanh κx tanh κa β a , < | x | < a . (5.11)Eqs.(5.6), (5.10), and (5.11), together with Eq.(2.3)provide a complete analytic solution for β ( x ) and throughEqs.(3.5) and (3.8) for the ψ α . In Fig.1 we show n M , n and n vs. x for Z = 1 and a = 1. We see thateach localized density spreads into the neighboring re-gion, and looks quite similar to an atomic density. Tosee the differences from isolated atomic orbitals, in Fig.2we make the distance smaller ( a = 0 . ψ ( x ) (solid line) and com-pare it with the pure exponential orbital ψ ( x ) of Eq.2.1(dashed line). The orbital ψ resembles ψ and tends toit for large a , but is distorted with respect to it for small a . Its maximum is still a cusp at x = a , but it also showsa second cusp at x = − a . Since κ > Z always (Eq.(2.3)),and either ψ or ψ is proportional to ψ M for | x | > a ,where β = β a is constant, the PT atomic densities andorbitals decay more rapidly than isolated atoms. Since n M ( x ) −1 1 x FIG. 1: Molecular density n M ( x ) (solid), and “atomic” den-sities n ( x ) and n ( x ) (dotted) for Z = 1 and a = 1. ) x ( ψ x FIG. 2: Right-side “atomic” orbital ψ ( x ) (solid) and pureexponential orbital ψ ( x ) (dashed) for Z = 1 and a = 0 . their normalization is the same, this in turn means en-hanced density between the ‘nuclei’, due to bonding. InFig.3, we show β ( x ) for Z = 1, and a = 0 .
1, 1, and 10.Qualitatively, from Eq.(5.11), β ( x ) ≃ β a xa , x < min(1 /κ, a )= β a , x > min(1 /κ, a )and if Za >> β a ≃ π/ Za << β a ≃ a . The interpretationof these results is given in terms of (3.8), outside the bondregion. If β a is small, both ‘atoms’ share the density ineach outside region. But if β a is close to π/
4, each atomdominates on its own side, consuming the entire densitythere. π /4 π /4 − x/a a aa =10 =0.1=1 β ( x ) −1 0 1 −1 0 1 FIG. 3: β ( x ) vs. x , as given by Eq.(5.11), for fixed Z = 1 and3 different values of a .
6. THE REACTIVITY POTENTIAL
The one-electron wave functions ψ ( x ) and ψ ( x ) arenot eigenstates of the part-Hamiltonians H and H ofEq.(4.2). The natural question arises: What are theyeigenstates of? The partition theory of refs. [1–3] dic-tates that they are eigenstates of the modified single-electron Hamiltonians H Rα = p / V α , α = 1 , (cid:18) p V α (cid:19) ψ α = µ M ψ α , α = 1 , . (6.1) V α = v α + v R (6.2)where the eigenvalue, regardless of the part α , is preciselyequal to the molecular chemical potential µ M of Eq.(2.8).The potential v R ( x ) is the reactivity potential that wenow construct explicitly. Summing over α and dividingby ψ + ψ yields a symmetric expression for v R , v R = µ M − ψ + ψ p ψ + ψ ) − v ψ + v ψ ψ + ψ . (6.3) ψ and ψ can be reexpressed in terms of ψ s and ψ a ,Eq.(3.5). Noting that n M = 2 ψ M , (6.4)using Eq.(3.7) for ψ s,a , and taking the δ -function char-acter of v α into account results in v R = µ M + 12 ψ M cos β d dx ( ψ M cos β ) −
12 ( v + v )(1 + tan β a ) . (6.5)The molecular wave function ψ M satisfies theSchr¨odinger equation, − d ψ M dx + ( v + v ) ψ M = µ M ψ M , (6.6)which can be used to transform Eq.(6.5) to v R = − ((cid:20) ψ M dψ M dx dβdx + d βdx (cid:21) tan β + (cid:18) dβdx (cid:19) ) + 12 ( v + v )(1 − tan β a ) . (6.7)Using Eq.(6.4), the Schr¨odinger-like equation for β ,Eq.(4.8), can be rewritten as − (cid:20) ψ M dψ M dx dβdx + d βdx (cid:21) + 12 ( v − v ) cos 2 β a = 0 . (6.8)Multiplying Eq.(6.8) by tan β , invoking the oddness of β and the δ -functions in v and v , and subtracting theresult from Eq.(6.7) yields for v R v R = − (cid:18) dβdx (cid:19) + 12 ( v + v ) [1 − (1 + cos 2 β a ) tan β a ] . (6.9) x / ax / ax / a v R ( x ) v R ( x ) v R ( x ) −0.6−0.3 0 0−0.4−0.2 0 −1 0 1−0.4−0.2 0 0−1 1−1 1 FIG. 4: Reactivity potential v R , Eq.(6.10) for fixed Z = 1 and3 different values of a : a = 0 . a = 1 (middle)and a = 10 (bottom). The δ -functions at ± a are indicated byarrows. Inserting our previous result for β ( x ), Eqs.(5.6) and(5.11) into (6.9) yields an explicit result for v R , v R = µ M β a tanh κa θ ( a − | x | )cosh κx + 12 ( v + v ) [1 − sin 2 β a ] , (6.10)where θ ( y ) = 0 for y <
0, 1 for y > v R ( x ) vanishes for | x | > a ,has attractive δ -functions at ± a whose weights increasemonotonically from 0 to Z as Za decreases from infinityto zero, and has an attractive inverse cosh ( x ) componentfor | x | < a . For the united atom case, Za ↓ v + v R = v + v R = 2 v simply reproduces the molecularpotential, and ψ = ψ = ψ M as they should. Figure4 displays v R vs. x for fixed Z = 1 and representativevalues of a . The reactivity potential is almost flat forsmall separations, a wide well in between the two atomsfor intermediate separations, and a narrow well that is farfrom both atoms at large separations. Figure 5 displaysthe weights of the δ -function components of v R dividedby Z vs. a .As shown in ref.[3], the Kohn-Sham (KS) HOMOeigenvalue of each part must be identical to the chem-ical potential of the whole in the added presence of v R .In our simple example, the KS potential of a part re-duces to the nuclear δ -function potential of one H atom.Adding v R to the nuclear potential must therefore trans-form the HOMO energy E , Eq.(2.2), of the isolatedatom to the more negative HOMO energy of the molecule E ( N = 1) = − κ /
2, which is its chemical potential(Eq.(2.8)). v R must be attractive to do that, which itis, from Eqs.(6.9) and (6.10). In our simple example, v R makes the delta function of the atom more negative,adds the attractive inverse cosh potential between the a δ a π /4 β ( a ) − f un c ti on s t r e ng t h −0.5−0.4−0.3−0.2−0.1 0 0 1 2 3 4 5 FIG. 5: Weights of the δ -function components of v R dividedby Z as a function of a for fixed Z = 1, from the second termof Eq.(6.10). The inset shows β a vs. a E a unoccupiedoccupied state −2−1.5−1−0.5 0 0 1 2 3 4 5 6
FIG. 6: Energy as a function of a , in atomic units, for thetwo lowest-energy solutions of Eq.(6.1). Z = 1 for this plot. atoms, and adds an attractive ghost delta function at theposition of the other atom to force the wave function todecay sufficiently rapidly outside the molecule.In the limit of infinite separation v + v R reduces to v and v + v R reduces to v , except for | x | < a , where theattractive potential v R ( x ) = π E
16 1cosh Zx , | x | < a , (6.11)persists. This potential has at least one additional boundstate, but with binding energy less than | E | . Thus itis unoccupied, and does not affect our results. The a -dependence of this state’s energy is shown for fixed Z inFig.6. For very large separation between the atoms, it islocalized at the center of the inverse cosh ( x ) componentof v R , but it rapidly delocalizes for smaller separations.In particular, for Z = 1, it is highly delocalized when a < ∼ .
4, where it vanishes into the continuum.
7. SUSCEPTIBILITY AND HARDNESS
Having found the reactivity potential, we now illus-trate the construction of reactivity indices. In the CRTof ref.[3], each part α is represented by an ensemble ofPPLB type containing contributions with only two in-teger electron numbers, p α and p α +1 . The principle ofelectronegativity equalization is expressed as the equal-ity of the chemical potential of each part in the presenceof the reactivity potential, µ Rα , to the chemical potentialof the molecule, µ M , µ Rα = µ M , ∀ α . (7.1)The µ Rα are defined as the difference between the groundstate energies of α for p α + 1 and p α electrons in thepresence of v R , µ Ra = E Rα ( p α + 1) − E Rα ( p α ) , (7.2)and similarly for µ M µ M = E M ( N M ) − E M ( N M − . (7.3)In our simple example, µ M is given in Eq.(2.8). Therelevant value of p α is zero, so that µ Rα is just E Rα (1), thelowest eigenvalue of H Rα = H α + v R , (7.4)with H α given by Eq.(4.2) and v R by Eq.(6.10). The ex-plicit construction of v R in Section 6, not possible in gen-eral, guarantees that Eq.(7.1) and therefore electroneg-ativity equalization holds. In the general case, a modi-fication of the Car-Parrinello scheme [14, 15] guaranteeselectronegativity equalization.The susceptibility of part α measures the response ofthe density of part α to a small change in the potential V α of Eq.(6.2): χ α ( x, x ′ ) = − δn α ( x ) δ V α ( x ′ ) . (7.5)For 2 electrons, it is simple to show that χ α ( x, x ′ ) = − ψ α ( x ) G α ( µ M ; x, x ′ ) ψ α ( x ′ ) , (7.6)where G α ( µ M ; x, x ′ ) is given by the E → µ M limit of: G α ( E ; x, x ′ ) = G α ( E ; x, x ′ ) − ψ a ( x ) ψ ( x ′ ) E − µ M , (7.7)and G α is the Green’s function for part α : G α ( E ; x, x ′ ) = (cid:20) E − (cid:18) p V α (cid:19)(cid:21) − ( x, x ′ ) . (7.8)Figure 7 shows the susceptibility of the right “atom” forvarious interatomic separations when the perturbing po-tential is added at x = 3 (the numerical calculationswere done as described in the Appendix). Electrons flowaway from x , building up a peak at x (positive becauseof the minus sign in the definition of χ α , Eq.(7.5)), anda negative peak at the closest maximum of the chargedensity, i.e. at a . With the analytic Green function of ( , ) x χ a =1 a =2 a =2.9 3.1 a =4 x x −0.04 0 0.04−2 0 2 4 6 −0.15 0 0.1 0.2−2 0 2 4 6−0.2 0 0.4 0.8−2 0 2 4 6 −0.15 0 0.1 0.2−2 0 2 4 6 FIG. 7: Susceptibility χ ( x , x ) of the right “atom” obtainedfrom Eqs.(7.6)-(7.8), as indicated in the Appendix, when x is set to 3 a.u. Each panel corresponds to a different valueof the internuclear distance, a . The lower-left panel shows χ ( x , x ) when a is just below (solid) and just above (dotted) x . an isolated “atom”[20] and Eqs.(7.6)-(7.7), χ α can be ob-tained analytically in the large-separation limit: χ α ( x, x ′ ) = 2 e − Z | x | (cid:26) e − Z | x − x ′ | − (cid:20)
12 + Z ( | x | + | x ′ | ) (cid:21) × e − Z ( | x | + | x ′ | ) o e − Z | x ′ | (7.9)We now construct the susceptibility of the whole sys-tem, χ R , by adding together the susceptibilities of theparts, χ R ( x, x ′ ) = X α χ α ( x, x ′ ) . (7.10)The inverse of χ R determines the hardness matrix η αβ asshown in refs.[2] and [3]: η αβ = Z Z dxdx ′ f α ( x ) χ − R ( x, x ′ ) f β ( x ′ ) , (7.11)where the Fukui function of part α , f α ( x ), f α ( x ) = dn α ( N a , x ) dN a , (7.12)is simply equal to ψ α ( x ) for 2 non-interacting electrons,since n α ( N α , x ) = N α ψ α ( x ) (see also ref.[16]). Thus, wehave η αβ = Z Z dxdx ′ ψ α ( x ) χ − R ( x, x ′ ) ψ β ( x ′ ) . (7.13)Figure 8 shows the self-hardness η αα for an isolated H-“atom”, as a function of Z . The constancy of the hard-ness for large Z can be understood qualitatively as fol-lows. The inverse susceptibility has units of energy timeslength squared. When Z is large, it establishes a lengthscale inversely proportional to Z , and an energy scaleproportional to Z , so the Z -dependence cancels out inthe inverse susceptibility. To obtain the hardness, wemultiply χ − R on the left and right by the Fukui function,which has the dimension of inverse length. Integratingover position on the left and right then cancels out the Z -dependence arising from the Fukui functions, and theresult is a Z -independent hardness.
8. CONCLUSIONS
Despite the extreme simplicity of the 1D-H2 modelanalyzed here – two non-interacting electrons movingin 1D under the influcence of two equivalent attractivedelta-function potentials – that model allows us to illus-trate the essential features of our partition theory and ofkey indices of our chemical reactivity via straightforwardanalysis and easy computations.We have shown that the electron density of themolecule can be decomposed exactly into a sum of atomicdensities, a rigorous solution of the “atoms-in-molecules”problem [17].Electronegativity equalization [18] is built into the par-tition by the symmetry of the problem, so this homonu-clear model does not illustrate that principle as well asa heteronuclear model would. Nevertheless, the currentexample does illustrate a key feature of the new CRT, thechemical context dependence of the reactivity indices,in this case the electronegativity of a part, introducedthrough the presence of v R in the Schr¨odinger equationfor ψ α , cf. Eq.(7.4). It also demonstrates that the reac-tivity potential remains finite as two atoms separate, buthas no effect on the partitioning after separation.Another serious shortcoming of the earlier formulationsof DFT-based CRT is the vanishing of the hardness. Wehave shown explicitly here that the self-hardness, as de-fined in [3], of an isolated “atom” is positive. Interest-ingly, the hardness saturates as the ionization energy ofthe “atom” increases, raising the very interesting ques-tion of whether such a saturation of hardness with ion-ization energy exists in real systems. For this model,a strong positive correlation between hardness and ion-ization energy exists only over the limited range of Z between 0.4 and 0.7.KB is supported by NSF CHE-0355405. APPENDIX: Numerical calculation of thesusceptibility
We first obtained G α ( E ; x, x ′ ) according to the well-known prescription [19]: G α ( E ; x, x ′ ) = 2 ψ α,L ( E, x < ) ψ α,R ( E, x > ) W [ ψ α,L , ψ α,R ] , (A.1) η Z FIG. 8: Self-hardness vs. Z in the separated-atom limit(atomic units). where x < = inf( x, x ′ ), x > = sup(x , x ′ ) W [ ψ α,L , ψ α,R ] = ψ α,L ( E, x ) ψ ′ α,R ( E, x ) − ψ ′ α,L ( E, x ) ψ α,R ( E, x ) , (A.2)and the orbitals ψ α,L and ψ α,R are solutions of (cid:20) p V α ( x ) (cid:21) ψ α,L,R ( E, x ) = Eψ α,L,R ( E, x ) (A.3)satisfying left and right-boundary conditions, respec-tively: | ψ α,L ( E, x ) | ↓ , x ↓ −∞ (A.4) | ψ α,R ( E, x ) | ↓ , x ↑ ∞ (A.5)The potential V α ( x ) of Eq.(A.3) is given by Eq.(6.2), withthe reactivity potential v R ( x ) of Eq.(6.10). The compu-tations of ψ α,L,R ( E, x ) were carried out at E = µ M ± ∆ E with ∆ E chosen for numerical convenience, i.e. largeenough so that sup x,x ′ | G α ( µ M ± ∆ E ) | does not becomeso large as to be inconvenient on the one hand, and smallenough so that [ G α ( µ M + ∆ E ) + G α ( µ M − ∆ E )] doesnot differ significantly from its limit at ∆ E ↓
0. We thencalculated G α of Eq.(7.7) as: G α ( µ M ; x, x ′ ) = 12 [ G α ( µ M + ∆ E ; x, x ′ ) + G α ( µ M − ∆ E ; x, x ′ )](A.6) [1] M.H. Cohen and A. Wasserman, Israel J. Chem. , 219(2003).[2] M.H. Cohen and A. Wasserman, J. Stat. Phys. , 1125(2006).[3] M.H. Cohen and A. Wasserman, J. Phys. Chem. A ,2229 (2007).[4] P. Hohenberg and W. Kohn, Phys. Rev. , 864(1964).[5] W. Kohn and L.J. Sham, Phys. Rev. , A1133 (1965).[6] M. Levy, Proc. Nat. Acad. Sci. USA , 6062 (1979).[7] E.H. Lieb, in Physics as Natural Philosophy , eds. A. Shi-mony, and H. Feshbach, MIT Press, Cambridge, p.111(1982).[8] J.P. Perdew, R.G. Parr, M. Levy, and J.R. Balduz, Jr.,Phys. Rev. Lett. , 1691 (1982).[9] J.P. Perdew, in Density Functional Methods in Physics ,ed. R.M. Dreizler and J. da Providencia, Plenum, NewYork, p.265 (1985).[10] R. Parr and W. Yang,
Density Functional Theoryof Atoms and Molecules , Oxford University Press,New York (1989).[11] P. Geerlings, F. De Proft, and W. Langenaeker,
Chem. Rev , , 1793 (2003).[12] R.G. Parr, R.A. Donnelly, M. Levy, and W.E. Palke, J.Chem. Phys. . , 3801 (1978).[13] R.G. Parr and R.G. Pearson, J. Am. Chem. Soc. ,7512 (1983).[14] R. Car and M. Parrinello, Phys. Rev. Lett. , 2471(1985).[15] M.H. Cohen and R. Car, unpublished.[16] P. Fuentealba, E. Chamorro, and C. C´ardenas, Int. J.Quant. Chem. , 37 (2007).[17] F.L. Hirshfeld, Theor. Chim. Acta , 129 (1977); R.F. W. Bader, Atoms in Molecules - A Quantum Theory ,Oxford University Press, Oxford, 1990.; R.F. Nalewa-jski and R.G. Parr, Proc. Natl. Acad. Sci. USA , 8879(2000).[18] R.T. Sanderson, Science , 670 (1951).[19] G. Arfken, Mathematical Methods for Physicists , Aca-demic Press, New York (1970).[20] A. Szabo and N.S. Ostlund,