Stationary Conditions of the Electron Density Along ... - ACS Publications

Subscriber access provided by Olson Library | Northern Michigan University

Article

On the Stationary Conditions of the Electron Density Along the Reaction Path: Connection with Conceptual DFT and Information Theory Carlos A. Gonzalez, Emilio Squitieri, Hector Julio Franco, and Luis C. Rincon J. Phys. Chem. A, Just Accepted Manuscript • Publication Date (Web): 28 Dec 2016 Downloaded from http://pubs.acs.org on December 28, 2016

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.

The Journal of Physical Chemistry A 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 55

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

The Journal of Physical Chemistry

On the Stationary Conditions of the Electron Density Along

the

Reaction

Path:

Connection

with

Conceptual DFT and Information Theory Carlos A. Gonzalez1,*, Emilio Squitieri2, Hector J. Franco2, Luis C. Rincon3,4 1. Materials Measurement Laboratory, Chemical Sciences Division, National Institute of Standards and Technology, Gaithersburg, Maryland 20899, USA. 2. Escuela de Química, Facultad de Ciencias, Universidad Central de Venezuela, Apartado 47102, Caracas 1041-A, República Bolivariana de Venezuela. 3. Departamento de Química, Universidad de Los Andes, Mérida, República Bolivariana de Venezuela. 4. Grupo de Química Computacional y Teórica, Universidad San Francisco de Quito, Colegio Politecnico de Ciencias e Ingeniería, Diego de Robles y Vía Interoceánica, 17-1200-841 Quito, Ecuador. * Corresponding author

ACS Paragon Plus Environment

1


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 55

ABSTRACT

The Kohn-Sham Density Functional Theory (DFT) formalism has been used to investigate the influence of the stationary behavior of the electron density ((; )) along a minimum energy path on the corresponding stationary conditions observed in the total potential energy of the reactive system, information theory measures (Shannon information entropy and Onicescu information energy), and chemical reactivity indexes (the chemical hardness). The theoretical treatment presented in this work, combined with DFT calculations on 3 different test reactions: H′ + H2, H ′ + CH4 and H- + CH4, suggest that for any reactive system, properties that can be cast as a functional of the electron density, must exhibit stationary points along the IRC path modulated by the corresponding stationary behavior of the electron density.

Introduction The use of modern computational chemistry methodologies in the prediction of molecular properties has become increasingly popular, mainly due to significant improvements in the algorithms, the accuracy of the methods, and the advent of powerful computer resources. This is particularly true in the area of thermochemistry where researchers in industry and academia perform quantum chemistry calculations on a routine basis. Due to this impressive success, there has been an increased interest in the application of state-of-the-art quantum chemistry methodologies in the study of physical and chemical phenomena in a wide variety of technological areas including biomedicine, chemical catalysis, environmental chemistry, combustion and alternative sources of energy, climate assessment, nanoelectronics, and the


2

Page 3 of 55

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


rational design of materials at the nano-scale. The field has reached such a state of maturity that modern quantum chemistry methods are currently used not only for the interpretation of experimental measurements but also to guide the experimental design for measurements of novel physical and chemical properties of a wide variety of chemical system and advanced materials. One of the areas where major advances in quantum chemistry methodologies have made a significant impact is in the quantitative description of mechanisms, energetics, and dynamics governing chemical reactions. Today, scientists routinely employ ab initio quantum chemical methodologies to map the potential energy surfaces (PES) of reactive systems in the vicinity of the minimum energy reaction path (MEP) to study the transformation of reactants into products. One of the most efficient ways to accomplish this entails the use of the intrinsic reaction coordinate (IRC) concept originally developed by Fukui1-5. Due to the availability of analytical gradient methods6-9, the implementation of reaction path following algorithms in the realm of ab initio quantum mechanical calculations has allowed scientists to study the energetics, mechanisms and electronic structure properties of reacting systems involving relatively large polyatomic systems10-16. In conjunction with these advances, the quantum chemistry community has also witnessed a major revolution due to the development and implementation of relatively computationally efficient quantum chemical methodologies based on Density Functional Theory (DFT). The seminal work of Hohenberg and Kohn17 in 1964 (Hohenberg-Kohn Theorems) followed by Kohn and Sham18 in 1965 (Kohn-Sham Equations), made it possible for computational scientists to perform calculations of physical and chemical properties along a minimum energy reaction path of relatively larger molecular systems at a computational cost significantly lower compared to traditional ab initio wavefunction-based electronic structure theory.


3


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 55

One of the most active areas of research in DFT has been the so-called “Conceptual DFT”19, that provides a theoretical sound basis describing concepts such as electronegativity, chemical hardness, and chemical softness, and their relations to chemical principles such as Hard-Soft Acid-Base (HSAB), Electronegative Equalization, and Maximum Hardness. The origins of Conceptual DFT can be traced back to the pivotal work of Parr20,21 who developed the theoretical formalism that establishes the direct relation between the chemical potential and the fundamental equations of DFT. Due to the development of computationally efficient DFT algorithms implemented in a variety of popular commercial and non-commercial quantum chemistry software packages, the literature has seen an explosion of applications employing Conceptual DFT. This is particularly true in the case of Chemical Reactivity Theory22-28 (CRT) that allows the connection between DFT concepts and the study of properties such as chemical bonding, reactivity, and dynamics. The scientific literature abounds with applications of Conceptual DFT in the study of changes of chemical reactivity indices (as well as other electronic structure properties) along reaction paths29-35. More recently, scientists have also made use of Information Theory (IT) concepts and performed calculations of the Fisher information, the Shannon information entropy, and the Kullback-Leibler information measures along the IRC on a series of chemical reactions36-44. The numerical results of these studies suggest that when a chemical reacting system moving along the IRC approaches the vicinity of the transition state (TS), properties such as the global hardness, the Shannon information entropy, dipole moment, and the Fisher Information reach a stationary point (maximum or minimum). Even though the existence of these stationary points can be interpreted in terms of properties such as bond breaking/formation, localization/de-localization of the electron density, as well as charge-transfer, to the authors’ knowledge, there does not seem


4

Page 5 of 55

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


to exist a solid theoretical argument to rationalize the origin and nature of the conditions leading to these stationary points. In this work, we investigate the conditions that lead to stationary points of the electronic density of a reacting system along the reaction path. We further show that the existence of these stationary points in the electron density can be used to explain the origin of the corresponding stationary points exhibited by quantities from IT such as the Shannon information entropy, and the Onicescu information energy (or Disequilibrium), as well as reactivity indexes such as the global hardness, as previously reported in the literature45.

Theoretical Basis Figure 1 depicts a typical IRC profile (total energy vs reaction coordinate s) corresponding to a MEP for a chemical reaction that describes the transformation of reactants (s -∞) into products (s +∞), connected through a transition state (TS) located at s = 0. As observed in Fig.1, the reactants and products are local minima on this curve, while the TS is a saddle point. According to the first Hohenberg-Kohn theorem17, the electron density, (), determines the ground state electronic wavefunction as well as any other electronic property (including the total energy) of any chemical system. Consequently, it follows that the electronic energy of a molecular system for a particular point s lying along the IRC should be a functional of the total electron density, which in turn is a function of the reaction coordinate s. This statement can be cast into the following relation:

→ (; ),

(1)

where the term (; ) has been used to explicitly indicate the dependence of the electron density on the reaction coordinate s, with indicating the coordinates of the electrons. ACS Paragon Plus Environment

5


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 55

For a particular point s along the IRC, the total potential energy (; ) is then given by the following expression:

(; ) = (; ) + ,

(2)

where is the classical nuclear-nuclear repulsion energy at point s: = ∑ ∑

.

(3)

In Equation (3), the indices α and β run over the total number of atoms ! , " are the nuclear charges, and # are the corresponding internuclear distances between atoms α and β. Taking the derivative of Equation (2) with respect to s leads to: $

%&' (();*) %+

, = - $

.&/ (();*) .(();*)

, $

For a stationary point so along the IRC, $ $

%&' (();*) %+

,

+5

= - $

.&/ (();*) .(();*)

%(();*)

, 0 1 + $

%+

%&' (();*) %+

, $

%(();*)

+5

%+

,.

(4)

, = 0, and Equation (4) becomes: +4

, 0 1 + $

%+

%233 (*)

+4

%233 (*) %+

,

+4

= 0.

(5)

The formalism presented in this work (Equations (4) and (5) above) provides the basis leading to a more fundamental understanding of this concept. In Eq. (5), the functional derivative $

.&/ (();*) .(();*)

,

+4

is the chemical potential46 7+4 evaluated at so. 7+4 = $

.&/ (();*) .(();*)

,

(6)

+4

The chemical potential is a measure of the escaping tendency of an electronic cloud from the molecular system. Taking into consideration that for any particular point s along the IRC 7+ is the same throughout all molecular space, and making use of Equation (5), it becomes clear that the following condition must be satisfied at any stationary point so along the IRC: 7+4 - $

%(();*) %+

, 0 1 + $ +4

%233 (*) %+

,

+4

= 0.

(7)


6

Page 7 of 55

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


Examination of Equation (7) indicates two possible cases for the existence of a stationary point along the reaction path: Case I: 7+4 - $

%(();*)

, 0 1 = $

%+

+4

%233 (*) %+

,

%233 (*)

+4

= 0,

(8)

Case II: -$

%(();*) %+

, 0 1 = − 9: < $ +4

;4

%+

, . +4

(9)

Equations (8) and (9) highlight how changes in the electron density of a reacting system lead to an extremal in the potential energy surface at any stationary point along the corresponding minimum energy reaction path. Thus, while in Case I a stationary point on the surface is reached when the first two terms in Equation (6) are identically equal to zero (implying an extremal condition in both, (; ) and VNN(s)), Case II leads to stationary points where these terms remain finite but cancel each other. Equations (4) – (9) relate three very important reactivity quantities at the transition state: the electronic chemical potential, 7+4 , the nuclear-nuclear reaction force, − $

%233 (*) %+

derivative of the electron density of the reactive systems along the reaction path, $

, , and the +4

%(();*) %+

, . To +4

our knowledge, the treatment discussed above recognizes for the first time that the interplay between these three reactivity indexes provides a set of rigorous conditions for the existence of transition states and in general, for the existence of any stationary point along the reaction path. Thus, the theoretical treatment discussed in this section offers a solid theoretical framework based on the behavior of the changes of the electron density along the reaction path (not previously examined in detail) that provides the basis for a more fundamental understanding of chemical reacting systems (and in particular the origin of a reaction’s potential energy barriers)


7


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 55

when compared to the conventional energy-based approach, in agreement with the density-based representation of reactivity proposed by Knoerr and Eberhart47. An alternative physical interpretation of the meaning of the derivative $

%(();*) %+

, can be found

by making use of the relation between the electron density and the functional derivative of the electronic energy with respect to the external potential (at constant N): (; s) = $

.&/ ();*)

,

.>

(10)

The derivative of Equation (10) with respect to s at the transition state follows the relation: $

%(();*) %+

.? (*)

/ , = − $ .>() , )

+4

(11)

where @ () is the electronic reaction force evaluated at the transition state: @ = − $

%&/ (*) %+

,

(12)

+4

The “reaction force” concept has been extensively used for the analysis of the IRC of many prototypical reactions48-51. According with Equations (10) – (12), $

%(();*) %+

, can then be interpreted as the local response

index of the electronic reaction force due to an external perturbation (reflected in changes in the geometry) in the vicinity of a particular point along the reaction path. Positive (negative) values in $

%(();*) %+

, indicates a decrease (increase) of the electronic reaction force due to an increment in

the external potential. This local index provides with a very solid framework to capture the nature of the stationary points along a reaction path from a DFT perspective. Finally, integration of $

%(();*) %+

, over all molecular space can be interpreted as the average response of the electronic

reaction force to changes in the molecular geometry of the reactive system along the reaction


8

Page 9 of 55

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


path. Thus, in the case that -

A %(();*) $ %+ , A +

4

= 0, the electronic reaction force becomes an

extreme at the transition state. Of particular interest in this work is the study of the behavior of the electron density in the vicinity of the transition state. According to Case I given by Equation (8) above, the integral -$

%(();*) %+

, 0 1 becomes zero at the TS. Evaluation of the derivative $ +5

%(();*) %+

,

+5

at B by

central differences can readily be used to demonstrate that for reaction paths involving symmetrical transitions structures (such as the ones found in a wide variety of thermo-neutral reactions), this derivative indeed becomes equal to zero at every point in the whole geometrical space spanned by the molecular system. In these cases, the corresponding transition structures exhibit an inversion plane perpendicular to the reaction path. Reflection through this plane leads to two identical geometries and consequently identical electron densities, which results in: $

%(();*) %+

,

+5

= 0.

(13)

with a similar stationary condition on the nuclear-nuclear repulsion energy (see Equation (3) above). Clearly, we still have to entertain the possibility of the existence of Case I for reactions exhibiting non-symmetrical transition structures, where - $ derivative itself $

%(();*) %+

%(();*) %+

, 0 1 becomes zero while the +5

, remains finite throughout some or all of the molecular space. So far, +5

our numerical calculations on 10 different reactive systems (including hydrogen abstractions, electrophilic additions, SN2 and Diels-Alder reactions, isomerizations and 1,2-hydrogen rearrangement reactions) have not been able to identify reactions illustrating these cases.


9


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 55

By contrast to Case I, reactions dominated by Case II stationary conditions, exhibit saddle points along the potential energy surface where the derivatives $

%(();*) %+

,

+5

(and their integral

along the whole molecular space) are finite and can only be cancelled out by the right side term in Eq. (9), suggesting the existence of transition structures lacking a symmetry plane of inversion similar to the one described above. It is quite possible to observe reacting systems where the integral - $

%(();*) %+

, 0 1 becomes zero in regions along the IRC that are not stationary points. +4

However, in these cases, the term − 9: < $ ;4

%233 (*) %+

, survives and the resulting derivative of the +4

potential energy with respect to the reaction coordinate evaluated at C is non-zero.

The Stationary Condition of the Electron Density and Shannon Entropy The past couple of decades have witnessed significant advances in the development of informational theoretical methodologies as also as applications in a wide range of scientific areas. In particular, the pioneer work by Sears, Parr and Dinur52, constitutes one of the first applications of information theory in quantum chemistry, and it has been followed by a plethora of investigations using both, wave- and density-based approaches53-56. Among the large number of IT methodologies currently available for the quantification of the amount of information associated with a probabilistic event, the ideas of information transmission and entropy introduced by Shannon57,58 in 1948 have become one of the most popular concepts adopted by scientists working in areas ranging from chemistry, biology and genetics, to linguistics and social sciences. The details of Shannon’s work and his remarkable contributions to IT and science in


10

Page 11 of 55

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


general are well beyond the scope of this work. Suffice it to mention that Shannon was able to connect probability distributions with the average information content of any particular system, leading to the following expression for what is called the “Shannon Entropy Functional”, S, or “Shannon information entropy” for short: D = − ∑F EF ln EF ,

(14)

where: E = IE , EK , ⋯ EM N is the probability distribution associated with an stochastic process with n possible outcomes (The natural logarithm in Eq. (14) implies that the information entropy is computed in units of nats). For an N-electronic molecular system reacting along a minimum energy path and characterized by a continuous probability distribution (normalized to one) described by the shape function, O(; ), the Shannon information entropy can be cast into the following form: DO(; ) = − -

A O(; ) ln O(; ) 0 1 , A

where s is the reaction coordinate, O(; ) =

(();+)

(15)

, (; ) is the electron density at s, and is

the 3N-dimensional vector containing the electrons’ Cartesian coordinates. The Shannon information entropy as presented in Eq. (15) provides a global metric for the delocalization or “spread” of the electron density of a particular molecular system. Alternatively, S can be interpreted as the average information loss due to the uncertainty in the electron position and consequently, S decreases as the electron density becomes more localized. Recent computational studies using DFT calculations on a series of simple chemical reactions indicate the existence of extremal values of the Shannon information entropy either in the vicinity of the TS or on particular points along the reaction path far from the TS30-33,35. This important finding suggests that a solid criterion for exploiting the properties of D(; ) could be used in the characterization of chemical reactivity along a particular MEP.


11


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 55

Taking the derivative of Equation (15) with respect to s, evaluating it at the TS (s = so), and imposing the stationary condition on DO(; C ) leads to: %P

$%+ ,

+5

= − -

A 1 A

+ ln O(; ) $

%R();+) %+

, +5 0 1 = 0.

(16)

As in the case with the potential energy previously discussed, Equation (16) leads to the following two possible stationary conditions: Case S.I: -

A %R();+) $ , A %+ +

5

0 1 = − -

A %R();+) ln O(; ) $ , A %+ +

5

0 1 = 0.

(17)

Case S.II: -

A %R();+) $ %+ , A +

5

0 1 = − -

A %R();+) ln O(; ) $ %+ , A +

5

0 1 .

(18)

Similarly to the result of the previous section, Case S.I is characteristic of reactions with transition structures exhibiting a symmetry inversion plane perpendicular to the reaction path, making the derivative $

%R();+) %+

,

+4

identically to zero at every point in the whole geometrical

space spanned by the molecular system. In reactions that entail transition structures lacking these symmetry inversion planes, the derivative $

%R();+) %+

,

+4

is finite and the extremal in Shannon

entropy can only be achieved if the two terms in Equation (18) cancel each other.

It is

interesting to notice how the nature of the stationary points in DO(; ) and (; ) are governed by similar behaviors of the electron density along the reaction path.

The Stationary Condition of the Electron Density and Onicescu’s Information Energy The Onicescu information energy or Disequilibrium, S) ()59, is a second-order entropic moment that provides a global measure of the information content in a quantum system. When


12

Page 13 of 55

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


applied to the electron distribution of a molecule, it measures the deviation of () from an equiprobability state or maximum Shannon entropy. Thus, while the Shannon entropy is used to measure “disorder” or “delocalization” of the electron distribution in a molecule, the disequilibrium provides a global measure of its “order” or “localization”. For an N-electronic molecular system reacting along a minimum energy path, S) (; ) can be expressed as (in volume units): S) (; ) = -

A O(; )K 0 1 , A

(19)

where as before, we have included the explicit dependence on the reaction coordinate s. Taking the derivative of Equation (19) with respect to s, and imposing the stationary condition on S) (; C ), leads to the following relation: $

%TU (();+) %+

,

+5

= 2 -

A %R();*) O(; B ) $ %+ , A +

5

0 1 = 2 〈$

%R();+) %+

, 〉 = 0 (20) +5

As with the potential energy and Shannon entropy, it is observed that in the case of a reacting systems involving transition structures with a symmetry inversion plane perpendicular to the reaction path, $

%TU (();+) %+

,

+5

= 0 given that the derivative

$

%(();+) %+

,

+4

=$

%R();*) %+

,

+5

is zero at

every point in the whole geometrical space spanned by the molecular system (Case D.I). For any other cases, the stationary condition on S) (; ) requires that the whole integral in Equation (20), (which provides a measure of the expectation value 〈$

%R();+) %+

, 〉) becomes zero (Case +5

D.II).

The Stationary Condition of the Electron Density and Chemical Hardness Hardness is one of the most popular concepts used in conceptual DFT. It was originally introduced by Pearson and collaborators in their study of the reactions between Lewis acids and


13


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 55

bases,24,60-62. Hardness is a measure of the stability of a molecule toward reactivity (i.e. larger values of hardness indicate a higher resistance of a molecular species to change its electron distribution and therefore, to react). Parr and Pearson defined the hardness, η, as the derivative of the chemical potential, 7, with respect to the number of electrons, N, at constant external potential, ν63:

Y = $ K

:

, = K $ >

Z&

, .

Z >

(21)

Assuming a quadratic relation between Y and N, and using a finite difference approximation, Equation (21) can be approximately computed by the following expression:

Y ≈

\ ] K

≈

&3^_ &3`_ K&3 K

,

(22)

where I and A are the ionization potential and the electron affinity, and the terms EN, EN-1 and EN+1 indicate the energies of the neutral, cationic and anionic species for the same geometrical structure. In a similar way, the chemical potential 7 is defined as:

7 ≈

\] K

≈

&3^_ &3`_ K

(23)

Given that the energy is a functional of the electron density, Equation (22) indicates that for a reactive system, the hardness is a functional of the electron density and that it varies with the reaction coordinate, Y(; ). Taking the derivative of Y(; ) with respect to s leads to: $

%a(();+) %+

, = -

A .a(();+) %(();+) $ .(();+) , $ %+ , 0 1 A

(24.a)

Making use of Equation (22) and working at constant external potential, the following expression for $

.a(();+) .(();+)

, is obtained: >


14

Page 15 of 55

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(();+) .(();+)

, ≈ K b7 + 7 – 2 7 d >

(24.b)

The extremal point in Y(; ) located at the TS, (s = s0), can then be found by imposing the following stationary condition: $

%a(();+) %+

,

+5

0 0 = K (7 + 7

0 )− 2 7

A %(();+) $ %+ , A +

5

0 1 = 0,

(24.c)

C C C where 7 , 7 , and, 7 , are the chemical potentials corresponding to the neutral, anionic and

cationic species evaluated at the TS. We find again that in reactions exhibiting highly symmetric transition structures the hardness becomes an extremal at the TS given that $

%(();+) %+

,

+5

= 0

C (Case H.I), while other extremals in Y(; ) can be found provided that either the (7 + C 7

C) − 2 7 factor or the integral -

A %(();+) $ %+ , A +

5

0 1 on the right side of the Equation (24.c)

vanish (Case H.II). The results obtained in the previous three sections suggest that for any reactive system, properties that can be cast as a functional of the electron density, exhibit stationary points along the IRC path defined by the corresponding stationary behavior of the electron density, in keeping with the first theorem introduced by Hohenberg and Kohn17. In particular, by formally studying the properties of the integral -

A %(();+) $ %+ , 0 1 , A

we are able to provide a unified explanation for

the existence of extreme values in different information theoretical measurements as well as the chemical hardness along the reaction path.

Methods and Computational Details


15


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 55

In order to numerically study the nature of the stationary points of the total energy, Shannon information entropy, the disequilibrium, and the hardness along the reaction path, DFT calculations were performed to compute the fully optimized geometries of reactants, transition structures and products of the following three model reactions: a) H ' + H2 H'-H + H (thermoneutral hydrogen abstraction reaction); b) H ' + CH4 H'-H + CH3 (non thermo-neutral -

-

hydrogen abstraction reaction), and c) H + CH4 CH4 + H (thermo-neutral SN2 reaction). All electronic structure calculations were carried out using a locally modified version of the Gaussian 09 (Revision A.02) suite of programs64,65. DFT calculations were performed with Becke’s three-parameter hybrid exchange functional66 and the correlation functional of Lee, Yang and Parr67 using the 6-311++G** basis set. The nature of the stationary points was confirmed by analytical frequency calculations on the optimized geometries using the same level of theory (one negative Hessian eigenvalue for TSs and all positive Hessian eigenvalues for minima). IRC calculations for the three reactions were also performed at the same level of theory using the second order algorithm developed by Gonzalez and Schlegel68 with a step size of 0.05 a0 u1/2. For each reaction, the hardness reactivity index was computed along the IRC using Equation (19) above. The Kohn-Sham orbitals obtained by the DFT calculations were used to construct the electronic density and obtain the Shannon information entropy and the disequilibrium along the IRC by performing numerical 3D integrations used by the DFT code as implemented in the Gaussian 09 program. Finally, the same 3D integrator was also used to compute the term -

A %(();*) $ %+ , A +

4

0 1 along the reaction path. To that effect, we make use of

following relation to compute the derivative of the electron density with respect to the reaction coordinate:


16

Page 17 of 55

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


%(();*)

$

%+

,

+4

= $

%(();*) , %e +

4

%e

$ %+ , .

(25)

+4

where the vector #e = If , fK , ⋯ , f1 N contains the Cartesian coordinates of all N atoms in the reacting system, and the rest of the terms have the same meaning as before. Given the basis set, a molecular geometry, and the corresponding molecular orbitals, the derivative $ e %

%(();*) , e % +

4

can be

obtained analytically, while the term $ %+ , is numerically computed by finite differences. +4

Results Thermo-neutral Hydrogen Abstraction Reaction One of the simplest radical reactions involves the thermo-neutral hydrogen abstraction from a hydrogen molecule by a hydrogen radical: H′ + H-H Η′-H + H . This reaction has been extensively studied using different levels of theory. The classical barrier of Figure 2.a. is approx. 18 kJ/mol, much lower than the MP2/6-31G* (74.06 kJ/mol) and CCSD(T)/6-31G* (62.34 kJ/mol) values previously reported in the literature69. A classical barrier of 47.2 kJ/mol was also found in the work of Esquivel using a CCSD(T)/6-311++G**. In addition, high-level MonteCarlo calculations predict a barrier of approx. 40 kJ/mol70. As expected for this thermo-neutral reaction, the potential energy profile is rather symmetric, with two identical minima (H′ + H-H and Η′-H + H ) connected through a symmetric transition structure (H ′---H---H). Figures 2.b and 2.c depict plots of -

e;s) A 0( 3 e $ ,0 A 0

and $: , $ ;

%233 (+) %+

, as functions of s, respectively (the

chemical potential 7+ was computed by Equation 22). We compute the derivative $

%233 (+) %+

, by

performing forward finite differences of the nuclear-nuclear repulsion energy given by Equation (3) above. Based on our previous discussion, it is clear that this reaction conforms with the


17


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 18 of 55

stationary conditions for Case I reactions, where at the TS (s = 0), both terms, +∞ %(();*)

7+ -−∞ $

%+

, 0 1 and − $

%233 (+) %+

, are equal to zero (see Fig. 2.b and 2.c). As discussed in the

Theoretical Basis section, given that the abstraction reaction H′ + H-H Η′-H + H exhibits a symmetric TS with a plane of inversion perpendicular to the reaction path, one would expect the integral -

e;s) A 0( 3 e $ 0 , 0 A

to vanish as a result of the derivative $

%(();*) %+

, becoming zero at every

point along the molecular space. In order to assess this possibility, central difference derivatives of the electron density were numerically computed at the TS and the results (not shown in this work) projected into a cube file with Gaussian format for further analysis. It was found that indeed, at each point in the cube, the corresponding value of the derivative was close to zero (within 10-8 accuracy), supporting the conclusion that for this reaction, $

%(();*) %+

, is zero for all

at the TS. Figures 3.a – 3.c show the corresponding profiles for the Shannon information entropy, the disequilibrium, and the global hardness.

The results depicted in Figure 3.a show that the

Shannon information entropy becomes a global minimum at the TS, indicating a minimum delocalization or “spread” of the electron density in the vicinity of the TS. We notice that the shape of this curve is somewhat different than the corresponding plot previously reported by Esquivel and collaborators for this reaction37, where the Shannon information entropy exhibits a local maximum at the TS, and two satellite minima symmetrically surrounding this maximum. The difference in entropy value between the satellite minima and the maximum at the TS is less than 0.02 nats. In the work reported by Esquivel and collaborators37 these satellite minima were proposed to be associated to the formation of a reactant (or product) complex (as has been observed in the study of many abstraction and substitution reactions), a feature that is not found


18

Page 19 of 55

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


in the IRC potential energy profile. In order to check the reliability of our results, we proceeded to compute the Shannon information entropy by a numerical 3D integrator developed in our laboratory based on the approach proposed by Becke71, where molecular functions to be integrated are decomposed into single-center components using a generalization of the Voronoi polyhedral scheme. Each single-center component is integrated in spherical polar coordinates around each nucleus using a combination of the Lebedev’s quadrature for the angular part (in our calculations we use a Lebedev’s quadrature of order 17 with a spherical grid of 434 angular points and 80 radial points on each atom) and the Gauss-Chebyshev of the second kind for the radial part. In addition, similar calculations of the Shannon information entropy were performed using the 3D integrator available in the multifunctional wavefunction analyzer software, Multiwfn v. 3.3.7

64,72

. Finally, integration of the Shannon entropy was also performed using

electron densities corresponding to the points along the IRC and dumped into Gaussian density cubes with a dimension of 110 x 110 x 160 (1,936,000 points). It all these cases, it was found that the resultant Shannon information entropy behaved in exactly the same way as predicted by our Gaussian calculations (Figure 3.a), suggesting that the differences between our results and the results reported by Esquivel et. al might be rooted in the different levels of theory used to compute the Shannon information entropy (B3LYP/6-311++G** in this work and QCISD(T)/6311++G**//MP2/6-311++G** in Esquivel’s work) rather than differences in the integration schemes. A study of the effect of the electron correlation on the prediction of Information Theoretical Indexes is being currently carried out in our Laboratory and the results will be reported in a future publication. The local minimum in the Shannon information Entropy located at the TS is consistent with the Case I stationary condition previously discussed.

Figure 3.b depicts the plot of the


19


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 55

Disequilibrium as a function of the reaction coordinate for the H′ + H-H Η′-H + H reaction. As expected, this plot shows a minimum located at the TS (s = 0), in accord with the stationary condition imposed by Case I, where $ %TU (();+)

$

%+

,

+5

%(();*) %+

,

+5

is zero for all , and consequently,

= 0. A similar trend can be found in the case of the Hardness (Figure 3.c),

where a minimum in Y (least stability towards reaction) is located at the transition structure. Non Thermo-neutral Hydrogen Abstraction Reaction Contrary to the previous case, the reaction H ' + CH4 H'-H + CH3 provides a prototypical example of a non thermo-neutral hydrogen abstraction characterized by a transition structure that does not exhibit a plane of inversion perpendicular to the reaction path. As can be seen from Fig. 4.a, this reaction is found to be slightly endothermic with a sizable energy barrier at the level of theory used in this study (Heat of reaction ≈ 3.44 kJ/mol, and classical barrier ≈ 39.56 kJ/mol, results not corrected by zero-point energy). Figures 4.b and 4.c depict plots of

and : $ ;

e;s) A 0( 3 e $ 0 , 0 A

%233 (+) %+

, as functions of s (the small fluctuations observed around the vicinity of the TS

are due to numerical noise in the calculations). The results of Fig. 4.b show that contrary to the abstraction reaction previously discussed, the integral -

e;s) A 0( 3 e $ 0 , 0 A

does not vanish at the TS,

indicating that the stationary condition of the potential energy at the TS does not conform to a Case I but rather to a Case II, where the terms -

e;s) A 0( $ , A 0 +

4

03 e and : $ ;

%233 (+) %+

, cancel +4

each other. Careful examination of Figures 4.b and 4.c, confirms that this is the case. Notice that despite the fact that the integral -

e;s) A 0( 3 e $ 0 , 0 A

becomes zero at around s = -0.62 a0 u1/2 (see


20

Page 21 of 55

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


Figure 4.b), the term

:;

$

%233 (+) %+

, remains finite at this point, indicating that s = -0.62 a0 u1/2 is

not a stationary point in the potential energy surface, given that $

%&' (();*) %+

, ≠ 0, in agreement

with the corresponding energy profile plotted in Fig. 4.a. The Shannon information entropy profile for this reaction is shown in Figure 5.a, while the corresponding disequilibrium and global hardness profiles are depicted in Figures 5.b and 5.c respectively. The results on Figure 5.a indicate the existence of a local maximum located at s = 0.05 a0 u1/2, surrounded by two local satellite minima asymmetrically located at around s = -0.62 a0 u1/2 and s = 0.25 a0 u1/2 (the difference between the minima and the central maxima is less than 0.01 nats). These stationary points are found to be lower in value compared to the corresponding entropy values at the reactants and products, indicating a larger localization of electron density in the vicinity of the TS. Keeping in mind the relation O(; ) =

(();*)

, it can be concluded that

the minimum in the Shannon information entropy profile observed at s = -0.62 a0 u1/2 (Figure 5.a) is the result of the integral 7+ -

A %(();*) $ %+ , 0 1 becoming A

zero at this point (see Figure 4.b),

in agreement with the stationary condition established by Equation (14).

The other two

stationary points found at s = -0.05 a0 u1/2 and s = 0.25 a0 u1/2 conform to the stationary condition given by Eq. (18), where the integral value of the integral -

A %R();+) $ %+ , 0 1 ≠ A

A %R();+) ln O(; ) $ %+ , 0 1 . A

0, is cancelled by the corresponding

In the case of the Disequilibrium (Fig. 5.b), a

local minimum (minimum localization of electron density) is observed located at around s = 0.27 a0 u1/2 and two maxima (maximum localization of electron density) are found at s = -0.58 a0 u1/2 and s = 0.92 a0 u1/2. The results shown in Figure 5.b confirm the lack of existence of a stationary point for S) (; ) at the TS, in agreement with the fact that this abstraction reaction


21


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

belongs to a Case II, where $

%R();+) %+

Page 22 of 55

, is expected to be finite at least partially along the whole +∗

molecular space spanned by . Furthermore, the existence of the three observed stationary points in S) (; ) indicates that it is the expectation value 〈$

%R();+) %+

,〉 that becomes zero at these

points, rather that the derivative itself, in accord with Equation (20). Figure 5.c shows the plot of Y vs s for this reaction. As can be observed in this plot, the hardness becomes a minimum very close to the TS (at s ≈ -0.02 a0 u1/2), indicating that, as in Figure 3c, the reacting system becomes pretty “soft” in this region, or in other words, the reacting system becomes less stable towards reactivity. The fact that the stationary condition for the hardness is not located exactly at the TS is once again consistent with the fact that this reaction belongs to a Case II as previously discussed. Furthermore, according to Eq. (24c), given that -

A %(();+) $ %+ , 0 1 ≠ A

0 at s = -0.02

a0 u1/2, indicates that the stationary condition in Y at this point could only happen if the following condition is satisfied: K

(7 + 7

− 2 7 ) = 0.

(26)

Use of the numerical values for 7 , 7 , and 7 , computed at s = -0.02 a0 u1/2, confirms this conclusion. SN2 Reaction -

-

The reaction H + CH4 CH4 + H constitutes one of the simplest examples of a reaction governed by a bimolecular, nucleophilic substitution mechanism, SN2, where one bond is broken and another one is formed in a synchronous manner, while an electron is transferred from the -

nucleophile group (H ) to the leaving group (H). Figure 6.a shows the corresponding energy


22

Page 23 of 55

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


profile for this reaction. As in the previous reaction, the computed classical barrier (20 kJ/mol) is lower than the calculated by Esquivel and coworkers,37 consistent with the fact the DFT tends to underestimate reaction barriers when compared to wavefunction based ab initio correlated methods. As in the case of the hydrogen abstraction H ' + H2 H'-H + H , the minimum energy -

path for this thermo-neutral reaction is highly symmetric with two identical minima (H + CH4 -

-

and CH4 + H ) connected through a symmetric transition structure [H---CH3---H] . Based on our previous discussions, one would expect this thermo-neutral hydrogen abstraction reaction to conform to a Case I, where the saddle point in the energy profile is the result of a stationary condition in the electron density (and consequently in VNN) at the TS (see Equation (8) above). The results shown in Figures 6.b and 6.c confirm this conclusion, and show that both, the integral -

e;s) A 0( 3 e $ 0 , 0 A

and the derivative

%233 %+

become zero at the TS (s = 0). Inspection of these two

plots also reveals that despite the fact that -

e;s) A 0( 3 e $ ,0 A 0

and

%233 %+

additional points (s = -1.46 a0 u1/2 and s = 1.41 a0 u1/2 in the case of 1.42 a0 u1/2 and s = 1.38 a0 u1/2 in the case of

%233 %+

%233 %+

e;s) A 0( 3 e, $ ,0 A 0

and s = -

), no additional stationary conditions in

(; ) are found at these points, given that whenever the derivative

become zero at two

e;s) A 0( 3 e = $ 0 , 0 A

0, the value of

is finite (and vice versa).

Figure 7.a depicts a plot of the Shannon information entropy as function of the reaction coordinate for this reaction. According to the results shown on this figure, a maximum is observed at the TS, which results from the integral -

e;s) A 0( 3 e $ ,0 A 0

becoming zero. In contrast to

the hydrogen abstraction reactions discussed above, this SN2 mechanism exhibits a maximum in the Shannon information entropy with a relatively larger value at the TS compared to the


23


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 55

corresponding values at the reactants and products. This maximum is bracketed by two minima symmetrically located at s = -1.85 a0 u1/2 and s = 1.85 a0 u1/2 respectively, with their corresponding stationary points defined by the stationary conditions established by Equation (15). A qualitatively similar profile for the Shannon entropy was found by Esquivel et al.37 The results in Figure 7.b show that the TS exhibits a maximum in the disequilibrium, S) (), bracketed by two other maxima symmetrically located at s = -1.5 a0 u1/2 and s = 1.5 a0 u1/2. These maxima indicate a relatively large degree of electron density concentration. The fact that a maximum in the Shannon information entropy is also observed at the TS suggests that the negative charge in the reacting system is mainly concentrated along the reaction coordinate with a significant spread along the H---C---H bond in the vicinity of the transition structure. As with any other symmetric thermo-neutral reaction conforming to Case I, the maximum in S) () lm

located at the TS is due to the stationary condition in the electron density ($ , ln

+5

= 0) as

discussed before. The other stationary points in the disequilibrium plot are the result of the expectation value 〈$

%R();+) %+

,〉 becoming zero. Finally, the plot depicted in Figure 7.c exhibits a

maximum in the hardness at the TS bracketed by two minima symmetrically located at approx. at s = -1.5 a0 u1/2 and s = 1.5 a0 u1/2. The maximum in the vicinity of the TS is once again consistent lo

with a Case I reaction mechanism, where $ ln ,

+5

= 0 and it indicates a relatively maximum

stability towards reactivity. Given that the value of the integral -

A %(();+) $ %+ , A +

4

0 1 evaluated at

the minima is finite, it indicates that these minima are the result of the stationary condition imposed by Equation (26). Conclusions


24

Page 25 of 55

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


The Kohn-Sham DFT formalism has been used to investigate the influence of the stationary behavior of the electron density along a minimum energy path on the corresponding stationary conditions observed in the total potential energy of the reactive system, information theory measures (Shannon information entropy and Onicescu information energy), and chemical reactivity indexes (the chemical hardness). The theoretical treatment presented in this work indicates that in the case of thermo-neutral reactions characterized by symmetric transition structures exhibiting an inversion plane perpendicular to the reaction path, an extremal in the potential energy (as well as for the other properties considered in this study) at this point results from a corresponding stationary condition in the electron density throughout the whole molecular space, $

%(();*) %+

,

+5

= 0, (Case I). For other cases (Case II), no such stationary behavior in

(; B ) is observed, and instead, the stationary behavior of the total potential energy at the TS results from the cancellation of the integrals - $

%(();*) %+

, 0 1 and 9 +4

:; 4

Stationary Conditions of the Electron Density Along ... - ACS Publications

Recommend Documents