Subscriber access provided by RMIT University Library
Article
Thermodynamic Justification for the Parabolic Model for Reactivity Indicators with Respect to Electron Number and a Rigorous Definition for the Electrophilicity: The Essential Role Played by the Electronic Entropy Marco Franco-Pérez, José L. Gázquez, Paul Woodson Ayers, and Alberto Vela J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00940 • Publication Date (Web): 21 Dec 2017 Downloaded from http://pubs.acs.org on December 26, 2017
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Thermodynamic Justification for the Parabolic Model for Reactivity Indicators with Respect to Electron Number and a Rigorous Definition for the Electrophilicity: The Essential Role Played by the Electronic Entropy Marco Franco-Pérez,*(a),(b) José L. Gázquez,(a) Paul W. Ayers,(b) and Alberto Vela(c). (a)
Departamento de Química, Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael
Atlixco 186, Ciudad de México, 09340, México. (b)
Department of Chemistry, McMaster University, Hamilton, Ontario, L8S 4M1.
(c)
Departamento de Química, Centro de Investigación y de Estudios Avanzados, Av. Instituto
Politécnico Nacional 2508, Ciudad de México, 07360, México.
*Corresponding author:
[email protected] 1 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 2 of 37
Abstract The temperature-dependence of the Helmholtz free energy with respect to the number of electrons is analyzed within the framework of the Grand Canonical Ensemble. At the zerotemperature limit, the Helmholtz free energy behaves as a Heaviside function of the number of electrons, however, as the temperature increases, the profile smoothens and exhibits a minimum value at non-integer positive values of the fractional electronic charge. We show that the exact average electronic energy as a function of the number of electrons does not display this feature at any temperature, since this behavior is solely due to the electronic entropy. Our mathematical analysis thus indicates that the widely-used parabolic interpolation model should not be viewed as an approximation for the average electronic energy, but for the dependence of the Helmholtz free energy upon the number of electrons, and this analysis is corroborated by numerical results. Finally, an electrophilicity index is defined for the Helmholtz free energy showing that for a given chemical species, there exists a temperature value for which this quantity is equivalent to the electrophilicity index defined within the parabolic interpolation of the electronic energy as a function of the number of electrons. Our formulation suggests that the convexity property of the energy versus the number of electrons together with the entropic contribution does not allow for an analogous nucleophilicity index to be defined. Keywords: density functional theory of chemical reactivity, Helmholtz free energy, grand canonical ensemble, electrophilicity index.
2 ACS Paragon Plus Environment
Page 3 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
1.0 Introduction Chemical reactivity theory based on density functional theory (CR-DFT) has proven to be a very useful tool to analyze several aspects associated to the reactivity of chemical species.1-8 From the CR-DFT perspective, some information regarding the chemical reactivity of the electronic systems can be captured by the chemical reactivity response indexes that quantify how a molecule or an atom is initially perturbed by an approaching reagent. Each of these response indexes has been connected with any of the partial derivatives of the electronic energy with respect to its natural variables, the number of electrons N and the external potential υ (r ) . The first order indexes are thus defined from the following exact differential, 1, 9-11
δE ∂E + dE = dN ∫ δυ ( r ) δυ ( r ) dr ∂N υ ( r ) N
,
(1)
where
∂E ∂N υ (r )
µ=
(2)
δE ρ (r) = δυ ( r ) N
(3)
is the electronic chemical potential and
is the ground state electronic density. In Eq. (1), E is the exact ground state electronic energy functional and thus, the partial or functional derivatives obtained from Eq. (1) at any order, can be considered as true responses of a chemical species in their corresponding
3 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 4 of 37
fundamental electronic states. It is expected that as more coefficients are evaluated for a given chemical species, a better understanding of its corresponding chemical reactivity is obtained. Nevertheless, the rigorous foundation of CR-DFT has been overshadowed by the problems coming from the exact dependence of the energy density functional with respect to the number of electrons, which is known to be a set of straight-lines interconnecting contiguous N-integer pure states.12-14 This dependence has brought difficulties in the evaluation of higher order response indexes with respect to N, since it is well documented that the corresponding second order response function, called chemical hardness, behaves as a Dirac delta function in the electron number.1, 14 As an alternative, it is common to consider interpolation models between the integer values of N, with the prerequisite of being N-differentiable. There exist several flavors of these models in the literature,15-20 each of them carrying its own issues,21 but the most widely used is the parabolic interpolation. The parabolic interpolation of the E vs N profile was initially proposed by Parr and Pearson,16 and has the following expression,
E N = E N − 12 (IP + EA) N − N 0 + 12 (IP − EA) N − N 0
2
,
(4)
0
where IP = E N −1 − EN and EA = E N − E N +1 are the vertical first ionization potential and 0
0
0
0
the vertical first electron affinity relative to the E N reference state. (We use E N to denote 0 the energy of the N-electron ground state at the geometry of the reference state.) From Eq.(4) one obtains the following results at first and second order,
∂E = − χ = − 12 (IP + EA) , ∂N υ (r )
µ=
(5)
4 ACS Paragon Plus Environment
Page 5 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
and
∂µ = η PP = IP − EA , ∂N υ (r )
η=
where µ = − χ
is the negative of Mulliken’s electronegativity9,
definition of chemical hardness proposed by Parr and Pearson.16,
(6)
22
and η = η PP is the
23-25
A few years later,
based on a report by Maynard et al.,26 Parr proposed a new reactivity indicator called the electrophilicity index,27-30 which corresponds to the change in the energy associated with the maximum amount of electronic charge that a chemical species can accept; i.e., the electrophilicity is the difference in energy between the system in its reference state and the system when it is saturated with electrons. The electrophilicity index therefore corresponds to the minimum energy value in the E vs N parabolic interpolation function,
∆E
parab min
µ2 =− , 2η
(7)
where µ and η are given in Eqs. (5) and (6), respectively, and we have used the superscript “parab” to indicate that this quantity has been defined within the E vs N parabolic model. The maximum charge accepted by a particular chemical species dictates the position of the minimum in the E vs N quadratic function and it can be obtained from, parab ∆N max =−
µ . η
(8)
In both of these equations, it is assumed that the chemical potential is negative. If it is positive, the reference system is supersaturated with electrons and, consequently, it is reasonable to define the electrophilicity using Eq. (7), but multipled by minus one.
5 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 6 of 37
Numerical results have proven that the electrophilicity index defined in Eq.(7) is closely related to experimental properties for a plethora of chemical phenomena.31-35 Furthermore, theoretical investigations have brought some insight about possible fundamental principles related to this quantity. The minimum electrophilicity principle asserts that a chemical transformation must evolve towards the state with the lowest electrophilicity value,36-38 whereas the electrophilicity equalization principle states that at equilibrium, all the constituents of a chemical species must display the same electrophilicity value than the species as a whole.39-40 Recently, an equivalency between the minimum electrophilicity principle and the Hard/Soft Acid/Base (HSAB) principle41 has been suggested, after analyzing the equilibrium conditions of the same double-exchange reaction scheme previously used by Ayers in his fundamental derivation of the HSAB principle.42-44 The electrophilicity index has been conceived as a measure of the electroaccepting “potency” of the chemical species and new concepts like the electrodonating and the electroaccepting powers have been developed by distinguishing between the potency of a reagent in donating (accepting) charge-transfer processes.45-46 Recently, Miranda-Quintana has highlighted several issues regarding the definition of the electrophilicity index given in Eq. (7).47 He discussed, for example, that for systems with very low electron affinities the “parabolic” electrophilicity index would be mainly determined by the ionization potential of the species, an inconsistent result since this quantity is related to an electron release and not to an electron attachment process. Given the mathematical similarities between ∆Emin (Eq. (7)) and ∆N max parab
parab
(Eq. (8)), this last
statement is also true for the maximum charge accepted by an electronic species. He also pointed out that Eq. (7) is only valid for systems close to electron saturation and cannot by 6 ACS Paragon Plus Environment
Page 7 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
applied for highly-charged cations. This analysis encouraged him to define the electrophilicity index within the framework of the Helmholtz free energy potential, here denoted as ” A ”, mainly focusing on the zero temperature limit case. He obtained a general zero temperature limit working equation in which the electrophilicity is expressed as the sum of all the possible electron affinities consistent with the number of electrons in the chemical species under consideration. This definition, called the thermodynamic electrophilicity, avoids issues that appear in the parabolic electrophilicity index, since it
does not depend on the ionization potential and can be applied for cations with any charge value. Nevertheless, it is not clear how this new definition can be used to obtain new information regarding the chemical reactivity of the species or how this exact result can be connected to the widely (and successfully) used parabolic electrophilicity index. In this work, we present a systematic and detailed analysis of the behavior of A as a function of the number of electrons at any temperature. It will be shown that A(N) resembles a quadratic function in the number of electrons even at relatively low temperature values. Through the minimization of the Helmholtz free energy with respect to the number of electrons we obtain an electrophilicity index which strongly depends on (but is not equal to) the electron affinity of the electronic species. Through this approach, we were able to find the temperature conditions under which the parabolic and our electrophilicity index are equivalent. As a final remark, we also show that within our formalism, a nucleophilicity index defined analogously to the electrophilicity index is not justifiable even at very high temperatures.
7 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 8 of 37
2.0 Theoretical development. 2.1 Helmholtz free energy in the Grand Canonical formulation.
It is more convenient to study the dependence of the Helmholtz free energy on the number of electrons within the framework of the Grand Canonical Potential (GCP), Ω , because then the language of fractional electron number is naturally introduced.48-49 In the GCP formalism molecules are treated as open reactive systems, able to exchange energy and electrons with external reservoirs with temperature T and chemical potential µ Bath . In several recent reports we have exploited the mathematical features of the GCP in the study of chemical reactivity of electronic systems.50-57 Within this framework we have found that the N-differentiability of the average energy density functional with respect to the average number of electrons is guaranteed, ensuring the existence of arbitrarily high-order response coefficients for both the average electronic energy51 and the average electron density50. Likewise, we have demonstrated that a plausible definition of “problematic” quantities, like the local hardness and the local chemical potential, can be formulated within this framework.54,
58
Using the GCP formalism it has also been possible to develop a new
reactivity index associated with the energetic interactions between electronic species at the very earlier stages of a chemical reaction, namely the global heat capacity of the electronic systems,53 providing a deeper understanding of how a chemical reaction proceeds.
() ()
At a given temperature T and modified potential u r = υ r − µ Bath , the grand
()
potential is a unique functional of the equilibrium electron density ρ r ,1, 59-60
Ω[ ρ (r)] = E[ ρ (r)] − T ST [ ρ (r)] − µ Bath N[ ρ (r)],
(9) 8
ACS Paragon Plus Environment
Page 9 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
where E[ ρ (r )] , ST [ρ (r)] and N[ ρ (r )] are the average electronic energy, the entropy and the average-electron number density functionals, respectively. (See Ref. [50] for detailed definitions of these quantities.) One of the most important features of the GCP electron density is that it has been designed to integrate to fractional values of the electron number,
∫ ρ (r) dr = N
0
+ ω [ ρ (r)] ,
(10)
where N 0 is an integer, N0 ≡ N [ ρ (r)] , and −1 ≤ ω[ ρ (r)] < 1 is the fractional part of the electrons number. Similarly, the GCP can be expressed in terms of the Grand Canonical partition function, Ξ ,
Ω[ ρ (r)] = − β −1T lnΞ = − β −1 ln ∑ exp β −1 (E N ,i − µ Bath N ) ,
(11)
N ,i
where we have introduced the usual notation for the thermodynamic beta, β = 1/ k BT ( k B is Boltzmann’s constant) and the electronic energy of state i of the N-electron system, E N ,i . Combining Eqs. (9) and (11) one gets,
E[ ρ (r)] − T ST [ ρ (r)] = µ Bath N[ρ (r)]− β −1 lnΞ
,
(12)
where one can identify,
A[ ρ (r)] = E[ ρ (r)]− T ST [ ρ (r)] = µ Bath N[ ρ (r)] − β −1 lnΞ ,
(13)
as the Helmholtz Free Energy.
9 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 10 of 37
2.1.1 Helmholtz free energy as a function of the chemical potential of the reservoir.
To obtain a practical expression for A as a function of µ Bath , it is necessary to define the size of our ensemble. In a recent study,55 we have shown that the three-state ensemble model, constituted by the ground states of the N 0 , N 0 − 1 and N 0 + 1 electron systems, provides a satisfactory description of the reactivity of the species at any temperature value of chemical interest. For this ensemble model, Eq. (11) can be expressed as
{
}
Ω[ ρ (r)] = E N − µ Bath N 0 − β −1 ln 1+ exp β ( EA + µ Bath ) + exp − β (IP + µ Bath ) , 0
(14)
Inserting Eq. (14) in Eq. (13) one obtains,
∆A[ ρ (r)] = ∆E[ ρ (r)] − T ∆ST [ ρ (r)]
{
= µ Bath ∆N[ ρ (r)] − β −1 ln 1+ exp β (EA + µ Bath ) + exp − β (IP + µ Bath )
}
,
(15)
where ∆X[ ρ (r)] = X[ρ (r)]− X N measures the deviation of the average property X[ ρ (r )] 0
from the value of the corresponding reference state ψ N0 . It is known that for the three state ensemble model under consideration,51
∆N[ ρ (r)] = ω =
exp β (EA + µ Bath ) − exp − β ( IP + µ Bath ) 1+ exp β (EA + µ Bath ) + exp − β (IP + µ Bath )
.
.
(16)
Inserting Eq. (16) in Eq. (15), one obtains
10 ACS Paragon Plus Environment
Page 11 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
exp β (EA + µ Bath ) − exp − β ( IP + µ Bath ) ∆A[ ρ (r)] = µ Bath 1+ exp β ( EA + µ Bath ) + exp − β (IP + µ Bath ) , − β −1 ln 1+ exp β ( EA + µ Bath ) + exp − β (IP + µ Bath )
{
(17)
}
which explicitly depends on the chemical potential of the reservoir, µ Bath , and on experimental properties of the reference electronic state. In the zero temperature limit, the second term at the right-hand side of Eq. (17) vanishes, and the first term shows the same behavior of ω in this limit, that is,50, 61 −1 1 − 2 lim ω = 0 T →0 1 2 1
for µ Bath < − I for µ Bath = − I for − I < µ Bath < − A .
(18)
for µ Bath = − A for µ Bath > − A
for the three-state ensemble model under consideration. Therefore, at the zero temperature limit, the Helhomltz free energy behaves as a Heaviside function in the chemical potential of the electrons reservoir. Nevertheless, when the temperature increases, the contributions from the second term on the right-hand side of Eq. (17) smoothen the A[ ρ (r )] vs µ Bath profile. In Figure 1a we display the A[ ρ (r )] vs µ Bath profile of the Carbon atom for a range of temperature values up to 30000 K. It can be observed that, as the system reaches a given temperature value (≈ 5000 K), the A[ ρ (r )] vs µ Bath profile smoothens and the curve displays a minimum. This minimum is located in the positive fractional charge interval
(µ
Bath
> − ( IP + EA) / 2
)
and, as it will be discussed in the next subsection, it is a
consequence of the electronic entropy. For higher temperatures, where the entropy
11 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 12 of 37
contribution is not negligible, the A[ ρ (r )] vs µ Bath profile exhibits clearly a parabolic shape (although A does not display precisely a quadratic dependence on µ Bath ) a)
b)
c)
12 ACS Paragon Plus Environment
Page 13 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Figure 1. a) A[ ρ (r)] vs µ Bath , b) A[ ρ (r)] vs ω and c) E[ ρ (r)] vs ω profiles for the
(
)
0 Carbon atom µe = − 12 ( IP + AE ) = 6.26 eV , at the indicated temperature ranges.
2.1.2 Helmholtz free energy as a function of the number of electrons. The main benefit of defining the Helmholtz free energy as a functional of the GCP equilibrium electron density is evident from Eq. (10). Now, for a fixed value of temperature
T and modified potential u(r ) , any equilibrium value of A can be characterized by a given value of the fractional electrons number ω . Nevertheless, in order to analytically write A as a function of the fractional charge, one must isolate µ Bath from Eq. (16), a task that is practical only for the three-state ensemble model mentioned above, since the order of the polynomial to be solved increase with the number of pure states considered to build the ensemble. The dependence of µ Bath against the fractional charge ω for the ensemble model under consideration, splits in the following two roots,51
(α − ω ) , 2(1+ ω )
(19)
(α + ω ) , 2(1− ω )
(20)
− µ Bath = − IP − β −1 ln
+ µ Bath = − EA + β −1 ln
13 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 14 of 37
where,
(
)
α = ω 2 + 4(1− ω 2 )exp[β (IP − EA)]
1
2
.
(21)
Note that Eqs. (19) and (20) provide results only for non-positive and non-negative values of ω , respectively, and consequently, both equations are needed in order to describe the chemical potential of the bath over the full range of fractional charge and temperature values. Thus, according to the results in Eqs. (19) and (20), one may have two different expressions for A as a function of the fractional charge, one for negative values, using Eq. (19) and another for positive values, using Eq. (20). Eqs. (19) and (20) can be equivalently expressed as,
− exp[− β (I + µ Bath )] =
− exp[β (A+ µ Bath )] =
α −ω , 2(1+ ω )
(22)
α +ω , 2(1− ω )
(23)
respectively. Eqs. (22) and (23) can be now inserted in Eq. (11) to obtain the following expression for the Grand Potential,
1+ α Ω = EN0 − µ Bath N 0 − β −1 ln . 2 1 − ω
(24)
Inserting Eq. (24) in Eq. (13) the Helmholtz free energy is thus expressed as follows,
1+α ∆A = µBath ω − β −1 ln , 2 1 − ω
(25)
14 ACS Paragon Plus Environment
Page 15 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
where ω is defined in Eq. (16). The following expressions for A for negative fractional charge,
(α − ω ) 1+ α + ln ∆A− = −ω IP − β −1 ω ln 1− ω 2 2(1+ ω )
,
(26)
and for positive fractional charge,
(α + ω ) 1− ω 2 + ln ∆A+ = −ω EA + β −1 ω ln 2(1− ω ) 1+ α
(27)
are obtained simply by inserting Eq. (19) or (20), respectively, into Eq. (25). For temperature values for which β −1 is small compared to the band gap ( IP − EA) , the expression,
µ Bath ≅ µe + β −1
(
)
ω α+ω , ln α 2(1− ω )
(28)
provides the same numerical results as Eqs. (19) and (20), in their corresponding interval values of fractional charge. In Eq. (28) ω represents the absolute value of the fractional charge whereas µe , the electronic chemical potential, is the response function of the average electronic energy, E , with respect to changes in the average number of electrons at constant temperature and external potential, which for the three state ensemble model, it is known to have the following mathematical expression,51 1 ω ∂E = − ( IP + EA) + ( IP − EA) . 2 α ∂N T ,υ (r )
µe =
(29) 15
ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 16 of 37
Thus, through Eq. (28) one can obtain a unique expression for A as a function of ω , valid for the mentioned temperature range,
(
)
ω 2 α + ω 1− ω 2 + ln ∆A ≅ −ω µe + β −1 ln . α 1+ α 2(1− ω )
(30)
Figure 1b shows the A vs ω profile for the Carbon atom, for a range of temperature values up to 30000 K. At low temperatures, the dependence of A[ ρ (r )] as a function of ω resembles a straight-line profile. Nevertheless, as the temperature increases, this profile smoothens and A (ω ) displays a minimum for temperatures beyond ≈ 5000 K. As in the case of the A vs. µ Bath , this minimum is located in the interval of positive values of the fractional charge. Which contribution to A is responsible for the minimum value observed in the A (ω ) ? As can be seen in Eq. (15), the Helmholtz free energy is constituted by the sum of
two terms, the average energy, E , and the entropy times the temperature, T S . Within the three states ensemble model under consideration, the average electronic energy is given by,51
∆E =
IP exp [ − β ( IP + µ Bath )] − EA exp [ β ( EA + µ Bath ) ]
1 + exp [ β ( EA + µ Bath )] + exp [ − β ( IP + µ Bath ) ]
.
(31)
To obtain an explicit and unique expression for the dependence of the average electronic energy as a function of the fractional charge, one inserts Eqs. (22) and (23) into Eq. (31), and after some algebraic manipulations one finds that,
16 ACS Paragon Plus Environment
Page 17 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
E = EN0 − 12 ( IP + EA)ω + 12 ( IP − EA)
α +ω2 . 1+ α
(32)
Notice that this would be exactly Parr and Pearson’s quadratic model for the energy except that α is a function of the fractional charge, ω . Eq. (32) constitutes the exact dependence of the (average) energy density functional with respect to the (average) number of electrons, valid for the three states ensemble model under consideration at any temperature value (even in the zero temperature limit). One may verify for example, that the exact result obtained for the finite temperature electronic chemical potential given in Eq. (29) can be reproduced if one takes the partial derivative of Eq. (32) with respect to the fractional charge ω , at constant temperature and external potential (this result was previously obtained by defining E + and E − through Eqs. (19) and (20) respectively; then, performing the corresponding partial differentiation, one recovers Eq. (29)51. In Figure 1c we report the E vs ω profile for the Carbon atom at temperature values up to 60000 K, using Eq. (32). Unlike the A vs
µ Bath
and the A vs ω profiles, the
E vs ω profile does not display a minimum, but the expected result ∆ E min = − EA , which is observed even at very high temperature values, implying that the minimum observed for the A[ ρ (r )] vs ω profile, cannot be attributed to the energy contribution.
To analyze the entropic contribution to the Helmholtz free energy, we conveniently consider the analytic and simpler dependence of A with respect to positive values of the fractional charge, Eq. (27). To simplify our treatment, we also consider the low temperature expression of Eq. (21), where β −1
( I − A) and therefore α ≈ ω . Under these conditions,
Eq. (27) can be simplified to 17 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
∆A = −ω EA + β −1 {ω ln ω + (1− ω ) ln(1− ω )} .
Page 18 of 37
(33)
The first term at the right-hand side of Eq. (33) can be identified as ∆E , i.e.,
∆E = −ω EA .51 The second term in Eq. (33) represents the entropic contribution to A . One can verify that this term displays the same mathematical expression than the entropy of a Fermi-Dirac distribution, 48-49 and therefore, ∆ ST + = k B {ω ln ω + (1 − ω ) ln(1 − ω )} .
Equation (34) has a minimum at ω =
1
2
(34)
and an extremum value for the entropic
contribution is thus expected. The presence of this extremum value explains the curvature of the A vs ω profile for temperature values where the entropy contribution is nonnegligible. When the temperature value is close to the zero temperature limit, the minimum in the A vs ω is located at ω = 1 or correspondingly, ∆ Amin = − EA , which is the same result obtained by Miranda-Quintana.47 As the temperature increases, the entropic term becomes more important and the minimum is gradually displaced towards the value ω =
1
2
. Recall that the entropic term, as a function of ω , has a linear dependence on the temperature only for relatively low temperature values (see Eq. (33)), where the approximation α ≈ ω is still valid. For situations where this approximation cannot be used, the temperature dependent term present in the expression of α (see Eq. (21)), must be also considered. Consequently, the position of the minimum in the A vs ω profile will depend on the temperature dependence of α at higher temperatures. To end this section we want to show that under the low-temperature conditions, the
A vs ω profile can be reasonably approximated by a polynomial in ω . To do so, we 18 ACS Paragon Plus Environment
Page 19 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
expand in series the two logarithm functions in Eq. (33). The expression in Eq. (33) has divergent first derivatives at the edges of the interval 0 ≤ ω ≤ 1 caused by the terms ω ln ω (when ω → 0 ) and (1− ω ) ln (1− ω ) (when ω → 1). Thus, a series expansion is not reliable. We will instead use the approximation (1− ω ) ln(1− ω ) ≈ ω 2 ln ω 2
∀ 0< ω < 1 .
(35)
which gives a new (approximate) expression for the entropy term present in Eq. (34), ∆ ST + ≈ k B {ω ln ω + ω 2 ln ω 2 } = k B
{(ω + 2ω ) ln ω} . 2
(36)
This expression still has a divergent first derivative as ω → 0 , but one can expand this function about ω = 1 , giving
.
(37)
As a second alternative, we have made a parabolic fit to the Helmholtz free energy by making use of the values obtained at the points ω = 0 , ω = 1 2 and ω = 1 , with Eq. (33). The resulting expression is,
(
) (
)
∆A+ ≈ − EA + β −1 ln16 ω + β −1 ln16 ω 2 ≈ −ω EA − β ω (1− ω )ln16 −1
.
(38)
Thus, the quadratic interpolation predicts “a parabolic entropy” which is equal to the second term in the right hand side of the second equality in Eq. (38). In Figure 3 we
19 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
compare the A vs ω profiles obtained with Eq. (37) using the second, third, fourth and fifth order expansion approximations, and the parabolic one obtained from Eq. (38), with the exact result obtained through Eq. (27), for the Carbon atom at T = 10000 K. Since the
(
)
quantity ω −1 is always fractional and because the expansion coefficients in Eq. (37) decrease, the contribution from high order terms quickly vanish, indicating that the A vs
ω profile can be reliably approximated by a low-order polynomial at relatively low temperature values (see Figure 2). This last statement is supported by the parabolic interpolation given in Eq. (38), where it can be observed that the best quadratic function has a similar performance to that of the highest order polynomial (n=5, Figure 2). Likewise, it can be observed that due to the divergent behavior of first derivatives of Eq. (33) mentioned above, Eqs. (37) and (38) are less accurate for integer values of the electronic charge, indicating that any polynomial interpolation will exhibit an small deviation from the exact thermodynamic profile, at low fractional charge values.
∆A (eV)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 20 of 37
0.00
Exact. (Eq. (27)) Eq. (38) n=2 n=3 n=4 n=5
-0.25 -0.50 -0.75 -1.00 -1.25 -1.50 0
0.2
0.4
0.6
0.8
1
ω 20
ACS Paragon Plus Environment
Page 21 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
Figure 2. Comparison between the exact A[ ρ (r )] vs ω profile, the polynomial expansion in Eq. (37), evaluated at several orders of approximation, and the parabolic interpolation given in Eq. (38), for the Carbon atom at T = 10000K.
2.2 The electrophilicity index at finite temperatures. 2.2.1 The Helmholtz electrophilicity. In this section we are going to establish a connection between the electrophilicity index and the extremum of A as a function of ω . This is motivated by the fact that A (ω ) has a minimum for a fractional charge, 0 < ω < 1 , at sufficiently high temperature, while E (ω ) does not.
The minimum in the A (ω ) can be determined by finding the value of the fractional charge, ω max , for which dA dω = 0 . In accord with standard statistical mechanics and, explicitly, from Eq. (25) one can verify that in equilibrium,48-49
0=
∂A = µ Bath ∂ω
,
(39)
Eq. (39) directly imposes the zero value for the thermodynamic chemical potential of the reservoir, consistent with the electron saturation condition. Thus, ω max is the maximum charge that an electronic system is able to accept at a given temperature, and can be obtained from Eq. (20), by finding the ω value for which,
(
)
0 = α + ω 1+ 2exp β EA − 2exp β EA .
(40)
21 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 22 of 37
Due to the mathematical expression of α (Eq. (21)), Eq. (40) must be solved numerically. Nevertheless, for temperature values where α ≈ ω , the solution of Eq. (40) provides the following simple result,
ω max = 12 + 12 tanh[ β2 EA].
(41)
It is worth noting that according to Eq. (41), the maximum charge accepted by an electronic species only depends on its corresponding EA value, avoiding in this way the dependence on the IP value exhibited by ∆N max
parabol
for chemical species with EA values
close to zero. Once ω max is determined either by solving Eq. (40) numerically or directly from the (accurate) approximate result in Eq. (41), one can thus obtain,
∆Amin = ∆A ω
,
(42)
max
that is, the Helmholtz free energy evaluated at ω max . We call this quantity the Helmholtz Electrophilicity Index (HEI). In Figure 3a we compare the ω max values obtained by solving Eq. (40) ( ω max ) with those directly determined from Eq. (41) ( ω max ) for the Carbon and exact
approx
Fluorine atoms, as well as for the diatomic molecule NO, at temperature values up to 40000 K, whereas in Figure 3b we compare ∆Amin evaluated at evaluated at
exact exact ω max , ∆Amin , with
∆Amin
approx approx ω max , ∆Amin , in the same temperature range for the three chemical species
mentioned above. Figure 3a shows that the approximation given in Eq. (41) provides exactly the same results than the ones obtained by the exact numerical solution of Eq. (40) for temperature values up to ≈ 20000 K for the Carbon and Fluorine atoms, and for the NO
22 ACS Paragon Plus Environment
Page 23 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
molecule both equations provide the same results for the whole temperature range under consideration. Figure 3b shows that the Helmholtz electrophilicity is even less sensitive to the approximation mentioned above, since the ∆Amin
approx
vs T profile exactly reproduce the
exact ∆Amin vs T profile for the three chemical species on the entire temperature interval studied.
23 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
ωmax
a)
1.0
F
0.9 0.8
C 0.7 0.6
NO 0.5 0.4 0
10,000
20,000
30,000
40,000 T (K)
b)
∆Amin (eV)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 24 of 37
0.0
NO -1.0
C -2.0
-3.0
F -4.0
-5.0 0
10,000
20,000
30,000
40,000 T (K)
Figure3. Comparison between a)
exact approx exact approx ω max and ω max , and b) ∆Amin and ∆Amin , values for
the Carbon and the Fluorine atoms and for the NO molecule, at temperatures up to 40000 K (continuous lines represent exact profiles whereas dotted lines correspond to the approximations). 24 ACS Paragon Plus Environment
Page 25 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
2.2.2 The temperature of the E vs N parabolic interpolation. As mentioned in our introduction, both ∆N max and ∆Emin indexes have been parab
parab
successfully used to describe many reactivity trends even though they are based on the ad hoc assumption of a parabolic energy model. The empirical utility of these indicators suggests that they may be implicitly, and certainly approximately, modelling some real chemical information about the system. On the other hand, one might consider that the temperature values where the A (ω ) has a minimum at fractional charge are far higher than the temperatures at which most chemical phenomena are observed. However, it is now well-accepted that the “temperature” that one uses when exploring chemical reactivity concepts is not the true temperature, but an effective temperature that models, in a simple but nonetheless useful way, the fact that the interactions between reacting systems force them away from the equilibrium state they occupy if they were isolated. There are several attempts to quantify this effect in the literature. The first was to define a local temperature due to the kinetic energy of electrons, which is very high.62-64 One can similarly define a local electronic temperature as the conjugate variable to a suitably defined local or information theoretic entropy.65-67 In these interpretations, approaching reagents’ electrons can impart significant kinetic energy to (or induce significant disorder in) the electrons of a substrate molecule, and so a high value for the “effective temperature” is justified. From an alternative perspective, the effect of an approaching reagent is to perturbatively mix in the wavefunctions of higher-energy states, just as temperature mixes in higher-energy wavefunction components.68-69 That is, the density matrix of a subsystem will not be a pure state even in the zero-temperature limit, and the occupation numbers for the subsystem density matrix can be (approximately) fit to a Fermi distribution function corresponding to 25 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 26 of 37
a substantial electronic temperature.68 This is similar to our observation about the electronic heat capacity, which indicates that at the onset of a chemical reaction, there is an electronic flux between the reagents that substantially elevates the temperature of the reactants, regardless the temperature value of the surroundings.53In the absence of a simple, practical, and convincing strategy for computing the effective electronic temperature and of the way it varies during a particular chemical reaction, a pragmatic approach seems advisable. Therefore we assume that ∆N max resembles the true ω max at high temperature values. This parab
gives one a rough approximation for the effective temperature, simply by solving parab ωmax (T ) = ∆Nmax . We call this temperature the parabolic temperature T parab and the
maximum fractional charge exchanged at this temperature is denoted as Helmholtz free energy evaluated at report T parab , ∆Emin , parab
parab ω max . Then, the
parab parab ω max and T parab is denoted as ∆Amin . In Table 1 we
parab parab ω max and ∆Amin values for six neutral atoms and for four simple
molecules. (When solving for T parab we used Eq. (40) and not the low-temperature approximation in Eq. (41).) To verify that the curvature of ∆Amin
parab
curvature of ∆Emin
parab
resembles to the
we are also comparing the band gap IP − EA with the inverse of the
thermodynamic softness, S , evaluated at T parab and
∂µ ∂2 A = Bath 2 ∂ω T parab ,ω parab ∂ω max
parab ω max , which is known to be,54, 61
= parab T parab ,ωmax
1 S
parab
α (1 − α )(1 − ω 2 ) =β 1−α 2 T
parab T parab ,ωmax
,
(43)
−1
parab
parab ,ωmax
26 ACS Paragon Plus Environment
Page 27 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
for the three-state ensemble model under consideration. Table 1 shows that the values of
T parab range from 9.3x103 K to 36.9x103 K, with an average value of 21.0 x103 K. It is worth to mention that in our previous report55 we have also shown that even at these temperature conditions the three states ensemble model is a reliable approximation to describe the chemical reactivity of electronic systems. It can be observed that T parab increases as the band gap IP − EA increases, indicating that a hard system requires larger amount of energy to react. It is worth mentioning that the same dependence of temperature on band gap was observed from our study of the global heat capacity of electronic systems.53 It is interesting to note that, in almost all cases, the ∆Amin
parab
to the ∆Emin
parab
values are very close
values and both of them follow exactly the same trend, R 2 = 0.981 . One can
also see that the band gap IP − EA and 1/ S parab strongly resemble each other. These results thus confirm that at T parab the widely used E vs N interpolation model is usually an accurate approximation of the exact A vs ω profile. The approximation is less accurate for the NO molecule, where the 1/ S parab value significantly differs from the corresponding band gap value, IP − EA . This discrepancy is a direct consequence of the very small EA value for this molecule, which causes the parabolic model to give a questionable description of both, the electron attachment and the electrophilicity index, as was pointed out by MirandaQuintana.47
27 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 28 of 37
Table 1. Comparison between some reactivity indexes obtained from the parabolic energy profile, with their Helmholtz free energy counterpart evaluated at T = T parab . Temperature units are given in Kelvin degrees. Energy values are given in eV. Experimental data were taken from ref [70]. Species
IP
EA
C Cs F K Ni Si SO2 S2 CH3I NO
11.26 3.89 17.42 4.34 7.64 8.15 12.3 9.36 9.5 9.26
1.27 0.47 3.4 0.5 1.15 1.39 1.1 1.66 0.2 0.026
parab parab parab − ( IP + EA) / 2 IP − EA ∆Emin T parab x103 ω max ∆Amin 1/ S parab
-6.26 -2.18 -10.41 -2.42 -4.39 -4.77 -6.70 -5.51 -4.85 -4.64
9.99 3.42 14.02 3.84 6.49 6.76 11.2 7.70 9.30 9.234
-1.96 -0.70 -3.87 -0.76 -1.49 -1.68 -2.00 -1.97 -1.27 -1.17
27.2 9.33 36.9 10.5 17.6 18.2 30.3 20.6 22.0 17.4
0.627 0.637 0.743 0.630 0.677 0.706 0.598 0.716 0.522 0.503
-2.52 -0.88 -4.49 -0.97 -1.83 -2.02 -2.61 -2.35 -1.52 -1.08
9.78 3.40 16.45 3.78 6.82 7.43 10.54 8.60 7.42 5.94
2.3 The inexistence of the nucleophilicity index. Unfortunately, it is impossible to define an analogous nucleophilicity index by extending the preceding arguments to negative fractional charge, −1 < ω < 0 . A (ω ) never has a minimum for ω < 0 because the equilibrium condition given in Eq. (39) can only be fulfilled for positive values of the fractional charge, with ω → 0 only for high temperatures. The Helmholtz free energy can therefore only be minimized for the acceptation process of electronic charge and not for the donation process. While there is no minimum in the free energy for electron-donation, there is a minimum in the electronic entropy. To show this, we use the α ≈ ω low temperature approximation in Eq. (26), obtaining in this way the following expression for the Helmholtz free energy change for negative values of fractional charge,
28 ACS Paragon Plus Environment
Page 29 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
α −ω + (1+ ω )ln(1+ ω ) , ∆A = −ω IP − β −1 ω ln 2
(44)
This last equation can be equivalently expressed as follow,
∆A = ω IP + β −1 { ω ln ω + (1 − ω ) ln(1 − ω )} .
(45)
where ω represents the absolute value of the fractional charge. Eq. (45) displays an entropic contribution numerically equivalent to Eq. (36) and thus, one might think that Eq. (26) would display a minimum for a certain temperature value. However, since the ionization potential is bigger than the electron affinity, the temperature required to make the entropy contribution significant, is expected to be much higher for the donation process than in the acceptation process. Nevertheless, when the temperature is high enough, the minimum exhibited by the exact entropy contribution (second term at the right-hand of Eq. (26)) is displaced towards ω = 0 . The temperature required for the entropy contribution to overcome the energetic one ( ∆E[ ρ (r )] = −ω IP ) is proportional to IP (with ω as the proportionality constant), whereas the minimum entropy value is displaced towards ω