Performance of the NOF-MP2 method in hydrogen abstraction reactions

TL;DR

This paper evaluates the NOF-MP2 method's effectiveness in modeling hydrogen abstraction reactions, demonstrating its accuracy in predicting dissociation energies comparable to CASPT2, especially for systems with multiconfigurational character.

Contribution

The study provides a comprehensive assessment of NOF-MP2's reliability in describing electron correlation in radical formation and bond cleavage, with validation against experimental and high-level computational data.

Findings

NOF-MP2 achieves quantitative agreement with experimental dissociation energies.

The method performs comparably to CASPT2 in systems with multiconfigurational character.

It reliably describes both dynamic and static correlation in radical reactions.

Abstract

The recently proposed natural orbital functional second-order M{\o}ller-Plesset (NOF-MP2) method is capableof achieving both dynamic and static correlation even for those systems with significant multiconfigurational character. We test its reliability to describe the electron correlation in radical formation reactions, namely, in the homolytic X-H bond cleavage of LiH, BH, CH4, NH3, H2O and HF molecules. Our results are compared with CASSCF and CASPT2 wavefunction calculations and the experimental data. For a dataset of 20 organic molecules, the thermodynamics of C-H homolytic bond cleavage, in which the C-H bond is broken in the presence of different chemical environments, is presented. The radical stabilization energies obtained for such general dataset are compared with the experimental data. It is observed that NOF-MP2 is able to give a quantitative agreement for dissociation energies, with a performance comparable to that of the accurate CASPT2 method.

Tables

Table 1. Table 1: Dissociation energies, in kcal/mol, calculated from single-point energies a . ZPVEs were added b to the experimental dissociation energies (Ervin and DeTuri, 2002 ) . PNOF6, CASSCF(2,2) and CASPT2(2,2) results are taken from Ref. (Lopez et al., 2015 ) .

	$\mathrm{LiH}\to {}_{}{}^{.}\mathrm{Li}+H^{.}$	$\mathrm{BH}\to {}_{}{}^{.}B+H^{.}$	${\mathrm{CH}}_4{\to }^{.}{\mathrm{CH}}_3+H^{.}$	${\mathrm{NH}}_3{\to }^{.}{\mathrm{NH}}_2+H^{.}$	$H_{2}O\to {}_{}{}^{.}\mathrm{OH}+H^{.}$	$\mathrm{FH}\to {}_{}{}^{.}F+H^{.}$	MAE
	cc-pVDZ
PNOF6c	42.8	78.4	104.6	102.4	106.7	113.2	14.6
PNOF6(3)d	48.5	89.4	112.9	111.2	114.7	119.8	9.1
NOF-MP2c	44.1	57.5	106.9	108.0	117.3	131.4	11.7
NOF-MP2(3)d	49.5	56.7	106.8	105.3	117.1	130.6	11.6
CASSCF(2,2)	42.8	78.2	97.3	95.1	99.4	107.8	19.2
CASPT2(2,2)	49.2	78.7	106.6	106.8	114.8	126.6	8.8
	cc-pVTZ
PNOF6c	44.1	80.7	105.2	104.8	110.4	119.5	11.8
PNOF6(3)d	51.3	92.6	113.3	114.4	120.1	127.8	6.5
NOF-MP2c	55.0	62.7	107.1	113.3	125.4	143.0	5.5
NOF-MP2(3)d	36.2	61.7	105.0	111.0	122.8	141.1	9.6
CASSCF(2,2)	44.0	81.1	98.0	97.3	102.7	113.4	16.5
CASPT2(2,2)	53.4	81.7	109.6	112.2	122.0	136.9	3.4
Exp.	58.0	81.5	113.0	115.9	126.0	141.1

Table 2. Table 2: C-H Bond Dissociation energies, in kcal/mol, for a dataset of 20 organic molecules. ZPVEs at the M062X/cc-pVTZ level of theory were added to the experimental dissociation energies (Henry et al., 2001 ; Blanksby and Ellison, 2003 ) . In case of CH 4 , this leads to a experimental D e of 112.7 kcal/mol, 0.3 kcal/mol lower than the value estimated in Table 1 . PNOF6(3), CASSCF(2,2) and CASPT2(2,2) data is taken from Ref. (Lopez et al., 2015 ) . Calculations carried out with the cc-pVDZ basis set, and considering the X-H distance of 5 Å as the dissociation limit.

		PNOF6(3)	CASSCF(2,2)	CASPT2(2,2)	NOF-MP2(3)	Exp
${\mathrm{CH}}_4$		112.8	97.3	106.6	104.5	112.7
${\mathrm{CH}}_3{\mathrm{CH}}_3$		111.2	96.0	104.2	103.4	109.7
${\mathrm{CH}}_3{\mathrm{CH}}_2{\mathrm{CH}}_3$a		110.8	94.5	102.5	103.9	106.9
${\mathrm{CH}}_3{\mathrm{CH}}_2{\mathrm{CH}}_3$b		112.7	96.5	104.8	105.8	108.5
${\mathrm{CH}}_{3}F$		111.4	97.9	103.0	108.5	108.7
${\mathrm{CF}}_2{H}_2$		118.3	100.3	105.4	116.1	111.8
${\mathrm{CF}}_{3}H$		121.2	101.7	107.3	117.6	113.5
${\mathrm{CH}}_3\mathrm{OH}$		113.4	94.5	99.1	106.7	103.2
${\mathrm{CH}}_3\mathrm{COH}$		110.6	91.7	97.4	101.3	100.8
$H_2\mathrm{CO}$		112.0	81.6	87.2	101.6	95.6
${\mathrm{CH}}_3{\mathrm{OCH}}_3$		113.7	95.6	100.7	108.0	102.5
${\mathrm{CH}}_3\mathrm{COOH}$		113.8	92.9	100.6	102.5	100.8
$C_2{H}_4$		128.9	103.7	115.8	116.9	119.3
$C_2{H}_2$		147.9	127.0	149.7	143.4	141.8
${\mathrm{CH}}_3\mathrm{CCH}$		112.1	92.1	94.5	101.7	97.1
$\mathrm{HCN}$		141.2	126.5	125.8	127.2	132.5
${\mathrm{CH}}_3\mathrm{CN}$		119.4	94.8	98.3	104.2	103.9
${\mathrm{CH}}_3{\mathrm{NO}}_2$		127.5	98.2	103.3	108.5	106.5
${\mathrm{CF}}_3{\mathrm{CH}}_3$		114.6	99.1	107.9	108.7	113.8
$C_6{H}_6$		129.8	101.6	114.6	124.3	120.5
$C_6{H}_6{\mathrm{CH}}_3$		111.9	89.6	103.7	103.4	96.1
MAE		9.0	11.1	5.0	3.7

Equations

$$E={\tilde{E}}_{hf}+E^{corr}={\tilde{E}}_{hf}+E^{dyn}+E^{sta}$$

$$\mathrm{XH}+{\mathrm{Y}}^{.}\to {\mathrm{X}}^{.}+\mathrm{YH}$$

$$\sum _{p\in {\Omega }_g}{n}_p=1,g=1,2,\dots ,\mathrm{N}/2$$

$$\begin{array}[]{c}{\displaystyle 2\sum\limits_{p\in\Omega}n_{p}=2\sum_{g=1}^{\mathrm{N}/2}}{\displaystyle\sum_{p\in\Omega_{g}}}n_{p}=\mathrm{N}\end{array}$$

$$\begin{array}[]{c}E=\sum\limits_{g=1}^{\mathrm{N}/2}E_{g}+\sum\limits_{f\neq g}^{\mathrm{N}/2}E_{fg}\\ \\ E_{g}=\sum\limits_{p\in\Omega_{g}}n_{p}\left(2\mathcal{H}_{pp}+\mathcal{J}_{pp}\right)+\sum\limits_{p,q\in\Omega_{g},p\neq q}\Pi_{qp}^{g}\mathcal{L}_{pq}\\ \\ E_{fg}=\sum\limits_{q\in\Omega_{f}}\sum\limits_{p\in\Omega_{g}}\left[n_{q}n_{p}\left(2\mathcal{J}_{pq}-\mathcal{K}_{pq}\right)+\Pi_{qp}^{s}\mathcal{L}_{pq}\right]\end{array}$$

$$\begin{array}[]{c}\Pi_{qp}^{g}=\left\{\begin{array}[]{cc}-\sqrt{n_{q}n_{p}}\,,&p=g\textrm{ or }q=g\\ +\sqrt{n_{q}n_{p}}\,,&p,q>\mathrm{N}/2\end{array}\right.\\ \\ \Pi_{qp}^{s}=-4n_{q}\left(1-n_{q}\right)n_{p}\left(1-n_{p}\right)\end{array}$$

$$\begin{array}[]{c}E^{sta}=\sum\limits_{g=1}^{\mathrm{N}/2}\sum\limits_{q\neq p}\sqrt{\Lambda_{q}\Lambda_{p}}\,\Pi_{qp}^{g}\mathcal{\,L}_{pq}+\sum\limits_{f\neq g}^{\mathrm{N}/2}\sum\limits_{p\in\Omega_{f}}\sum\limits_{q\in\Omega_{g}}\Pi_{qp}^{s}\mathcal{\,L}_{pq}\end{array}$$

$$E^{dyn}=\sum _{g,f=1}^{\mathrm{N}/2}\sum _{p,q>N/2}^M⟨gf|pq⟩[2T_{pq}^{gf}-T_{pq}^{fg}]$$

$$\begin{array}[]{c}C_{p}^{tra}=\begin{cases}\begin{array}[]{c}\begin{array}[]{c}1-4h_{p}^{2}\end{array}\\ 1-4n_{p}^{2}\end{array}&\begin{array}[]{c}p\leq\mathrm{N}/2\\ p>\mathrm{N}/2\end{array}\end{cases}\\ \>C_{p}^{ter}=\begin{cases}\begin{array}[]{c}\begin{array}[]{c}1\end{array}\\ 1-4h_{p}n_{p}\end{array}&\begin{array}[]{c}p\leq\mathrm{N}/2\\ p>\mathrm{N}/2\end{array}\end{cases}\end{array}$$

$${\tilde{\mathcal{F}}}_{pq}=\{\begin{array}{ll}{C}_p^{tra}{C}_q^{tra}{\mathcal{F}}_{pq},& p,q\in {\Omega }_g\\ C_p^{ter}{C}_q^{ter}{\mathcal{F}}_{pq},& otherwise\end{array}$$

$$\widetilde{⟨pq|rt⟩}=\{\begin{array}{ll}{C}_p^{tra}{C}_q^{tra}{C}_r^{tra}{C}_t^{tra}⟨pq|rt⟩,& p,q,r,t\in {\Omega }_g\\ C_p^{ter}{C}_q^{ter}{C}_r^{ter}{C}_t^{ter}⟨pq|rt⟩,& otherwise\end{array}$$

$$\widetilde{⟨ab|ij⟩}+({\mathcal{F}}_{aa}+{\mathcal{F}}_{bb}-{\mathcal{F}}_{ii}-{\mathcal{F}}_{jj})T_{ab}^{ij}+$$

$$\sum _{c\ne a}{\tilde{\mathcal{F}}}_{ac}{T}_{cb}^{ij}+\sum _{c\ne b}{T}_{ac}^{ij}{\tilde{\mathcal{F}}}_{cb}-\sum _{k\ne i}{\tilde{\mathcal{F}}}_{ik}{T}_{ab}^{kj}-\sum _{k\ne j}{T}_{ab}^{ik}{\tilde{\mathcal{F}}}_{kj}=0$$

$$\mathrm{XH}\to {}^{.}\mathrm{X}+{\mathrm{H}}^{.}$$

$$\mathrm{RS}{\mathrm{E}}^{\mathrm{XY}}={\mathrm{D}}_{\mathrm{e}}^{\mathrm{XH}}-{\mathrm{D}}_{\mathrm{e}}^{\mathrm{YH}}$$

Full text

I Introduction

Natural orbital functional (NOF) theory (Piris, 2007) is being configured as an alternative formalism to both DFT and wavefunction methods, by describing the electronic structure in terms of the natural orbitals (NOs) and their occupation numbers (ONs). Various functionals have been developed in the last years, a comprehensive review can be found in Refs. (Piris and Ugalde, 2014; Pernal and Giesbertz, 2016). Recently (Piris, 2017), a single-reference global method for the electron correlation was introduced taking as reference the Slater determinant formed with the NOs of an approximate NOF. In this approach, called natural orbital functional - second-order Møller–Plesset (NOF-MP2) method, the total energy of an N-electron system can be attained by the expression

[TABLE]

where $\tilde{E}_{hf}$ is the Hartree-Fock energy obtained with the NOs, the dynamic energy ($E^{dyn}$) is derived from a modified MP2 perturbation theory, while the non-dynamic energy ($E^{sta}$) is obtained from the static component of the employed NOF.

In fact, NOF theory is a particular case of the one-particle reduced density matrix (1RDM) functional theory (Gilbert, 1975; Levy, 1979; Valone, 1980), in which the spectral decomposition of the 1RDM is assumed. In this representation, restrictions on the ONs to the range $\left[0,1\right]$ represent the necessary and sufficient conditions for ensemble N-representability of the 1RDM (Coleman, 1963) under the Lowdin’s normalization. The exact functional in terms of the 1RDM has been an unattainable goal so far, and we really work with approximations. Approximating the energy functional implies that the functional N-representability problem arises (Piris, 2018a). To date, only NOFs proposed by Piris and coworkers (Piris, M, 2013) rely on the reconstruction of the two-particle reduced density matrix (2RDM) subject to ensemble N-representability conditions.

The success of the NOF-MP2 method is determined by the NOs used to generate the reference. The functional PNOF7s proved (M. Piris, 2018) to be the functional of choice for the method. The "s" emphasizes that this interacting-pair model takes into account only the static correlation between pairs, and therefore avoids double counting in the regions where the dynamic correlation predominates, already in the NOF optimization. Moreover, the correction $E^{dyn}$ is based on the orbital-invariant formulation of the MP2 energy (Saebo, 2002).

In the present paper, we analyze the performance of NOF-MP2 in the description of X-H bond dissociations, important process in biological (Valko, M.; Rhodes C. J.; Moncol J.; Izakovic, M.; Mazur, M., 2006; Valko, M.; Leibfritz, D.; Moncol, J.; Cronin, M.T.D.; Mazur, M.; Telser, J., 2007) and organic chemistry (Breher, F., 2007). Firstly, we evaluate the dissociation energy for the X-H bonds in LiH, BH, CH4, NH3, H2O and HF molecules. Results are compared to our previous calculations (Lopez et al., 2015), and the experimental data.

The proper description of the X-H homolytic bond dissociation curves is a fundamental step for the accurate characterization of the electronic structure of these important species (Basch and Hoz, 1997; Coote, 2004; Temelso et al., 2006; Vandeputte et al., 2007). This requires the appropriate treatment of strong correlation effects since a single Slater determinant wavefunction leads to incorrect results. We need to include several determinants that lead to computationally demanding methods. An alternative is the density functional theory (DFT), however, it suffers from methodological problems to treat strong electron correlation or near-degeneracy effects (Shao et al., 2003; Krylov, 2006; Ess and Cook, 2001). It is worth noting that cost-effective bond dissociation energies can be obtained in the context of spin-dependent DFT, but at the price of obtaining solutions with breaking symmetry (Menon et al., 2007; Menon and Radom, 2008). Valence bond theory has also been used for this type of systems (Lai et al., 2012).

The formation of radicals by hydrogen abstraction is a fundamental step to explain the oxidation of hydrocarbons (Carstensen et al., 2007; Huynh et al., 2008, 2010), lipid-peroxidation (Valko, M.; Leibfritz, D.; Moncol, J.; Cronin, M.T.D.; Mazur, M.; Telser, J., 2007), formation of reactive oxygen species (Stadtman and Levine, 2003), Fenton chemistry (Prousek, 2007) and DNA damage (Balasubramanian et al., 1998). Due to this widespread interest on the thermodynamic stability of organic radicals, we analyze secondly the cleavage of the C-H bond in a dataset of 20 organic molecules, previously designed in our group (Lopez et al., 2015). As a measure of radical stability we employ the bond dissociation energy ($\mathrm{D_{e}}$) which has often been used in the literature (Vereecken and Peeters, 2001; Blanksby and Ellison, 2003). Based on $\mathrm{D_{e}}$, we estimate the radical stabilization energy (RSE) for a variety of hydrogen abstraction reactions of the type:

[TABLE]

RSE is equivalent to the difference in bond dissociation energies of XC-H and Y-H species.

II Theory

In this work, we address only singlet states, so we adopt the spin-restricted theory in which a single set of orbitals is used for $\alpha$ and $\beta$ spins. We shall use PNOF7s (M. Piris, 2018), which is a NOF based on the electron-pairing approach in NOF theory (Piris, 2018b).

Consider the orbital space $\Omega$ is divided into N/2 mutually disjoint subspaces $\Omega{}_{g}$, so each orbital belongs only to one subspace. Each subspace contains one orbital $g$ below the level N/2, and $\mathrm{N}_{g}$ orbitals above it, which is reflected in additional sum rules for the ONs,

[TABLE]

Taking into account the spin, each subspace contains only an electron pair. The Lowdin’s normalization condition is automatically fulfilled,

[TABLE]

Coupling each orbital $g$ below the N/2 level with only one orbital above it ($\mathrm{N}_{g}=1$) leads to the orbital perfect-pairing approach. In general, we fix $\mathrm{N_{g}}$ to the maximum allowed value determined by the basis set used in calculations. It is important to note that orbitals satisfying the pairing conditions (3) are not required to remain fixed throughout the orbital optimization process (Piris and Ugalde, 2009).

The energy of PNOF7s can be conveniently written as

[TABLE]

where

[TABLE]

$\mathcal{J}_{pq}$, $\mathcal{K}_{pq}$, and $\mathcal{L}_{pq}$ are the usual direct, exchange, and exchange-time-inversion two-electron integrals. The first term of the energy in Eq. (5) draws the system as independent N/2 electron pairs, whereas the second term contains the interactions between electrons belonging to different pairs. PNOF7s provides the reference NOs to form $\tilde{E}_{hf}$ in the NOF-MP2 method, Eq. (1).

$E^{sta}$ is the sum of the static intrapair and interpair electron correlation energies:

[TABLE]

where $\Lambda_{p}=1-\left|1-2n_{p}\right|$ is the amount of intra-pair static correlation in each orbital as a function of its occupancy.

$E^{dyn}$ is obtained from the second-order correction $E^{\left(2\right)}$ of the MP2 method. The first-order wavefunction is a linear combination of all doubly excited configurations, and their amplitudes $T_{pq}^{fg}$ are obtained by solving the equations for the MP2 residuals (Saebo, 2002). The dynamic energy correction takes the form

[TABLE]

where $M$ is the number of basis functions, and $\left\langle gf\right|\left.pq\right\rangle$ are the matrix elements of the two-particle interaction.

In fact, $E^{dyn}$ is the modified $E^{\left(2\right)}$ in order to avoid double counting of the electron correlation. It is divided into intra- and inter-pair contributions, and the amount of dynamic correlation in each orbital $p$ is defined by functions $C_{p}$ of its occupancy, namely,

[TABLE]

According to Eq.(9), fully occupied and empty orbitals yield a maximal contribution to dynamic correlation, whereas orbitals with half occupancies contribute nothing. Using these functions as the case may be (intra-pair or inter-pair), the modified off-diagonal elements of the Fock matrix ($\tilde{\mathcal{F}}$) are defined as

[TABLE]

as well as modified two-electron integrals:

[TABLE]

where the subspace index $g=1,...,\mathrm{N}/2$. This leads to the following linear equation for the modified MP2 residuals:

[TABLE]

[TABLE]

where $i,j,k$ refer to the strong occupied NOs, and $a,b,c$ to weak occupied ones. It should be noted that diagonal elements of the Fock matrix ($\mathcal{F}$) are not modified. By solving this linear system of equations the amplitudes $T_{pq}^{fg}$ are obtained, which are inserted into the Eq. (8) to achieve $E^{dyn}$.

All calculations have been carried out using the DoNOF code developed by M. Piris and coworkers. The procedure is simple, showing a formal scaling of $M^{5}$ ($M$: number of basis functions). However, our implementation in the molecular basis set requires also four-index transformation of the electron repulsion integrals, which is a time-consuming step, though a parallel implementation of this part of the code has substantially improved its performance. As a result, the possibility of addressing large systems opens up.

III Results and Discussion

Results are organized as follows. First, the X-H bond dissociation energies for LiH, BH, CH4, NH3, H2O and HF molecules are studied using NOF-MP2. Next, we analyze the performance of NOF-MP2 for describing C-H bond cleavage in a variety of 20 organic molecules. Finally, radical stabilization energies are calculated based on the calculated C-H bond dissociation energies. In all calculations, recall that the maximum value allowed by the basis set used is assumed for $\mathrm{N}_{g}$ by default. In NOF-MP2(3) calculations, only three orbitals ($\mathrm{N}_{g}=3$) above the N/2 level in each electron pair are considered.

Geometries are taken from our previous publication (Lopez et al., 2015), which were obtained at the M06-2X level of theory (Zhao and Truhlar, 2008). The dissociation limit is calculated by considering a frozen X-H distance of 5Å, and optimizing the rest of internal coordinates. At these geometries, single-point energies are evaluated at the NOF-MP2 level of theory. The correlation-consistent valence double-$\zeta$ (cc-pVDZ) or triple-$\zeta$ (cc-pVTZ) basis sets developed by Dunning et al. (Dunning and Dunning Jr., 1989) are used. The zero point vibrational energies (ZPVEs) were taken from the NITS Computational Chemistry Comparison and Benchmark Database (CCCBDB) (Johnson III, 2018), and corresponds to CCSD(T)/cc-pVTZ values.

We also provide the PNOF6 (Piris, 2014) and wavefunction-based calculations obtained in Ref. (Lopez et al., 2015). For the latter, an active space was defined by the distribution of two electrons in two molecular orbitals, CASSCF(2,2) (Roos et al., 1980; Siegbahn et al., 1980). The dynamic correlation effects were included through complete active space second-order perturbation theory calculations, CASPT2(2,2) (Andersson et al., 1992). MOLCAS 7.0 suite of programs (Aquilante et al., 2009) was used in Ref. (Lopez et al., 2015), for these wavefunction-based calculations.

III.1 X-H homolytic bond cleavage

X-H bond dissociation energies were calculated according to the following reaction:

[TABLE]

with X = Li, B, CH3, NH2, OH, F. The results are presented in Figure 1 and Table 1. The different hydrides considered expand a wide range of dissociation energies, from 58.0 kcal/mol for LiH to 141.1 kcal/mol for FH. The ordering in dissociation energies is $\mathrm{LiH}<\mathrm{BH}<\mathrm{CH_{4}}<\mathrm{NH_{3}}<\mathrm{H_{2}O}<\mathrm{FH}$. In general, NOF-MP2 reproduces satisfactorily these trends.

Let us focus our attention, for example, on the delicate case of the CH4/NH3 ordering. The difference in experimental dissociation energies for these two molecules is very small, only 2.9 kcal/mol with $\mathrm{NH_{3}}$ having a higher dissociation energy. NOF-MP2 is able to reproduce the correct ordering $\mathrm{CH_{4}<NH_{3}}$, except for NOF-MP2(3)/cc-pVDZ. It should be noted that CASSCF(2,2) and PNOF6 gives the reverse order, whereas CASPT2(2,2) recovers the right trend.

It is well known that to reach the experimental values we must go to the complete basis set limit. Therefore, taking into account the moderate basis sets used here, we can say that a good semi-quantitative agreement has been achieved with the experimental data by the NOF-MP2 method. In general, NOF-MP2 shows an intermediate performance between the CASSCF(2,2) and CASPT2(2,2) methods, and a significant improvement with respect to the previously tested PNOF6.

For the six reactions considered, a mean absolute error (MAE) of 5.5 kcal/mol is obtained at NOF-MP2/cc-pVTZ level of theory. NOF-MP2(3)/cc-pVTZ, leads to a higher MAE, namely 9.6 kcal/mol, but this is mainly due to LiH case. For the latter, only one effective pair appears so more $\mathrm{N}_{g}$ orbitals are needed in order to describe properly the dominant intra-pair electron correlation (Piris, M, 2013) in this system. On the other hand, for $\mathrm{CH_{4},NH_{3},H_{2}O}$ and FH, NOF-MP2(3) yield very reasonable results. Therefore, we can say that $\mathrm{N}_{g}=3$ is a good compromise for the characterization of the electron pairs, except for small systems like LiH and BH.

Comparing the performance of NOF-MP2 with wavefunction methods, it is clear that NOF-MP2 and NOF-MP2(3) show a better performance than CASSCF(2,2) (MAE=16.5 kcal/mol with the cc-pVTZ basis set). Introduction of dynamical electron correlation at the CASPT2(2,2) level of theory, reduces the MAE to 3.4 kcal/mol, however, if we reduce the set to CH4, NH3, H2O and FH molecules, there is a similar performance of NOF-MP2 with respect to CASPT2(2,2) method.

III.2 Hydrogen Abstraction in a Dataset of 20 Organic Molecules

III.2.1 Dissociation Energies

As in our previous work (Lopez et al., 2015), we have considered a dataset of 20 organic molecules to evaluate the performance of NOF-MP2 for the C-H bond dissociation energy ($\mathrm{D_{e}^{CH}}$). The selected set covers a wide range of $\mathrm{D_{e}^{CH}}$ values, from 95.6 kcal/mol (H2CO) to 141.8 kcal/mol (C2H2), showing the sensitivity of the C-H bond to different chemical environments. We have considered functional groups with different degree of electron withdrawing/donating ability (-F, -OH, -NO2, -CN, -CH3, …), aromaticity (-C6H5), variety of C-X bonds (HCN, H2CO, CH3NO2, CH3CF3, …), different chain lengths (CH4, CH3CH3, CH3CH2CH3) and different C-C bond orders, single (as in CH3CH3), double (as in C2H4) and triple (as in C2H2). We have decided to use the cc-pVDZ basis set due to the large number of compounds to be treated.

The results can be found in Table 2 and Figure 2. The agreement between NOF-MP2(3) and experimental values is remarkable, with a MAE of 3.7 kcal/mol, even smaller than the MAE for the very accurate CASPT2(2,2) method, namely 5.0 kcal/mol. Notice that previously tested PNOF6(3) method has a MAE of 9.0 kcal/mol, slightly better than CASSCF(2,2), 11.1 kcal/mol. Thus, NOF-MP2(3) method allows for a quantitative description of these dissociation energies, with a similar degree of accuracy as CASPT2(2,2).

Specifically, NOF-MP2(3) is able to reproduce important trends in C-H bond energies. For instance, the experimental $\mathrm{D_{e}^{CH}}$ increases in the following order (Blanksby and Ellison, 2003): CH3CH3 (109.7) < C2H4 (119.3) < C2H2 (141.8) . NOF-MP2(3) is able to reproduce properly this trend, namely, CH3CH3 (103.4) < C2H4 (116.9) < C2H2 (143.4).

The effect of aromaticity can be inferred from the comparison of these dissociation energies with that of the phenyl C-H bond. C6H6, with a formal 1.5 C-C bond order, shows a high dissociation energy (120.5 kcal/mol) even slightly larger than that observed (119.3 kcal/mol) in C2H4, with a formal bond order of 2. This is a clear signature of aromaticity in C6H6, partially lost upon hydrogen abstraction and radical formation. NOF-MP2(3) yields larger values of bond dissociation energies for benzene than for ethene, with values of 124.3 kcal/mol and 116.9 kcal/mol, respectively.

In the case of the benzylic C-H bond (C6H6CH2-H), the effect of the aromaticity works in the opposite direction. In this case, the C-H cleavage does not break the aromaticity, furthermore, the radical itself is stabilized by the aromatic character of the phenyl ring, and consequently, one obtains a much lower $\mathrm{D_{e}^{CH}}$ than for C6H6, namely 96.1 kcal/mol versus 120.5 kcal/mol. NOF-MP2(3) correctly describes this effect, $\mathrm{D_{e}^{CH}}$ for the benzylic C-H bond (103.4 kcal/mol) is also much lower than for the phenyl C-H bond (124.3 kcal/mol) at a magnitude very similar to the experimental value. It is remarkable the right description of aromatic radical stabilization by the NOF-MP2 method, since aromatic stabilization is key to describe radical stability in chemistry.

The chain length is also a factor influencing the C-H bond strength (Vereecken and Peeters, 2001; Carstensen et al., 2007; Huynh et al., 2010). It is known that a larger chain stabilizes the resulting radical: observe the first 3 lines of Table 2. However, NOF-MP2(3) exhibits a poorer sensitivity of radical stability towards chain-lengths with a similar $\mathrm{D_{e}^{CH}}$ for these three molecules. On the other hand, if we consider the same alkane, CH3CH2CH3, and measure both possibilities for hydrogen abstraction, namely, from the central -CH2- or from the terminal -CH3 group, NOF-MP2(3) correctly reproduces the more favorable hydrogen abstraction from the central carbon by 1.9 kcal/mol.

There is also a sizable effect in hydrogen abstraction upon the inclusion of electron withdrawing groups. For instance, fluorination (Korchowiec, 2002) and oxidation (Blanksby and Ellison, 2003) of methane tend to alter the dissociation energy of the C-H bond. Regarding fluorination, a decrease of $\mathrm{D_{e}^{CH}}$ is observed upon the inclusion of a first flour, from 112.7 kcal/mol (CH4) to 108.7 kcal/mol in (CH3F). However, upon higher degree of fluorination in the fluoromethane, $\mathrm{D_{e}^{CH}}$ increases again, 111.8 kcal/mol in CF2H2 and 113.5 kcal/mol in CF3H. NOF-MP2(3) yields a higher $\mathrm{D_{e}^{CH}}$ for CH3F (108.5 kcal/mol) than for CH4 (104.5 kcal/mol). Nevertheless, NOF-MP2(3) describes the proper trend in increasing $\mathrm{D_{e}^{CH}}$ with the degree of fluorination in fluoromethane, namely, CH3F (108.5 kcal/mol) < CH2F2 (116.1 kcal/mol) < CHF3 (117.6 kcal/mol).

With respect to the oxidation of a methyl group, NOF-MP2(3) gives the right trend. For instance, in going from CH3OH to H2CO, there is an important reduction in C-H bond strength, from 103.2 kcal/mol to 95.6 kcal/mol. NOF-MP2(3) yields a similar, although more discrete reduction, of $\mathrm{D_{e}^{CH}}$ from 106.7 kcal/mol to 101.6 kcal/mol.

In general, we can conclude that NOF-MP2(3) represents an accurate balance between dynamical and non-dynamical electron correlation for this set of molecules, yielding $\mathrm{D_{e}^{CH}}$ values that are of the CASPT2(2,2) quality.

III.2.2 Radical Stabilization Energies

RSEs are defined as the energy change in the isodesmic reaction for hydrogen abstraction (Wood et al., 2005; Menon et al., 2007; Menon and Radom, 2008) of Eq. (2). Thus, the RSE for a pair X,Y is defined as

[TABLE]

For the dataset of 21 dissociation energies of Table 2, there are 210 possible combinations of RSEs. It provides with an extensive dataset for the determination of the suitability of a given method to estimate the effect of the substituents on the radical stability in organic molecules. The results for NOF-MP2(3) are summarized in Fig. 3, compared to the performance of wavefunction methods such as CASSCF(2,2) and CASPT2(2,2). In general, there is a reasonable agreement with experimental RSEs for NOF-MP2(3) with a MAE of 5.5 kcal/mol. Slightly better values are obtained for CASSCF(2,2) (4.4 kcal/mol) and CASPT2(2,2) (4.1 kcal/mol) levels of theory.

Another way to compare the results with respect to experimental values is to calculate the linear fit of the theoretical versus the experimental values, and determine the correlation coefficient ($r$). In this sense, NOF-MP2(3) shows a similar performance to the CASSCF(2,2) and CAS2PT2(2,2) methods with an $r$ of 0.9314 versus a value of 0.9482 for both wavefunction methods. In summary, taking into account the large number of hydrogen abstraction reactions considered, the correlation between NOF-MP2(3) and experimental data is highly satisfactory, yielding a quantitative agreement with respect to well established wavefunction methods such as CASPT2(2,2), and providing results close to chemical accuracy.

IV Conclusions

The recently proposed parameter-free natural orbital functional second-order Møller–Plesset (NOF-MP2) method has been applied to the description of radical formation reactions, a delicate problem in quantum chemistry. The application of NOF-MP2(3) to the calculation of the C-H bond dissociation energy in a dataset of 20 organic molecules, and the estimation of the corresponding radical stabilization energies support the use of NOF-MP2(3) as a quantitative theory for the description of these important set of reactions. Comparison of NOF-MP2 with experimental data reveals a similar performance of NOF-MP2 to well-established wavefunction methods such as CASPT2 for these type of problems. We conclude that NOF-MP2 is capable of recovering both dynamical and non-dynamical electron correlation effects in this type of systems. NOF-MP2 is a global electron correlation method for the description of radical stability, which provides results close to chemical accuracy as the widely used and well-established CASPT2 wavefunction method.

Acknowledgements.

Financial support comes from Ministerio de Economía y Competitividad (Ref. CTQ2015-67608-P). The authors thank for technical and human support provided by IZO-SGI SGIker of UPV/EHU and European funding (ERDF and ESF).