Atomic decomposition of conceptual DFT descriptors: application to proton transfer reactions.

PCCP View Article Online

Published on 19 June 2015. Downloaded by Nanyang Technological University on 19/06/2015 14:21:34.

PAPER

Cite this: DOI: 10.1039/c5cp01515d

View Journal

Atomic decomposition of conceptual DFT descriptors: application to proton transfer reactions Ricardo Inostroza-Rivera,a Meziane Yahia-Ouahmed,b Vincent Tognetti,*b a Laurent Joubert,b Ba´rbara Herrera*a and Alejandro Toro-Labbe ´ In this study, we present an atomic decomposition, in principle exact, at any point on a given reaction

Received 15th March 2015, Accepted 26th May 2015 DOI: 10.1039/c5cp01515d

path, of the molecular energy, reaction force and reaction flux, which is based on Bader’s atoms-inmolecules theory and on Penda ´s’ interacting quantum atoms scheme. This decomposition enables the assessment of the importance and the contribution of each atom or molecular group to these global properties, and may cast the light on the physical factors governing bond formation or bond breaking.

www.rsc.org/pccp

The potential use of this partition is finally illustrated by proton transfers in model biological systems.

1. Introduction The course of a chemical reaction, which can be considered as a very complex process at the molecular scale, can be described, within the Born–Oppenheimer approximation and the static approach of reactivity, by the molecular energy profile Emol(x), where the reaction coordinate x parameterizes the reaction path on the potential energy surface (PES) and notably locates the relevant points (reactants, products, transition states. . .) along it.1 From this profile, other meaningful quantities can be extracted, which may provide deeper understanding of the reaction process and of the physico-chemical factors that drive it. Among them are the so-called reaction force and reaction ´ flux, which were derived and extensively studied by Toro-Labbe and co-workers.2–15 The first one relies on the fact that the general shape of the reaction profile for any individual chemical step is universal and enables us to divide the energy profile into significant regions associated to well-characterized physical events. The second one is based on conceptual density functional theory (CDFT),16,17 which is a theoretical interpretative subfield where relevant chemical information is extracted from the electron density, and which has been also proven to be a quantitative predictive theory. We will in particular focus on the electronic chemical potential, which is a molecular property that is

actually18 the opposite of Mulliken’s electronegativity in terms of molecular electron affinity and ionization potential19 when the number of electrons is varied by one unit. Its evolution along the reaction path is thus expected to unravel the electronic activity inside the systems studied. All these mentioned descriptors are global ones, in the sense that they do not depend on real space position. They accordingly do not give insight into local behavior and regiospecificity. However, in general, chemists focus on the role of functional groups or moieties, they seek to deduce their role in the chemical properties of a molecular species and are used to determine which atoms contribute the most to the ‘‘chemistry’’ they develop.

a

Nucleus Millenium Chemical Processes and Catalysis (CPC), Laboratorio de Quı´mica Teo´rica Computacional (QTC), Facultad de Quı´mica, Pontificia Universidad Cato´lica de Chile, Santiago, Chile. E-mail: [email protected] b Normandy Univ., COBRA UMR 6014 & FR 3038, Universite´ de Rouen, INSA Rouen, CNRS, 1 rue Tesnie´re 76821 Mont St Aignan, Cedex, France. E-mail: [email protected]

This journal is © the Owner Societies 2015

Scheme 1 Chemical equations for the proton transfers in formamide (up) and 2-pyridone (down). Atoms whose properties are followed along the reaction path are coloured in orange.

Phys. Chem. Chem. Phys.

View Article Online


Paper

PCCP

In addition, experimentalists usually test various substituents and spend considerable time fine-tuning them in order to improve their processes. The primary aim of this study is to present an additive atomic decomposition of such global descriptors, which can translate the information conveyed by these theoretical tools into the more common language of reactive sites and functional groups effects. Such a decomposition scheme will be built on Bader’s atoms-in-molecules theory (QTAIM),20,21 which shares the same fundamental ingredient as DFT and CDFT, namely, the electron density, and provides an unambiguous and physically rooted partition of a molecule into physical atoms. As an example of its application, it will be used to investigate proton transfer in two model systems: formamide and 2-pyridone (as represented in Scheme 1), because proton transfer is an ubiquitous fundamental process, particularly in biological systems.22–29 In addition, both molecules possess an amide functional group, which is the basic building block of proteins and enzymes, and they can also be considered as models for nucleic bases. Furthermore, a comparison between formamide and 2-pyridone might disclose information on the effect of aromatic delocalization.

2. Atomic decomposition of global descriptors We first mathematically define the global descriptors that will be decomposed into atomic contributions. The reaction force is the opposite of the derivative of the molecular energy with respect to the reaction coordinate:2 FðxÞ ¼

dEmol ðxÞ : dx

(1)

In the energy canonical representation within CDFT, one inspects the partial derivatives of the electronic energy (at variance with the reaction force that involves the molecular energy) with respect to the total number of electrons N and the external potential (generated by the nuclei) n(r). The electronic chemical potential is then defined as follows: @Eel ðxÞ mðxÞ ¼ : (2) @N vð~ rÞ

To distribute these quantities over all atoms, a partition of the molecule itself is required. As stated in the introduction, we will advocate the use of QTAIM, which divides the molecule into non-overlapping atomic basins (domains in real space) separated by interatomic surfaces, according to the electron density gradient vector field. As for the molecular energy, an additive atomic decomposition already exists, based on the virial theorem (where the Xa denote the Cartesian nuclei coordinates):31,32 Emol ¼ Tel

X

Xa

a

@Emol : @Xa

(4)

At a stationary point on the PES (the reactants, products and transition states), the Hellmann–Feynman forces acting on the nuclei all vanish, so that (since the electronic kinetic energy operator is monoelectronic): X Emol ¼ Tel ¼ TðAÞ; (5) A

where the atomic kinetic energy T(A) is well defined. On the contrary, when the system is out of equilibrium, the partition of the molecular energy into atomic components becomes a less trivial issue in this framework, because the second term in eqn (4) also depends itself on the total molecular energy. As we are interested in following the evolution of atomic properties along a whole reaction path, such a virial approach seems to be problematic. As an alternative, one can use the interacting quantum ´s and atoms (IQA) energy decomposition developed by Penda co-workers,33–36 which is an unarbitrary way of partitioning the energy of a molecule into intra-atomic energies E self IQA(A) and inter-atomic terms, these last ones exactly reducing to pair interactions, E inter IQA (A, B): Emol ¼

X

self EIQA ðAÞ þ

A

1 X X inter E ðA; BÞ: 2 A BaA IQA

(6)

! 1 X inter þ E ðA; BÞ : 2 BaA IQA

(7)

1 X inter E ðA; BÞ; 2 BaA IQA

(8)

Eqn (6) is equivalent to: Emol ¼

X

self EIQA ðAÞ

A

If we define: The reaction electronic flux (REF) was designed in analogy with thermodynamics according to7 JðxÞ ¼

dmðxÞ : dx

(3)

It offers a method of characterizing the electronic activity taking place during the chemical reaction. Positive values of the REF indicate a spontaneous electronic activity, in other words, the reaction is mainly driven by bond formation/strengthening. Negative values of the REF, conversely, reflect a nonspontaneous electronic activity, which means that the reaction is mainly driven by bond breaking/weakening.7–9,11,15,30


self ðAÞ þ EIQA ðAÞ ¼ EIQA

we get 37 Emol ¼

X

EIQA ðAÞ:

(9)

A

Importantly, this partition is always valid whether the system is taken at a PES stationary point or not. (However, note that at a PES stationary point, it will not in general give the same values as the virial approach).


View Article Online

PCCP

Paper

One can now extend this procedure to other related quantities. This is straightforward for the reaction force: X dEIQA ðAÞ : (10) FðxÞ ¼ dx A |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl}


FðAÞ

For the chemical potential, we first use finite difference linearization (FDL)38 to write: EelNþ1 ðxÞvð~rÞ EelN1 ðxÞvð~rÞ ; (11) mFDL ðxÞ ¼ 2 EN1 el

denotes the electronic energy of the examined where molecule with N 1 electrons. As this energy difference is evaluated at constant external potential (geometrical relaxation is not considered upon addition/removal of electrons), eqn (11) is equivalent to 1 0 Nþ1 N1 E ðAÞ E ðAÞ IQA X B IQA vð~ rÞ vð~ rÞ C FDL m ðxÞ ¼ (12) A: @ 2 A |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} mx ðAÞ

With eqn (12), we have thus assigned a contribution for each atom to the total molecular chemical potential. Hence X dmx ðAÞ : (13) JðxÞ ¼ dx A |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}

intermolecular ones on the other hand, we can define the atomic mintra, inter(A) quantities for the electronic chemical potential (evidently, this can be done similarly for the energy, the reaction force and the reaction electronic flux): 0 self Nþ1 self N1 E ðAÞ E ðAÞ B IQA X B IQA vð~ rÞ vð~ rÞ FDL m ðxÞ ¼ B @ 2 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} A mintra ðAÞ

þ

inter Nþ1 ðA; BÞ X EIQA

vð~ rÞ

inter N1 EIQA ðA; BÞ

1 C

vð~ rÞ C C

C 4 A |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

:

BaA

minter ðAÞ

(14) ¨rling–Levy perturbational theory,39 to Finally, it is possible, using Go exactly decompose the interaction energy into classical (gathering the nuclei-nuclei repulsion, the attraction energy of the electrons by the nuclei and the pure Coulombic part of the electron-electron repulsion) and exchange–correlation contribution (collecting all ‘‘quantum effects’’): cl xc Einter IQA (A, B) = EIQA(A, B) + EIQA(A, B).

(15)

Therefore, minter(A) (and correlatively, the reaction electronic flux) can be split into classical and exchange–correlation components, denoted as minter cl(A) and minter xc(A), respectively:

JðAÞ

This FDL approach also calls for some theoretical comments. Consider for instance a stationary point (at xSP) on the PES (reactants, products, or transition states) for the N-electrons species. In general, the related geometry will not correspond to a stationary point on the PES for the same molecule with N 1 electrons. As a consequence, even at these xSP points, the use of the virial approach (eqn (5)) to evaluate the atomic energy of these last species would still remain ambiguous, at variance with our proposed approach based on the IQA decomposition. dE N1 Furthermore, as mol ðxSP Þa0, J(xSP) may not vanish in general, dx even at the beginning or at the end of the reaction studied. It can also be added that eqn (12) can be furthermore directly translated in terms of the molecular electronegativity (that is the opposite of the chemical potential in Mulliken’s definition), which can also be accordingly decomposed in terms of atomic contributions. This would enable us to identify, inside a molecule, which atom or which substituent group (note that, as Bader’s atoms do not overlap, a moiety quantity is simply the sum of those for the atoms that constitute the functional group) is responsible for the rather high or low electronegativity of the whole molecular system studied. In addition to partitioning global molecular quantities in terms of atomic contributions, we can explore each mono-atomic contribution in terms of its intra-atomic and inter-atomic com1 P inter self ponents. Replacing EIQA(A) by EIQA ðAÞ þ E ðA; BÞ, 2 BaA IQA and gathering together self-contributions on one hand and


minter(A) = minter cl(A) + minter xc(A).

(16)

In other words, the interatomic value of all these descriptors can be traced back to their electrostatic and covalent origins, which will enable us to assess to what extent the nature of given interactions impact on the global properties. It is worthy noticing that this scheme is general and can be applied to any global descriptor that depends on the energy in a linear way. For instance, one can use finite difference linearization to decompose the molecular hardness (that is the second derivative of the energy with respect to the electron number) in atomic contributions, which could be useful in optoelectronics since the molecular hardness is closely related to the molecular gap: X Nþ1 N1 N ZFDL ðxÞ ¼ EIQA ðAÞ þEIQA ðAÞ 2EIQA ðAÞ : vð~ rÞ vð~ rÞ vð~ rÞ A |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ZðAÞ

(17)

3. Computational details All stationary points were fully optimized at the DFT level using the B3LYP hybrid exchange–correlation functional with the 6-311++G(d,p) basis set and the Gaussian 09 program.40 In addition, frequency calculations were carried out to check their nature. The minimum energy reaction pathways connecting the reactants and the products through the transition state (TS) were computed at the same level of theory using the intrinsic reaction coordinate (IRC) formalism40–43 (x will be reported in


View Article Online

Paper


atomic units throughout this study). This IRC procedure gave a set of 108 geometries for formamide and 117 points for 2-pyridone. Then, a full QTAIM/IQA treatment (featuring topological analysis as well as energy decomposition) was performed on each structure using the AIMAll software.44 Let us recall that ~ rÞ ¼ ~ critical points are points in real space where rrð~ 0, and that they can be classified according their rank and signature. When r2r(r) is positive (or negative), the electron density is locally depleted (or concentrated). Such a local analysis is complemented by the so-called integrated one (basin populations for instance) and by the IQA treatment. Note that for a detailed discussion about the use of IQA within the Kohn–Sham (in KS) approach, we refer the interested reader to our recent studies.45–48 It must be remarked that in KS-IQA, the IQA atomic energies do not exactly add up to the total molecular energy (the difference being very small). We thus decided to uniformly scale our original IQA atomic energies to make them exactly add up to the SCF energy given by Gaussian 09. More details about this procedure are given in the Appendix section. Note that the very same trick has been routinely used for decades in the DFT virial approach (since KS does not afford the true kinetic energy, but only that of the noninteracting system).

PCCP

Finally, we recall that the reaction force enables us to divide the reaction path into three regions.49,50 Indeed, for any elementary chemical step, F(x) starts from 0 (reactants at xR), goes through a minimum located at x1 o 0, is once again equal to 0 at the TS, and then becomes positive, reaching a maximum at x2 4 0, before vanishing at the product (located by xP). In all of the presented graphs, two vertical red lines will indicate the positions of x1 and x2. The x1 o x o x2 range can be properly coined the ‘‘transition state region’’ from the reaction force point of view. The question whether it coincides with such a region from a QTAIM perspective will be also addressed in the following section.

4. Results and discussion With these tools in hand, we now focus our attention on the proton transfer in the formamide and 2-pyridone molecules (Scheme 1). Note that we arbitrarily chose to consider the endothermic reaction where the transferred proton moves from the nitrogen atom (negative values for x) to the oxygen atom (positive values). Of course, the REF sign does not depend on the retained direction because it has an intrinsic physical meaning (a bond that would be considered as stretching from a x point of view would be contracting from the x one).

Fig. 1 Electron density (upper part) and laplacian (lower part) profiles at the selected (3,1) and (3,+1) critical points with respect to the reaction coordinate describing proton transfer in formamide and 2-pyridone. Red lines delimit the transition state region as defined by the reaction force profile.



View Article Online

PCCP


4.1. Atomic interactions from the QTAIM point of view: electron density topology and delocalization indexes We will first have a brief look at the evolution (in the context of catastroph theory) of the electron density topology along the reaction path, as done for instance in ref. 51–53. For both compounds, the same behavior was in fact observed. At the beginning, the transferred proton is only bonded to nitrogen, one (3,1) critical point (CP) being located between the two nuclei. The r(3,1)(N–H) electron density first slowly decreases (see the upper part of Fig. 1). However, when arriving inside the TS area (as defined by the reaction force profile), it drops significantly. However, before it vanishes, another (3,1) CP is created between the hydrogen atom and the oxygen. The electron density at this second CP rapidly increases until the system gets out of the transition state area (where r(3,1)(O–H) becomes more and more constant), where the N–H (3,1) CP is being annihilated. Thus, inside the TS region, both (N–H) and (H–O) CPs ´–Hopf relaco-exist. Then, as the consequence of the Poincare tionship, one (3,+1) CP (typical of a topological ring) must be simultaneously present. The corresponding density at this (3,+1) CP begins to increase until reaching the TS where it takes its maximal value, before decreasing until annihilation at the frontier of the TS zone. The lower part of Fig. 1 now shows the evolution of the electron density laplacian at the mentioned CPs, which serves as a way of characterizing the nature of the formed bonds. r2r(3,1)(N–H) is first negative until the system reaches the TS, which is in agreement with the shared interaction nature (i.e. covalent) of the N–H interaction. At the TS, r2r(3,1)(N–H) becomes positive, which means the nature of the N–H interaction has changed to a closed shell (i.e. mainly electrostatic (ionic) one). As intuitively expected, the exact opposite behavior is observed for r2r(3,1)(H–O), and as already observed for r(3,+1), r2r(3,+1) is maximal at the transition state. It should be added that the TS is also characterized by the concomitant crossing of the (3,1) CP density and laplacian curves for the forming and breaking bonds: such a topological

Paper

characterization in terms of crossing point is an appealing feature, which is complementary to the pure energetic one. Lastly, it is an important result that the ‘‘transition state region’’, as defined by the reaction force (marked by the vertical red lines in Fig. 1), almost corresponds to that of the presence of the (3,+1) critical point (delimited by the apparition and annihilation of (3,1) CPs). The fact that the reaction force theory and QTAIM provides a consistent definition of this zone is, from our viewpoint, remarkable. Such a conclusion is also fully supported by the inspection of the QTAIM delocalization indexes (DIs, which measure the number of electron pairs delocalized between two atoms) that serve as a type of bond order,54–56 and whose shape as a function of the internuclear distance clearly reflects the interaction nature (exponential for repulsive/non bonded interactions, sigmoidal for formation/breaking of chemical links).57 They are depicted in Fig. 2. As expected, they only vary slightly in the reactants and the products regions, conversely exhibiting large variations in the TS region. Interestingly, they provide a new QTAIM characterization of the transition state: for the two systems, it reveals to be the point where the delocalization indexes for the N–H and O–H bonds (note that they are really bonds from Bader’s orthodox viewpoint since these two (3,1) N–H and O–H CPs are exactly bond critical points at a PES stationary point) intersect. 4.2.

Energy profiles

We now look at the relative energy profiles. Fig. 3 shows, for both systems, the DFT molecular energy profile (dashed curve) as well as the contribution of each of the four main atoms (transferred proton, O, N and the carbon atom bonded to O and N) to the DFT energy, according to eqn (9). It immediately appears that the carbon contribution is negligible, contrary to those for O, H and N (one can notice that the curves shapes are similar for both systems). Interestingly, inside the reactants and products regions, all these atomic components (with the possible exception of oxygen) vary only slightly, whereas inside the TS area, large variations

Fig. 2 Delocalization indices profiles with respect to the reaction coordinate describing proton transfer in formamide and 2-pyridone. Red lines delimit the transition state region as defined by the reaction force profile.



View Article Online


Paper

PCCP

Fig. 3 Molecular and atomic energetic profiles with respect to the reaction coordinate describing proton transfer in formamide and 2-pyridone. Red lines delimit the transition state region as defined by the reaction force profile.

are observed. More specifically, the H and N energies increase and thus raise the activation barrier (destabilizing contributions), whereas oxygen stabilizes the systems once the transfer is initiated, because its energy drops to about a very negative value of 43 kcal mol1. It is noteworthy that the hydrogen energy curve goes through a maximum value at the TS point. It is therefore possible to define an atomic activation energy for this atom: we found it is equal to 35 kcal mol1 in formamide and 32 kcal mol1 in 2-pyridone, representing about 70% of the total molecular activation barrier. This primacy of hydrogen atom is fully consistent with the preponderance observed in a previous study ¨hringer-Martinez and Toro-Labbe ´.58 by Vo One can go further by scrutinizing the decomposition in terms of intra-atomic and inter-atomic terms as shown in Fig. 4. For the carbon atom (upper part of Fig. 4), the interatomic component (Einter) as well as the intra-atomic one (Eintra) contribute to the same extent to the overall atomic energy (about 100 kcal mol1 in absolute value), but they are of opposite signs and thus counterbalance each other; Eintra is stabilizing, whereas Einter has a destabilizing effect. This can be accounted for by following Bader’s atomic charges along the reaction path. Indeed, as shown in Fig. 5, the one for carbon decreases, so that the carbon atom gains electrons (about 0.20 e in both cases), making the intra-atomic contribution more and more negative. We note that this energetic stabilization is logically higher in the aromatic system because the electron gain is also slightly more important. The situation is exactly opposite for the hydrogen atom: Einter is stabilizing, whereas Eintra is destabilizing. This last one is thus the main factor responsible for the existence of the aforementioned atomic activation barrier for hydrogen. Once more, it can be explained in terms of the electronic population. Indeed, along the proton transfer, its charge increases, so that there are less electrons inside the hydrogen basin, and so its self energy is less negative (basin less stabilized). These intra-atomic and inter-atomic contributions can also be summed up to reveal the inter atomic and intra atomic


contributions at the molecular level: the already observed trend inter still holds true: Eintra mol and Emol tend to have opposite signs along the whole reaction path. 4.3.

The reaction force

As shown in Fig. 6, the reaction force profile is characterized by the two already mentioned extremal points (denoted as x1 and x2) corresponding to a minimum and a maximum for the reaction force, respectively. Previous studies gave evidence that in each of these regions, well-defined and differentiated chemical events tend to prevail. For instance, the reactants region (xR o x o x1) mainly involves structural rearrangements (such as rotational, conformational changes. . .) that prepare the transformation itself. In the TS region (x1 o x o x2), intense electronic reordering occurs with bond formation and bond breaking, so that most of the electronic activity takes place inside this range. Lastly, in the products region (i.e. x2 o x o xP), structural arrangements occur again allowing relaxation towards the final geometry of the products. By analogy with classical mechanics, such evolutions along the reaction path can be quantified by evaluating the reaction force work between two given points xa and xb, which can be given by the following integral: ð xb (18) Wab ¼ FðxÞdx ¼ Emol ðxa Þ Emol ðxb Þ: xa

Four particular works, associated to each substep of the considered elementary chemical step, can thus be defined (using obvious notations): WR1, W1TS, WTS2, W2P, which can all be decomposed into atomic contributions using eqn (9) in order to better grasp their physical nature. Interestingly, they are linked to the direct activation barrier and the reaction energy by: DE act ¼ WR1 þ W1TS : (19) DE react ¼ WR1 þ W1TS þ WTS2 þ W2P From Table 1, we can infer that the structural work WR1 is the main contribution to the molecular activation energy DEact (since it represents 67% and 60% of DEact in formamide and


View Article Online


PCCP

Paper

Fig. 4 Intra-atomic and inter-atomic component profiles for the C and H atomic energetic profiles with respect to the reaction coordinate describing proton transfer in formamide and 2-pyridone. Red lines delimit the transition state region as defined by the reaction force profile.

Fig. 5 Atomic charges profiles with respect to the reaction coordinate describing proton transfer in formamide and 2-pyridone. Red lines delimit the transition state region as defined by the reaction force profile.

in 2-pyridone, respectively), hinting that structural rearrangements drive the proton transfer, whatever be the direction. Indeed, in the reverse direction, W2P is about twice WTS2. However, from an atomic point of view, the situation is more intricate. For instance, in the forward direction, WR1 seems to be mainly due to the hydrogen atom (with a non-negligible difference between the two molecules: 64% in formamide and


74% in the 2-pyridone). On the contrary, there is no more discrepancy between the two systems for W2P, and the hydrogen contribution now only amounts to about 35%. In the transition state region, nitrogen and oxygen play opposite roles: the N contribution being positive, whereas the O one is negative for both molecules (as also shown in Fig. 6). It is also instructive to look at the W1TS + WTS2 sum: for instance, for formamide,


View Article Online


Paper

PCCP

Fig. 6 Molecular and atomic reaction force profiles with respect to the reaction coordinate describing proton transfer in formamide and 2-pyridone. Red lines delimit the transition state region as defined by the reaction force profile.

Table 1 Atomic decomposition of the reaction force works. All values in kcal mol1

Reaction

Atom

DEact DEreact WR1 W1TS

WTS2

W2R

Formamide C1 2.3 1.1 3.1 0.8 0.4 0.8 N2 24.6 37.9 7.6 17.0 11.6 1.8 O4 11.4 43.3 0.8 12.1 16.6 15.3 H6 34.9 24.2 20.0 14.0 2.7 8.0 C + N + O + H 45.7 17.7 26.1 19.7 7.4 20.7 Molecule 48.5 13.0 32.7 15.9 12.1 23.4 2-Pyridone

C5 N11 O10 H12 C + N + O+H Molecule

0.2 0.0 2.3 21.5 30.8 5.1 8.9 43.0 1.0 31.2 19.7 17.1 43.6 7.4 21.0 38.0 0.9 22.9

2.1 0.2 0.1 16.4 10.4 1.1 9.9 19.2 15.0 14.0 3.3 8.2 22.7 11.9 24.3 15.1 14.0 23.0

it equals 28.5 kcal mol1 (N) and 28.7 kcal mol1 (O), to be compared with the N and O contributions to the total reaction energy: 37.9 kcal mol1 (N) and 30.8 kcal mol1 (O). This means that both atoms are mainly involved in strong electronic activity.

4.4.

The reaction electronic flux

The REF profiles are shown in Fig. 7. We recall that the REF measures the (opposite) variation rate of the chemical potential along the reaction path, and therefore the variation of the molecular electronegativity. In the case of formamide, we clearly see that the molecular REF (dashed curve) increases from negative values to positive ones, the change in sign occurring at the transition state. This means that the reaction is first driven by a non-spontaneous event and that bond weakening/breaking initially dominates. At the TS, both processes are balanced so the total molecular reaction electronic flux vanishes. Then, bond strengthening/ forming overcomes bond weakening/breaking until the end of the reaction. Note that as the products and reactants geometries do not in general correspond to stationary points on the PESs of the ionic N 1 species, the molecular REF does not actually vanish at these points as already discussed at the end of Section 2. The REF atomic components for formamide can now be considered. In the reactants region, the carbon REF, as well as the nitrogen REF, is negative, while the proton’s REF is positive,

Fig. 7 Molecular and atomic reaction electronic flux profiles with respect to the reaction coordinate describing proton transfer in formamide and 2-pyridone. Red lines delimit the transition state region as defined by the reaction force profile.



View Article Online


PCCP

thus opposing to the general trend. To interpret these signs, one must look at the evolution of each bond inside this region: the N–H bond starts breaking as the proton leaves the nitrogen atom, which weakens the CQO bond while the O–H interaction increases in absolute value. Thus, for each atom, there is a competition between a forming bond and a breaking one. The signs that we obtained for the atomic REF suggest that the O–H strengthening overcomes the N–H weakening, and that the CQO weakening process (from double to single bond in the pure Lewis picture) overcomes the strengthening of the CQN bond (from a single to a double bond in the idealized Lewis view). These facts also explain the negative values for nitrogen’s REF in this reactants region. In the TS region, where most of the electronic activity takes place, we observe two large contributions to the total REF, stemming from the N and O basins, which is concordant with the previously noted fact this is the region where the work of the reaction force is the most important for these two atoms (with very large absolute values for the atomic W1TS and WTS2). The nitrogen’s REF is positive, meaning that the CQN bond is spontaneously formed. As for the oxygen REF, it is however negative, meaning that the partial destruction of the CQO bond overcomes the creation of the O–H bond. These opposite behaviors are actually reminiscent of the behaviors observed for the atomic decomposition of the W1TS and WTS2 molecular quantities, since the contributions of N and C atoms are of opposite signs in all cases. Lastly, the carbon REF changes sign from negative to positive, so that, at this point, the strengthening of the CQN bond overcomes the weakening of the CQO bond. Then, in the products region, all atomic components of the REF decrease to low positive values. This means that bond breaking and forming have finished, being replaced by geometrical relaxation effects, and this is in agreement with the usual interpretation of the reaction force. We now study the 2-pyridone system. While no important discrepancies (with respect to formamide) were observed for the energetic and reaction force profiles, the reaction electronic fluxes considerably differ from each other. In that sense, the reaction electronic flux thus acts as a subtle indicator of the electronic effects. For 2-pyridone (Fig. 7), it can be noticed that the molecular REF profile goes from positive values (suggesting predominant bond strengthening/forming before the TS) to negative ones (indicating that bond breaking processes dominate this part of the reaction path). From a structural point of view, the main difference between 2-pyridone and formamide is evidently the presence of the aromatic cycle with delocalized p electrons over the whole cycle. When the proton moves away from the nitrogen atom (negative flux, i.e. overall bond breaking for nitrogen), the latter (although being already negatively charged) does not transfer electronic charge to the carbon because the charges are delocalized over the cycle. This can be confirmed by looking at the charge versus reaction coordinate plots (Fig. 5). At the beginning of the reaction with 2-pyridone (first half of reactants zone), nitrogen gains electrons (either coming from the departing


Paper

proton and/or from the cycle), but the charge of the carbon atom does not vary. As a consequence, the CQO bond is considerably less weakened at the very first moments of the reaction than in formamide. In the TS region, the same mechanism as in formamide occurs. The CQO bond partly breaks and the CQN bond strengthened. The amplitude of the oxygen and the nitrogen fluxes are, however, reduced. Partly breaking the CQO bond is thus a process that is less favorable than reinforcing the CQN bond, giving the negative REF value observed for the carbon atom. Such differences are not unexpected and are due to the extended delocalization on the whole cycle. Once more, there is a pure QTAIM descriptor that is able to recover similar information. Indeed, it must be noticed that the DI between the carbon and nitrogen does not increase in 2-pyridone (from 0.9 to 1.3) as much as that found in formamide (from 1 to 1.6), which reflects the different variations in the atomic REF values. One can even go further in the analysis of the aromaticity changes in the 6-membered ring along the reaction path. To this aim, we inspected the so-called ‘‘para-delocalization index’’ (PDI)59,60 defined as the average value for the three DIs involving atomic pairs in the para position. Interestingly, we found that it features an almost perfect sigmoidal shape according to (where all quantities are expressed in atomic units and with R2 = 0.997): PDIðxÞ ¼ 0:0935

0:0284 : 1 þ e2:330x

(20)

This implies that the aromaticity of the six-membered ring monotonically increases during the proton transfer, but it also provides another valuable characterization of the transition state since it exactly corresponds to an inflexion point (sharp dPDIðxÞ maximum for ) in the PDI profile. It is also noteworthy dx that we have a similar trend for the evolution of the corresponding ring critical point density value (a quantity that has been already used to describe heterocycles61,62) since it is linearly correlated (R2 = 0.99) to PDI. Finally, one can wonder whether this QTAIM analysis has its conceptual DFT counterpart. In particular, it has been early suggested that aromaticity is also related to hardness.63 We actually found that the molecular hardness is higher in the product than in the reactant, which is in line with a higher aromaticity. Interestingly, the transition state does not correspond to a minimal hardness value (neither to a minimal value for the HOMO–LUMO gap) as it has been encountered in several reactions (see ref. 64 for a thorough discussion on this point). Once more, QTAIM and conceptual DFT, both based on the electron density, offer complementary and consistent approaches of the same physicochemical property.

5. Conclusions In this work, we proposed an atomic decomposition, based on Bader’s atoms-in-molecules theory and on the interacting quantum atoms approach, of the molecular energy, the reaction force,


View Article Online


Paper

PCCP

and of any global conceptual density functional theory descriptor (in particular of the reaction electronic flux) at any point on a reaction path. This scheme enables us to follow the contribution of each atom to the reactivity of the whole system along the reaction path. Besides, these atomic components can also be split into intra-atomic contributions and inter-atomic ones, the latter involving electrostatic and exchange–correlation (or covalency) effects. As a proof of concept of this method, we applied it to proton transfers in model systems of biological interest. On these two examples, the complementarity and consistency of the various approaches (reaction force, conceptual DFT, reaction electronic flux, Bader’s atoms-in-molecules theory) were highlighted, and new characterization of the activation barriers (through their decomposition into atomic activation energies) and of the transition states identified. In particular, QTAIM and the reaction force theory provide a consistent definition of the transition state region, embodying the strong relationships between energy and electron density topology. Such links are actually supported both from the mathematical point of view through the known structural homeomorphism65 between the virial and electron density field and from a more physicochemical perspective via the concept of privileged exchange channels.66 These atomic descriptors notably enabled us to show that even a simple proton transfer is actually a more complex process at the atomic scale than suggested by the simple view of atomic moves. Applications to other reactions will be reported in due course.

This decomposition is not implemented in AIMAll. Nevertheless, only in the cases of LDA, B3LYP and M06-2X, the following scheme is available in the latest AIMAll versions when a nondefault flag (‘‘/MODEL’’) is added to the wfx file. By definition, this scheme exactly sums up to the DFT molecular energy and includes full exact non-local exchange for the inter-atomic terms (exactly as in post-HF theories): ! X a 1 X XX hybrid self;simpl Emol ¼ Exc ðAÞ þ E ðA; BÞ 2 BaA x A 1 X X cl EIQA ðA; BÞ þ ExXX ðA; BÞ 2 A BaA

þ

X

E self;modelAIMAll ðAÞ þ

A

(A3) For other functionals (or when the ‘‘/MODEL’’ flag is absent), the atomic self-energy, Eself,defaultAIMAll(A), is approximated by the only atomic classical and XX exchange contributions. As we seek to provide a methodology able to recover the molecular energy for any functional, we proposed the following general scaling strategy: 1X X X hybrid cl Emol ¼ g EIQA ðAÞ þ ExXX ðAÞ þ gE inter;AIMAll ðA;BÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} 2 A BaA A E self;defaultAIMAll ðAÞ

X

E self;thiswork ðAÞ þ

A

1 X X inter;thiswork E ðA;BÞ: 2 A BaA (A4)

Appendix: details on KS IQA schemes In this appendix, we discuss in more detail the main possible decomposition of the molecular energy within the Kohn–Sham framework. For any global hybrid functional characterized by its exact exchange (XX, Hartree–Fock (HF) type) mixing parameter a and the local (in the mathematical sense) exchange– correlation functional (belonging to the LDA, GGA or metaGGA families), the molecular energy reads: hybrid Emol ¼

X

Ts ðAÞ þ

A

X

cl EIQA ðAÞ þ

A

1 X X cl E ðA; BÞ 2 A BaA IQA

ðð 3 ~ r þ a eXX ð~ þ eloc ð r Þd r;~ r0 Þd3 r d3 r0 : xc ð

(A1) The simplest exact energy partition would be (note that EXX x (A, B) is of course calculated using KS orbitals and not the HF): X hybrid cl loc Emol ¼ Ts ðAÞ þ EIQA ðAÞ þ Exc ðAÞ þ aExXX ðAÞ A

þ

1 X X cl EIQA ðA; BÞ þ aExXX ðA; BÞ 2 A BaA X

E self;simpl ðAÞ þ

A

1 X X inter;simpl E ðA; BÞ: 2 A BaA (A2)


1 X X inter;AIMAll E ðA; BÞ: 2 A BaA

Acknowledgements We gratefully acknowledge the CRIHAN computational center for providing HPC resources. VT and LJ thank the Centre National de la Recherche Scientifique (CNRS) for a ‘‘Chaire d’Excellence’’ at the University of Rouen and the Labex SynOrg (ANR-11-LABX-0029). MYO acknowledges the CNRS and the ´gion Haute-Normandie’’ for PhD funding. This work was ‘‘Re also supported by FONDECYT through projects No. 1120093 and No. 1130072, and Grant ICM No. 120082. RIR wishes to thank Grant ICM No. 120082.

References 1 E. G. Lewars, The concept of the potential energy surface, Computational Chemistry, Springer, 2nd edn, 2011. ´, J. Phys. Chem. A, 1999, 103, 4398–4403. 2 A. Toro-Labbe ´, J. Chem. Phys., 2004, 121, 3 B. Herrera and A. Toro-Labbe 7096–7102. ń, P. Jaque and A. Toro-Labbe ´, J. Phys. Chem. A, 4 E. Rinco 2006, 110, 9478–9485. ´, S. Gutie ´rrez-Oliva and J. S. Murray, 5 P. Politzer, A. Toro-Labbe Adv. Quantum Chem., 2012, 64, 189–209.


View Article Online


PCCP

´, 6 P. Politzer, J. R. Reimers, J. S. Murray and A. Toro-Labbe J. Phys. Chem. Lett., 2010, 1, 2858–2862. ´, J. Phys. Chem. A, 2008, 112, 7 E. Echegaray and A. Toro-Labbe 11801–11807. ń, E. Echegaray, S. Gutie ´rrez-Oliva, B. Herrera and 8 M. L. Cero ´ A. Toro-Labbe, Sci. China: Chem., 2011, 54, 1982–1988. ´rrez-Oliva, E. Silva and A. Toro9 P. Flores-Morales, S. Gutie ´ Labbe, THEOCHEM, 2010, 943, 121–126. ´, J. Phys. Chem. A, 2007, 111, 10 B. Herrera and A. Toro-Labbe 5921–5926. ˜ ez, F. Lund and 11 S. Giri, E. Echegaray, P. W. Ayers, A. S. Nun ´, J. Phys. Chem. A, 2012, 116, 10015–10026. A. Toro-Labbe ¨hringer-Martinez and A. Toro-Labbe ´, J. Comput. Chem., 12 E. Vo 2010, 31, 2642–2649. ´, 13 J. V. Correa, P. Jaque, P. Geerlings and A. Toro-Labbe Int. J. Quantum Chem., 2012, 112, 2142–2153. ´, J. Chem. Phys., 2009, 14 S. Vogt-Geisse and A. Toro-Labbe 130, 244308. 15 F. Duarte and A. Toro-Labbe, J. Phys. Chem. A, 2011, 115, 3050–3059. 16 P. Geerlings, F. de Proft and W. Langenaecker, Chem. Rev., 2003, 103, 1793–1874. 17 H. Chermette, J. Comput. Chem., 1999, 20, 129–154. 18 R. G. Parr, R. A. Donnelly, M. Levy and W. E. Palke, J. Chem. Phys., 1978, 68, 3801. 19 R. S. Mulliken, J. Chem. Phys., 1934, 2, 782. 20 R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, 1990. 21 P. L. A. Popelier, Atoms in Molecules An Introduction, Prentice Hall, 2000. 22 M. L. Bender, Mechanisms of Homogeneous Catalysis from Protons to Proteins, John Wiley and Sons, New York, USA, 1971. 23 T. Boutis, Proton Transfer in Hydrogen Bonded Systems, Plenum, New York, USA, 1992. 24 T. J. Zielinski and R. A. Poirier, J. Comput. Chem., 1984, 5, 466–470. ´rrez-Oliva, B. Herrera, A. Toro-Labbe ´ and H. Chermette, 25 S. Gutie J. Phys. Chem. A, 2005, 109, 1748–1751. 26 V. Tognetti and C. Adamo, J. Phys. Chem. A, 2009, 113, 14415–14419. 27 X.-C. Wang, J. Nichols, M. Feyereisen, M. Gutowski, J. Boatz, A. D. J. Haymet and J. Simons, J. Phys. Chem., 1991, 95, 10419–10424. 28 A. Luna, B. Amekraz, J. Tortajada, J. P. Morizur, M. Alcamı´, ´ and M. Yan ˜ ez, J. Am. Chem. Soc., 1998, 120, O. Mo 5411–5426. ´, M. Yan ˜ ez and J.-Y. Salpin, Phys. Chem. 29 A. Eizaguirre, O. Mo Chem. Phys., 2011, 13, 18409–18417. ´, Phys. 30 R. Inostroza-Rivera, B. Herrera and A. Toro-Labbe Chem. Chem. Phys., 2014, 16, 14489–14495. 31 J. C. Slater, J. Chem. Phys., 1933, 1, 687. 32 S. K. Ghosh and R. G. Parr, J. Chem. Phys., 1985, 82, 3307. ´s, M. A. Blanco and E. Francisco, J. Chem. Phys., 33 A. M. Penda 2004, 120, 4581.


Paper

´s and E. Francisco, J. Chem. Theory 34 M. A. Blanco, A. M. Penda Comput., 2005, 1, 1096–1109. ´s, M. A. Blanco and E. Francisco, J. Chem. Theory 35 A. M. Penda Comput., 2006, 2, 90–102. ´s, M. A. Blanco and E. Francisco, J. Chem. Phys., 36 A. M. Penda 2006, 125, 184112. ´s, E. Francisco and M. A. Blanco, Faraday 37 A. M. Penda Discuss., 2007, 135, 423–438. 38 F. Zielinski, V. Tognetti and L. Joubert, Chem. Phys. Lett., 2012, 527, 67–72. ¨rling and M. Levy, Phys. Rev. A: At., Mol., Opt. Phys., 39 A. Go 1994, 50, 196. 40 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, ¨ . Farkas, J. B. Foresman, J. V. Ortiz, A. D. Daniels, O J. Cioslowski and D. J. Fox, Gaussian 09, Revision D.01, Gaussian, Inc., Wallingford CT, 2009. 41 K. Fukui, J. Phys. Chem., 1970, 74, 4161–4163. 42 K. Fukui, Acc. Chem. Res., 1981, 14, 363–368. 43 H. P. Hratchian and H. B. Schlegel, J. Chem. Phys., 2004, 120, 9918. 44 T. A. Keith, AIMAll (Version 13.11.04), TK Gristmill Software, 2013. 45 V. Tognetti and L. Joubert, J. Chem. Phys., 2013, 138, 024102. 46 V. Tognetti and L. Joubert, Chem. Phys. Lett., 2013, 579, 122–126. 47 V. Tognetti and L. Joubert, Phys. Chem. Chem. Phys., 2014, 16, 14539–14550. 48 M. Yahia-Ouahmed, V. Tognetti and L. Joubert, Comput. Theor. Chem., 2015, 1053, 254–262. ´, S. Gutie ´rrez-Oliva, B. Herrera, 49 P. Politzer, A. Toro-Labbe P. Jaque, M. Concha and J. S. Murray, J. Chem. Sci., 2005, 117, 467–472. ´, S. Gutie ´rrez-Oliva, J. S. Murray and 50 A. Toro-Labbe P. Politzer, J. Mol. Model., 2009, 15, 707–710. 51 O. O. Brovarets’ and D. M. Hovorun, J. Comput. Chem., 2013, 34, 2577–2590. 52 O. O. Brovarets’, R. O. Zhurakivsky and D. M. Hovorun, J. Comput. Chem., 2014, 35, 451–466. 53 M.-L. Bonnet and V. Tognetti, Chem. Phys. Lett., 2011, 511, 427–433. 54 X. Fradera, M. A. Austen and R. F. W. Bader, J. Phys. Chem. A, 1999, 103, 304–314.


View Article Online


Paper

55 J. Cioslowski and S. T. Mixon, J. Am. Chem. Soc., 1991, 113, 4142–4145. ńgya ń, M. Loos and I. Mayer, J. Phys. Chem., 1994, 98, 56 J. G. A 5244–5248. 57 M. Garcıá-Revilla, P. L. A. Popelier, E. Francisco and ´s, J. Chem. Theory Comput., 2011, 7, 1704–1711. A. M. Penda ¨ ´, J. Phys. Chem. A, 58 E. Vohringer-Martinez and A. Toro-Labbe 2012, 116, 7419–7423. `, Chem. – Eur. J., 59 J. Poater, X. Fradera, M. Duran and M. Sola 2003, 9, 1113–1122. `, Chem. Soc. Rev., 60 F. Feixas, E. Matito, J. Poater and M. Sola DOI: 10.1039/C5CS00066A, in press.


PCCP

61 A. A. Ebrahimi, R. Ghiasi and C. Foroutan-Nejad, THEOCHEM, 2010, 941, 47–52. 62 M. K. Griffiths and P. L. A. Popelier, J. Chem. Inf. Model., 2013, 53, 1714–1725. 63 Z. Zhou and R. G. Parr, J. Am. Chem. Soc., 1989, 111, 7371–7379. ´, Phys. Chem. 64 V. Labet, C. Morell, A. Grand and A. Toro-Labbe Chem. Phys., 2010, 12, 4142–4151. 65 T. A. Keith, R. F. W. Bader and Y. Aray, Int. J. Quantum Chem., 1996, 57, 183–198. ´s, E. Francisco, M. A. Blanco and C. Gatti, Chem. – 66 A. M. Penda Eur. J., 2007, 13, 9362–9371.


Proton transfer reactions of methylglyoxal synthase.

Non-adiabatic couplings and dynamics in proton transfer reactions of Hn (+) systems: Application to H2+H2 (+)→H+H3 (+) collisions.

Assessing the importance of proton transfer reactions in DNA.

Exploring the metric of excited state proton transfer reactions.

Conceptual DFT Descriptors of Amino Acids with Potential Corrosion Inhibition Properties Calculated with the Latest Minnesota Density Functionals.

Proton transfer reactions and hydrogen-bond networks in protein environments.

The proton transfer reactions catalyzed by yeast pyruvate kinase.

Photoelectric conversion based on proton-coupled electron transfer reactions.

A DFT study on proton transfers in hydrolysis reactions of phosphate dianion and sulfate monoanion.

TDDFT investigation of the excited state proton transfer reaction of fisetin chromophore.

Folding of Protein Ions in the Gas Phase after Cation-to-Anion Proton-Transfer Reactions.

TD-DFT assessment of the excited state intramolecular proton transfer in hydroxyphenylbenzimidazole (HBI) dyes.

Local self similar descriptors: comparison and application to gastroenterology images.

Animal reactions to oncoming vehicles: a conceptual review.

Proton radiography and tomography with application to proton therapy.

Uncovering the Geometry of Barrierless Reactions Using Lagrangian Descriptors.

Theoretical Study of Proton Coupled Electron Transfer Reactions: The Effect of Hydrogen Bond Bending Motion.

Reactions of proton-bound dimers.

Branching Ratio between Proton Transfer and Electron Transfer Channels of a Bidirectional Proton-Coupled Electron Transfer.

Conceptual DFT: chemistry from the linear response function.

Photo-induced proton-coupled electron transfer reactions of acridine orange: comprehensive spectral and kinetics analysis.

Competition and interplay of various intermolecular interactions in ultrafast excited-state proton and electron transfer reactions.

Mechanisms for Two-Step Proton Transfer Reactions in the Outward-Facing Form of MATE Transporter.

Operando NMR spectroscopic analysis of proton transfer in heterogeneous photocatalytic reactions.