Effect of Sintering Temperature on the Photocatalytic Activity of Flower-Like Bi₂WO₆

1.   Introduction

The concepts of electrophilicity and nucleophilicity are cornerstones of our understanding of chemical reactivity ^{1, 2}. Usually, they are merely considered in a qualitative way, together with the notions of electron accepting/donating groups, to rationalize chemical reactivity trends. However, their remarkable utility has also inspired researchers to find ways to quantify the electrophilic and nucleophilic character of a given molecule. These efforts have deepened our understanding of what makes a compound a better or worse electrophile/nucleophile, how these properties could be tuned, and how they relate to other quantities. We view the work of Mayr et al. ^3-7 as especially important. These authors have established electrophilicity and nucleophilicity scales by means of extensive experimental tests. The defining element of these studies is the following expression:

where k_{20 ℃} is the rate constant of the reaction (at 20 ℃) between the nucleophile and the electrophile, while s and N_Mayr are parameters specific to the former, and E_Mayr to the latter. Analyzing several combinations of electrophile-nucleophile pairs, the corresponding parameters can be determined. These definitions of the nucleophilicity and electrophilicity is clearly based on the concept that strong electrophiles react with strong nucleophiles rapidly and vigorously, while weak electrophiles and weak nucleophiles react slowly and reluctantly, if at all.

Parallel to these experimental determinations, there is steadfast interest in finding mathematical equations to quantify the electrophilicity of a molecule ^8-10. Perhaps the best example of this is the early work of Maynard et al. ¹¹, which proposed a simple expression for the electrophilicity in terms of electron attachment/removal energies. Parr et al. ¹² then provided a theoretical justification for Maynard's proposal within the framework of conceptual density functional theory (c-DFT) ^13-23. This landmark contribution gave chemists a simple, yet powerful, tool that has been successfully applied to the analysis of regioselectivity of cycloaddition reactions ^24-27, the study of leaving groups ^{28, 29}, transition states ^30-32, and redox reactions ^33-35. This c-DFT electrophilicity measure has also been used as a descriptor in several quantitative structure activity relationship (QSAR) studies ^36-39. But arguably the biggest impact of the precise quantitative definition of the electrophilicity is that it provided a way to prove the minimum electrophilicity principle (MEP) ^40-44. Expanding on Parr's work, and also within the c-DFT formalism, Gázquez et al. ⁴⁵, and Chamorro et al. ^{46, 47} have generalized this index. Recently, new working expressions for this index has been proposed, which are free from some inconsistencies found in previous versions ^{48, 49}.

The proliferation of theoretical electrophilicity measures makes it desirable to assess their performance. This is main goal of this work. As a reference for comparison, we will use Mayr's parameter. Since the theoretical and experimental measures of electrophilicity refer to different effects, we do not expect them to be exactly the same. However, we believe that a good theoretical measure of electrophilicity must capture the trends from Mayr's work. In addition to testing existing theoretical definitions, we will also attempt to use machine learning techniques to find a mathematical expression that reproduces the Mayr electrophilicity.

2.   Conceptual DFT electrophilicity measures

Parr et al. ¹² defined the electrophilicity, ω, as the stabilization energy associated with the electronic saturation of the system (e.g., molecule). In other words, how the energy of a system changes when it reaches equilibrium with a perfect electron donor. Evaluating this definition requires a model for how the electronic energy, E, changes with respect to the number of electrons, N. If we neglect the changes in the external potential (e.g., molecular geometry) during the electron saturation process, we can estimate energy variation in terms of the following (truncated) Taylor expansion ⁵⁰:

Within the c-DFT framework is customary to identify the coefficients of this expansion with different reactivity descriptors. For example, Parr et al. ⁵¹ proposed to identify the coefficient of the linear term with the (electronic) chemical potential of the system, μ, which is in turn defined as the additive inverse of the electronegativity, χ:

In a similar way, Parr and Pearson ⁵⁰ identified the coefficient of the quadratic term with the chemical hardness, η, which is in turn defined as the multiplicative inverse of the chemical softness, S:

In c-DFT, a perfect electron donor is defined as a system with zero electronegativity (e.g., zero chemical potential). This means that Parr's electrophilicity definition can be written as:

where μ₀ is the chemical potential of the isolated system.

The standard way to approximate this index is to notice that, according to Sanderson's electronegativity equalization principle ⁵², the system will reach electronic saturation when:

Substituting Eq. (2) in this expression we can calculate the maximum number of bound electrons, △N_max:

From this it is now easy to see that:

Notice that we have explicitly indicated that this expression corresponds to the quadratic energy model used in Eq. (2). (It is possible to generalize this expression to differentiate between the electrophilicity of spin up and spin down electrons ^53-60, but this will not be considered here.)

To obtain a working expression for the electrophilicity we need to express the chemical potential and chemical hardness in terms of quantities from electronic structure theory that can be computed using standard quantum chemistry software. This can be done by viewing Eq. (2) as a parabolic interpolation of the electronic energy of the system with N₀, N₀ + 1, and N₀ -1 electrons or, equivalently, using a three-point finite differences approximation with △N = 1 ⁶¹:

here, I and A are the ionization energy and electron affinity, respectively. Substituting these expressions in Eq. (8) we obtain:

Motivated by the success of the parabolic model, Gázquez, Cedillo, and Vela (GCV) proposed a two-parabolas model that allowed them to describe electron transfer processes in a more realistic way ⁴⁵. The key insight in their work was to distinguish between the electron-accepting and the electron-donating processes. (Since we will be only concerned here with the ability of a system to incorporate electron, we will neglect the latter.) In their formulation, the electron-rich and electron-poor regimes are described by different parabolas (e.g., different E vs N models) with the same curvature. Imposing different conditions on this simple model, they obtained the expressions for the chemical potentials of the system when is gaining, μ_GCV⁺ electrons:

Their model is completed with the assumption that the chemical hardness for both electron-accepting and electron-donating processes is equal:

These new expressions can be substituted in the general working equation obtained by Parr (cf. Eq. (8)), resulting in the following electroaccepting, ω_GCV⁺, power⁴⁵:

Inspired by the work of GCV, Chamorro, Duque, and Pérez (CDP) ^{46, 47} rewrote Eq. (8) using the electron accepting, μ₊, chemical potential obtained from the piecewise linear interpolation of Perdew, Parr, Levy, and Balduz ⁶²:

Therefore, they obtained:

Additionally, these authors also proposed a simpler ansatz, where the chemical hardness of a system gaining electrons could be related to the chemical potential of the inverse process. In other words, making ^{46, 47}:

From here results:

Recently, it has been pointed out that there are several problems with the classical way to calculate the electrophilicity, Eq. (11) ⁴⁹. Most notably, this formula gives unreasonable results for polycations (where △N_max > > 1) and for species with small electron affinities (where △N_max should be close to zero, but may not be). These, and other, problems could be traced back to the fact that while the parabolic interpolation is mathematically consistent ⁶³, it lacks a rigorous physical justification ^{64, 65} (except for some cases) ^{66, 67}. This situation motivated the analysis of the electrophilicity index using a finite temperature formulation of the grand-canonical ensemble (an approach that has provided rigorous foundations to several c-DFT results) ^{48, 68-75}. Within this formulation, it is easy to obtain an exact working expression that preserves the physical meaning of Parr's definition. This thermodynamic electrophilicity, ω_TD, is simply given by ⁴⁹:

where ^M0A_k is the kth electron affinity of the system with M₀ particles, and θ(x) is Heaviside (step) function. If we are only considering the same three states used in the parabolic model (e.g., "neutral", "cationic", and "anionic" states), this result is reduced to:

This merely takes the electrophilicity to be equal to A when A > 0, and zero in the other cases. This is the simplified expression that we will be using in the following.

Finally, it has been also shown that a consistent description of chemical reactivity within c-DFT requires the descriptors to incorporate the effects of the molecular surroundings ^75-79. The simplest way to do so is using the insights obtained from the Klopman-Salem frontier molecular orbital treatment ^80-86, which allows us to obtain model expressions for the perturbed chemical potential, μ_P, and the perturbed chemical hardness, η_P ⁷⁶:

where γ and ζ are non-negative parameters that model the interaction with the environment. Including these modified descriptors in Eq. (8) we obtain a more general, perturbed electrophilicity, ω_P:

In the following, we will analyze how the previously discussed electrophilicity measures correlate with Mayr's electrophilicity parameter. To do this, we will work with a simple linear correlation:

Since the effect of the $\frac{1}{{2\xi {{(1 + \gamma )}^2}}}$ factor in Eq. (23) can be included in the slope, m, of the linear regression, we will work with a simplified version of the perturbed electrophilicity, namely:

3.   Electrophilicity measures from symbolic regression (SR)

The indices presented in the previous section are chemically motivated and, as such, they help us determine the factors which govern the acceptance of electrons by an electrophile. While this is a valuable approach, sometimes we are more interested in obtaining an accurate estimate of the electrophilicity, regardless of the insight, it might provide. An attractive alternative in this respect is the use of machine learning (ML) tools ^87-91.

Here, we will work with a very powerful data mining tool: symbolic regression (SR) ^92-100. SR is a supervised method, which means that we provide the algorithm with a "training set" in order to find correlations between different descriptors and a property of interest. The goal of SR is to use a training set to obtain a mathematical function that relates the desired properties to the data. In this way, SR gives the best functional form (of a given pre-specified complexity), along with the best numerical parameters, for a given training set. In more technical terms, given a set of data values $\{ \vec x;y\}$ (training set), with $\vec x$ representing a vector of independent variables (e.g., descriptors), and y the dependent variable (which will be predicted from $\vec x$), SR finds the function $\bar y = f(\vec x)$ that best describes the overall data. The optimization process takes place over a function space, expanded by mathematical operators and constants ¹⁰⁰. If the function space includes all possible functional forms, that is, the length of the mathematical expression can go up to infinity and there is no restriction on the operators used, the method will find many different functions which fit all the data exactly. However, this might cause problems when describing points/observations outside the training set. It is therefore common to enforce some constraints during the SR procedure (e.g., the length of the functions, the number of operands and numerical parameters used, etc.)

Once we restrict the accessible function space, an intuitive choice would be to select a set of basis functions to expand it. However, while this will turn the SR problem into a more familiar linear regression, the size of the basis needed makes this impractical in most cases. This demands the use of different optimization methods, with evolutionary algorithms being a popular choice. Here we will focus on two such algorithms: genetic programming (GP) ^100-103 and grammatical evolution (GE) ^104-106.

3.1   Genetic programming (GP)

Genetic programming (GP) ^100-102 draws inspiration from the mechanism of DNA replication, with new "generations" (e.g., functional forms) leading towards better solutions. On this algorithm, each possible solution is referred to as an "individual". In the case of an SR, these are functions proposed to fit the data. These individuals are commonly represented as a tree (see Fig. 1). The leaves represent constants or data variables, and the nodes represent the allowed mathematical operators. It is important to remark that it is possible if desired, to use basically any function as an operator, and not only the four basic arithmetic operators.

图 1

图 1  Tree representation of Parr's electrophilicity equation.

Figure 1.  Tree representation of Parr's electrophilicity equation.

下载: 全尺寸图片 幻灯片

The algorithm starts with a random generation of possible solutions (trees) that fulfill the restrictions imposed (e.g., the operators allowed and the maximum length of the solution). Then, a new generation of solutions is created, following three steps: breeding, evaluation, and selection. The breeding involves the creation of a new set of alternative solutions. This can be done in several ways ¹⁰⁷: copying a solution to the new generation; combining two solutions to produce two new ones by exchanging two randomly selected subtrees between the parents (crossover) ¹⁰⁸; changing part of a solution (subtree) for a new, randomly generated, one (mutation) ¹⁰³ (see Figs. 2 and 3). The crossover operator ensures that the individuals of the new generation share characteristics of their parents; while the mutation process prevents premature convergence of the algorithm by adding randomness. These changes occur in an uncontrolled way. Then, the evaluation step consists in the evaluation of the fitness function for each of the newly created candidate solutions. The last step consists of the selection of some of the newly created solutions in order to conform the population of the new generation. The probability to select an individual increases with the quality of the fitness calculated in the evaluation step. This process is repeated until a termination criterion is fulfilled (typically given by the number of iterations or a fitness threshold).

图 2

图 2  Effect of the mutation operator on a tree.

Figure 2.  Effect of the mutation operator on a tree.

下载: 全尺寸图片 幻灯片

图 3

图 3  Effect of the crossover operator on two trees.

Figure 3.  Effect of the crossover operator on two trees.

下载: 全尺寸图片 幻灯片

3.2   Grammatical evolution (GE)

Grammatical evolution (GE) ^104-106, a more recent evolutionary algorithm, uses a grammar to map an integer array into a program (in the case of an SR, into a function). To do the mapping, GE uses a context-free grammar ^109-111, which along with an alphabet of terminal symbols, a set of non-terminal symbols, production rules, and a start symbol, is capable of deriving a model from an integer string. These integer strings are then optimized by an evolutionary technique. At every step, the composition rules are used to map from the array to the corresponding model, in order to evaluate the fitness of the latter. This resembles the natural process in which the genotype is separated from the phenotype. The optimization procedure can mimic the steps used in the GP case, but other algorithms such as particle swarm optimization ¹¹² can also be used.

4.   Model systems and computational tools

The different electrophilicity measures were tested over a set of 58 molecules for which Mayr's E_Mayr parameter is known¹¹³ (see Table 1). This global set includes several families of electrophiles, grouped following Mayr's "star classification system" ¹¹⁴. According to this criterion, we have a group of 3 1-star electrophiles, 21 3-stars electrophiles, 22 5-stars electrophiles, and 12 neutral electrophiles.

表 1

表 1  Studied molecules, their classification according to Mayr's system, their experimental EMayr values ¹¹³, and their calculated I and A (given in eV).

Table 1.  Studied molecules, their classification according to Mayr's system, their experimental EMayr values ¹¹³, and their calculated I and A (given in eV).

下载: 导出CSV

Classification	Molecule	EMayr	I	A
1-star	methyl(phenyl) methyleneammonium	-5.17	13.399	5.145
1-star	vilsmeier ion	-5.77	16.255	5.233
1-star	Dimethylmethyl eneammonium	-6.69	17.321	5.273
3-stars	tropylium ion	-3.72	15.627	5.784
3-stars	2-phenyl-1.3-dithiolan-2-ylium	-5.91	13.050	5.718
3-stars	2-phenyl-1.3-dithianylium	-6.43	12.806	5.431
3-stars	1.1 -bis(4-dimethylaminophenyl)-3-phenylallylium	-8.97	9.480	4.749
3-stars	1-(4-chlorophenyl)cyclopent-2-enylium	3.20	12.716	6.145
3-stars	(pfp)2CH+	5.40	12.716	6.356
3-stars	1 -phenylcyclopent-2-enylium	2.89	13.155	6.068
3-stars	N-ethylidenecarbazolium	2.41	11.750	5.264
3-stars	methoxy-(4-methylphenyl) methylium	1.90	13.652	5.649
3-stars	1.3-diphenylpropyn-1 -ylium-Cr(CO)3	1.07	10.641	6.100
3-stars	methoxy-(4-methoxyphenyl) methylium	0.14	12.916	5.306
3-stars	1.1 -dianisyl-3-phenylallylium	-2.67	10.622	5.373
3-stars	methoxy-phenylmethylium	2.97	13.995	5.851
3-stars	flavylium	-3.45	12.094	5.570
3-stars	methoxy-flavylium	-4.95	11.238	5.265
3-stars	acridinium	-7.15	11.886	5.321
3-stars	1.1.3-tris(4-dimethylphenyl) allylium	-9.84	8.901	4.468
3-stars	pop(tol)CH+	2.16	11.156	5.663
3-stars	1.1.3-triphenylallylium	0.98	11.349	5.925
3-stars	1.3-dithianylium	-2.14	14.606	5.674
3-stars	fc(Me)CH+	-2.57	11.699	4.755
5-stars	(ani)2CH+	0.00	11.186	5.429
5-stars	(fur)2CH+	-1.36	10.873	5.307
5-stars	(mfa)2CH+	-3.85	10.224	5.090
5-stars	(mor)2CH+	-5.53	9.719	4.806
5-stars	(dma)2CH+	-7.02	9.809	4.691
5-stars	(pyr)2CH+	-7.69	9.533	4.533
5-stars	(thq)2CH+	-8.22	9.456	4.489
5-stars	(pcp)2CH+	6.02	12.205	6.417
5-stars	(ind)2CH+	-8.76	9.456	4.506
5-stars	(jul)2CH+	-9.45	9.149	4.354
5-stars	(lil)2CH+	-10.04	9.160	4.294
5-stars	benzhydrylium ion	5.90	12.785	6.289
5-stars	pfp(Ph)CH+	5.60	12.750	6.321
5-stars	tol(Ph)CH+	4.59	12.392	6.107
5-stars	(tol)2CH+	3.63	12.038	5.945
5-stars	ani(Ph)CH+	2.11	11.915	5.801
5-stars	ani(tol)CH+	1.48	11.595	5.665
5-stars	ani(pop)CH+	0.61	10.843	5.437
5-stars	(mpa)2CH+	-5.89	9.597	4.643
5-stars	pop(Ph)CH+	2.90	11.493	5.859
5-stars	(dpa)2〇H+	-4.72	8.969	4.752
5-stars	fc(Ph)CH+	-2.64	10.952	5.323
neutral	p-(methoxy)	-10.80	8.443	1.437
neutral	p-(dimethylamino)	-13.30	7.589	1.149
neutral	ani(Br)2QM	-8.63	6.790	1.488
neutral	5-methoxyfuroxano[3.4-d]	-8.37	9.358	2.072
neutral	benzylidenemalononitril	-9.42	9.085	1.626
neutral	benzaldehyde-boron	1.12	9.413	2.127
neutral	ani(Ph)2QM	-12.18	6.595	1.470
neutral	dma(Ph)2QM	-13.39	5.773	1.515
neutral	tol(t-Bu)2QM	-15.83	6.848	1.269
neutral	ani(t-Bu)2QM	-16.11	6.582	1.069
neutral	dma(t-Bu)2QM	-17.29	5.941	1.01 7
neutral	jul(t-Bu)2QM	-17.90	6.438	1.144

The structures of all the molecules were optimized at the HF/6-31g(d) level of theory. Then, single-point calculations were performed for these geometries after adding one electron, and after removing one electron. In this way, we obtain the necessary data to calculate the vertical ionization potentials and electron affinities (see Table 1). All the calculations were performed using Gaussian 09 ¹¹⁵.

The SR calculations were performed using GP and GE algorithms, with different "tree lengths". The mutation probability was set to 15%, and the maximum number of generations was set to 50. For each case, three runs were performed, as a way to check the stability of the SR results. The final expressions were simplified, when possible. These calculations were performed using the HeuristicLab-3.3.12 software package ¹¹⁶.

5.   Results

5.1   Conceptual DFT results

In Table 2 we show the values of the fitted parameters for every tested model. We have performed 5 different fits, one for each of the 4 separate families of compounds, and one that takes into account all the molecules. It is interesting to note that in the case of the 1-star family the slope of the linear regression (cf. Eq. (24)) is negative when we use the ω_parabolic, ω_TD, and ${\tilde \omega _{\rm{P}}}$ indices. This is a surprising result since these are precisely the electrophilicity measures that have the strongest theoretical justification. This implies that these indices can even give qualitatively wrong results, namely, they are in some cases incapable of predicting the relative order of electrophilicity observed experimentally. This is a reminder that these simple electronic structure-based indices miss many of the factors that govern the tendency of a system to gain electrons during the course of an actual chemical reaction. Nonetheless, it is reassuring that when we work with bigger datasets, all the studied models have positive slopes. This means that, overall, these indices are able to recover the rough trends in Mayr's electrophilicity parameter.

表 2

表 2  Fit parameters for all the c-DFT-based electrophilicity measures studied.

Table 2.  Fit parameters for all the c-DFT-based electrophilicity measures studied.

下载: 导出CSV

Electrophilicity Measure	Group	m	b	γ
${{\omega }_{\text{parabolic}}}$	1-star	-16.74420	81.9794	-
	3-stars	3.98655	-25.6185	-
	5-stars	6.39720	-39.5342	-
	neutral	13.71310	-34.4945	-
	all	2.78903	-18.0825	-
$\omega _{{\rm{GCV}}}^{\rm{ + }}$	1-star	2.09010	-18.2650	-
	3-stars	1.71621	-15.3898	-
	5-stars	4.26748	-36.4802	-
	neutral	12.84160	-30.6268	-
	all	1.76362	-15.6848	-
$\omega _{{\rm{CDP1}}}^{\rm{ + }}$	1-star	1.44274	-9.7232	-
	3-stars	1.41243	-8.4555	-
	5-stars	5.87687	-32.0067	-
	neutral	26.43710	-21.5261	-
	all	2.33240	-13.1769	-
$\omega _{{\rm{CDP2}}}^{\rm{ + }}$	1-star	3.55743	-12.1201	-
	3-stars	5.80997	-16.4733	-
	5-stars	13.70820	-37.5049	-
	neutral	37.74420	-22.8075	-
	all	5.09763	-14.5782	-
${\omega _{{\rm{TD}}}}$	1-star	-10.99610	51.4882	-
	3-stars	6.62658	-38.3932	-
	5-stars	7.75751	-42.8484	-
	neutral	12.82170	-30.4143	-
	all	3.01701	-17.8939	-
${{\tilde \omega }_{\rm{P}}}$	1-star	-1.69823	69.3941	1.03686
	3-stars	0.04583	-35.7190	5.28113
	5-stars	0.04319	-42.3104	6.16113
	neutral	3.16158	-33.4227	0.66305
	all	0.08921	-19.9208	2.59209

With the parameters given in Table 1, it is easy to evaluate the quality of the different fits (see Table 3). In general, good correlations are obtained in all cases, with the exception of the 3-stars family. This is probably because the 3-star family is the most structurally diverse group. For example, here we can find heterosubstituted carbocations (like flavylium and 1, 3-dithianylium), methal-stabilized carbocations (like the 1, 3-diphenylpropyn-1-ylium-Cr(CO)₃), benzhydril cations (like pop(tol) CH⁺ and (pfp)₂CH⁺), cyclic conjugated carbocations (like the tropylium ion), and other π-delocalized carbocations (like 1, 1, 3-triphenylallylium and 1-(4-chlorophenyl) cyclopent-2-enylium). In this case, even the best fit (given by ω_TD) still results in an RMSD of 2.950, and a correlation coefficient lower than 0.6. This is a strong indication that to predict the electrophilicity of a given compound we should use a model trained on a set of molecules with similar characteristics. An alternative, in the case we are especially interested in structurally diverse sets, could be to include information regarding condensed local reactivity descriptors ^117-120.

表 3

表 3  Average error (AE), mean absolute error (MAE), root mean square deviation (RMSD), and correlation coefficient (R²) for all the c-DFT-based electrophilicity measures studied.

Table 3.  Average error (AE), mean absolute error (MAE), root mean square deviation (RMSD), and correlation coefficient (R²) for all the c-DFT-based electrophilicity measures studied.

下载: 导出CSV

Electrophilicity measure	Group	AE (×105)	MAE	RMSD	R2
${\omega _{{\rm{parabolic}}}}$	1-star	-430.40	0.043	0.048	0.999933
	3-stars	-3.29	3.025	3.348	0.472271
	5-stars	-3.89	0.654	0.809	0.979628
	neutral	-3.93	1.660	2.289	0.968319
	all	2.02	2.993	3.542	0.775577
$\omega _{{\rm{GCV}}}^{\rm{ + }}$	1-star	-6.24	0.263	0.300	0.997419
	3-stars	-3.37	3.303	3.717	0.349430
	5-stars	-2.81	0.907	1.143	0.959311
	neutral	7.12	1.898	2.368	0.966093
	all	2.63	3.076	3.698	0.755331
$\omega _{{\rm{CDP1}}}^{\rm{ + }}$	1-star	0.13	0.246	0.277	0.997809
	3-stars	-0.87	3.430	3.973	0.256914
	5-stars	-1.16	1.266	1.646	0.915585
	neutral	2.10	2.760	3.250	0.936163
	all	0.04	3.356	3.967	0.718476
$\omega _{{\rm{CDP2}}}^{\rm{ + }}$	1-star	1.72	0.238	0.266	0.997971
	3-stars	1.24	3.289	3.683	0.361174
	5-stars	-9.99	0.792	0.983	0.969887
	neutral	0.95	2.655	3.083	0.942530
	all	-5.20	3.122	3.778	0.744672
${\omega _{{\rm{TD}}}}$	1-star	8.78	0.188	0.204	0.998807
	3-stars	4.06	2.470	2.950	0.590157
	5-stars	4.29	0.515	0.647	0.986951
	neutral	3.49	1.900	2.363	0.966248
	all	3.83	3.101	3.644	0.762501
${{\tilde \omega }_{\rm{P}}}$	1-star	-169.20	0.002	0.002	1.000000
	3-stars	-2.97	2.643	3.038	0.565367
	5-stars	1.09	0.514	0.639	0.987296
	neutral	-4.03	1.690	2.256	0.969225
	all	-1.26	2.978	3.499	0.781001

This insight is supported by the very good fits obtained for the 5-stars and neutral groups (particularly, when we consider the values of R²). These families are mostly formed by structurally similar compounds: benzhydryl cations in the 5-star case, and quinone methides (with a few alkenes substituted with electron-withdrawing groups) in the neutral case. The structural variety of the neutral family is reflected instead in the RMSD values, which are significantly higher than those found in the 5-star case. Overall, all the electrophilicity measures studied have very similar performances over these sets. Nonetheless, in terms of RMSD and R², the ω_CDP1⁺ and ω_CDP2⁺ indices are outperformed by others. These trends are preserved when we perform the fits including all the studied molecules. Finally, we should note that the simple model based on ω_TD consistently ranks among the best. Recall that is simply equal to the electron affinity; it is not surprising that the electrophilicity and the electron affinity are closely related.

5.2   SR results

In Tables 4 and 5 we show the expressions obtained using GE and GP for all families, respectively. In cases where the optimization algorithm resulted in more than one expression, we will only consider the two that have the highest correlation coefficient values. All the solutions with tree lengths up to 4 are explicitly presented. Additionally, we will also discuss some examples of more complex expressions obtained with greater tree lengths.

表 4

表 4  Equations obtained using grammatical evolution.

Table 4.  Equations obtained using grammatical evolution.

下载: 导出CSV

Group	Tree length	Equation code	Equation
1-star	3	ge_13a	$22.89987 - 1.05725{A^2}$
	4	ge_14a	$164.68588 - 30.71200\left( {A + \frac{A}{{\rm{I}}}} \right)$
3-stars	3	ge_33a	$- 20.61898 + 0.61253{A^2}$
	4	ge_34a	$- 14.64513 + 0.07463{A^3}$
5-stars	3	ge_53a	$- 42.84836 + 7.75751A$
	4	ge_54a	$- 46.17069 + 7.67312\left( {A + \frac{A}{{\rm{I}}}} \right)$
neutral	3	ge_N3a	$- 30.41432 + 12.82174A$
	4	ge_N4a	$- 29.00367 - 13.69916\left( { - A + \frac{A}{{\rm{I}}}} \right)$
all	3	ge_A3a	$- 14.29564 + 0.42896{A^2}$
	4	ge_A4a	$- 13.01861 + 0.06894{A^3}$
The equation codes reflect the method (ge: grammatical evolution), the group ("1": 1-star, "3": 3-stars, "5": 5-stars, "N": neutral, "A": all), the three length (3 or 4), and an identifier that indicates if the equation corresponds to the best ("a") or second best ("b") value of R2.

表 5

表 5  Equations obtained using genetic programming.

Table 5.  Equations obtained using genetic programming.

下载: 导出CSV

Group	tree length	equation code	Equation
1-star	3	gp_13a	$22.89987 - 1.05725{A^2}$
	4	gp_14a	$266.68256 - 40.14199\left( {1.20624A + \frac{{1.47614A}}{{\rm{I}}}} \right)$
		gp_14b	$- 4.81713 - \frac{{0.08337{\rm{I}}}}{{11.35373 - 0.61098{\rm{I}}}}$
3-stars	3	gp_33a	$- 20.61898 + 0.61253{A^2}$
	4	gp_34a	$- 2.56271 + \frac{{30.02478\left( { - 15.66098 + 2.92280A} \right)}}{{ - 0.59686A + 1.48008{\rm{I}}}}$
		gp_34b	$- 17.79382 + \frac{{6.97869A}}{{5.84563 - 0.61188A}}$
5-stars	3	gp_53a	$- 41.95494 + 2.60385\left( {3.44526A - 0.26105{\rm{I}}} \right)$
		gp_53b	$39.26578 - \frac{{213.70576}}{A}$
	4	gp_54a	$23.56579 + 6.52492\left( { - 0.51781A + 0.82428{\rm{I}} - \frac{{4.91576{\rm{I}}}}{A}} \right)$
		gp_54b	$75.61898 - \frac{{908.67326}}{{5.50625 + 1.18802A}}$
neutral	3	gp_N3a	$- 33.14915 - 10.36197\left( { - 1.01822A - 0.07852{\rm{I}}} \right)$
		gp_N3b	$- 33.14997 + 4.31176\left( {2.44675A + 0.18877{\rm{I}}} \right)$
	4	gp_N4a	$- 20.89515 - 9.87138\left( { - 1.08308A + \frac{{4.69108}}{{\rm{I}}}} \right)$
		gp_N4b	$- 20.89563 + 7.80050\left( {1.37063A - \frac{{5.93618}}{{\rm{I}}}} \right)$
all	3	gp_A3a	$- 14.29564 + 0.42896{A^2}$
	4	gp_A4a	$- 19.18065 + \frac{{86.80017}}{{15.67561 - 1.92879A}}$
		gp_A4b	$- 19.19082 - \frac{{29.28176}}{{ - 5.28135 + 0.64965A}}$
The equation codes reflect the method (gp: genetic programming), the group ("1": 1-star, "3": 3-stars, "5": 5-stars, "N": neutral, "A": all), the three length (3 or 4), and an identifier that indicates if the equation corresponds to the best ("a") or second best ("b") value of R2.

Unsurprisingly, given the choice to select between I or A to construct a function to fit E_Mayr, both algorithms choose the latter. However, here we see much more general functional forms including the electron affinity. For example, now we have a greater variety of terms with different powers of A (ranging from -1 to +3). In all these cases (except in the pathological 1-star family), the final equations are monotonically increasing functions of A.

The results of the different statistical measures corresponding to the SR fits can be found in Table 6. Overall, the performance of these methods is very similar to the c-DFT electrophilicity measures. For example, for the 3-stars family, the RMSD and R² values obtained both by GE and GP algorithms are virtually identical to those corresponding to the best c-DFT indices (ω_TD and ${\tilde \omega _{\rm{P}}}$). For this family, the quality of the fit remains essentially unaltered if we increase the tree length using the GE algorithm (increasing the tree length up to 7 only increases R² up to 0.539182). On the other hand, increasing the tree length does improve the fit quality when we use GP. In this case, a tree length of 7 gave an R² of 0.703165. The corresponding expression (which, in our notation, would be equation gp_73a) is:

表 6

表 6  Average error (AE), mean absolute error (MAE), root mean square deviation (RMSD), and correlation coefficient (R²) for all the SR-based electrophilicity measures studied.

Table 6.  Average error (AE), mean absolute error (MAE), root mean square deviation (RMSD), and correlation coefficient (R²) for all the SR-based electrophilicity measures studied.

下载: 导出CSV

Group	Equation Code	AE (×105)	MAE	RMSD	R2
1-star	ge_13a	7.82	0.187	0.202	0.895135
	gp_13a	7.82	0.187	0.202	0.895135
	ge_14a	442.77	0.093	0.010	0.974507
	gp_14a	15.73	0.000	0.000	1.000000
	gp_14b	-11.74	0.000	0.000	1.000000
3-stars	ge_33a	-0.38	2.455	2.937	0.533955
	gp_33a	-0.38	2.455	2.937	0.533955
	ge_34a	-56.10	2.442	2.931	0.535784
	gp_34a	-4.19	2.463	2.898	0.546226
	gp_34b	15.04	2.428	2.926	0.537606
5-stars	ge_53a	0.29	0.515	0.647	0.985249
	gp_53a	-9.50	0.507	0.622	0.986389
	gp_53b	0.41	0.480	0.631	0.985974
	ge_54a	-1.49	0.499	0.619	0.986486
	gp_54a	13.15	0.401	0.520	0.990477
	gp_54b	-14.91	0.417	0.579	0.988212
neutral	ge_N3a	-0.31	1.896	2.363	0.778481
	gp_N3a	-2.90	1.688	2.263	0.796833
	gp_N3b	-9.76	1.688	2.263	0.796833
	ge_N4a	-0.27	1.761	2.287	0.792417
	gp_N4a	-4.24	1.656	2.242	0.800661
	gp_N4b	-2.88	1.656	2.242	0.800661
all	ge_A3a	-8.68	2.689	3.325	0.718186
	gp_A3a	-8.68	2.689	3.325	0.718186
	ge_A4a	16.78	2.355	3.119	0.752041
	gp_A4a	-4.98	2.305	3.037	0.764869
	gp_A4b	-8.51	2.304	3.037	0.764869

Nonetheless, requiring an extremely complicated expression to get a semi-quantitative fit is a reminder of the difficulty of predicting experimental electrophilicities with a single expression over a structurally diverse set of compounds.

For the 5-stars family the SR results are excellent, and increasing the tree length in the latter case does not produce a noticeable improvement in the quality of the fit. On the other hand, it is interesting to note that when we work with all the molecules at the same time, the c-DFT indices perform slightly better than the SR expressions. However, when we increase the tree length both GP and GE algorithms provide significantly better solutions. For example, the equation ge_A7a would be:

and it corresponds to an R² of 0.85897. Its GP analog, namely, equation gp_A7a is:

which in turn has an R² of 0.87667.

These results indicate that we need to increase the flexibility of SR methods when the complexity of the system under study increases. This is to be expected since these methods need to "learn" the physical insights that are already built in the c-DFT descriptors.

6.   Conclusions

In this work we have assessed how different theoretical electrophilicity measures can be used to predict the experimental values of this index determined by Mayr et al. The different c-DFT indices used provided acceptable results, as long as they were applied to molecules with similar structures. From a quantitative point of view, it is reassuring that a measure as simple as the electron affinity, ω_TD, provides results of comparable (if not better) quality than models based on more complex indices. This means that we can predict experimental electrophilicity values without computing the ionization energy, which reduces the computational effort.

Perhaps the biggest appealing of the c-DFT approach is that one uses simple expressions, with chemically-motivated terms. In this way, we gain valuable insights into the factors that determinate the electrophilicity at a molecular level. However, we noticed that there are cases when the c-DFT measures were unable to provide the correct qualitative ordering corresponding to the experimental electrophilicity. This serves as an indication that still more effort is needed to improve the c-DFT descriptors, so they can provide a more realistic description of the electron uptake processes. For example, local reactivity indicators are probably needed in order to compare the electrophilicities of diverse families of molecules that include several different types of electrophilic reactive sites (A similar observation has been made for quantitative studies of the quality of leaving groups ²⁸).

To the best of our knowledge, we are reporting the first application of symbolic regression techniques to the prediction of experimental electrophilicity values and the first study of symbolic regression within the conceptual DFT framework. Symbolic regression has the advantage that the resulting models can be improved systematically, by removing restrictions from the optimization process. We showed that this is a key factor for making adequate predictions over relatively complex sets. The simplest models once again proved that the best way to predict the electrophilicity, which often only requires calculating the electron affinity. In general, the expressions obtained using GE are simpler than those obtained using GP, but in most cases the corresponding fits are equivalent. Finally, it should be remarked that, while attractive, these models could be prone to overfitting.

In general, when used carefully, the different theoretical methods studied here provide adequate predictions of the experimental electrophilicity. Further studies are necessary, however, to assess their utility over a broader range of compounds, including their validation using independent "test" sets of molecules. It is also interesting to check whether using the exact electrophilicity formula given in Eq. (19) will improve the quantitative predictions.

References

1. [1]

   (1) Zhang, A. P.; Zhang, J. Z. Acta Phys. -Chim. Sin. 2010, 26, 1337. [张爱平, 张进治. 物理化学学报, 2010, 26, 1337.] doi: 10.3866/PKU.WHXB20100533
   

2. [2]

   (2) Lin, X.; Lü, P.; Guan, Q. F.; Li, H. B.; Li, H. J.; Cai, J.; Zou, Y. Acta Phys. -Chim. Sin. 2012, 28, 1978. [林雪, 吕鹏, 关庆丰, 李海波, 李洪吉, 蔡杰, 邹阳. 物理化学学报, 2012, 28, 1978.] doi: 10.3866/PKU.WHXB201205172
   

3. [3]

   (3) Liu, Y. F.; Ma, X. G.; Yi, X.; Zhu, Y. F. Acta Phys. -Chim. Sin. 2012, 28, 654. [刘艳芳, 马新国, 易欣, 朱永法. 物理化学学报, 2012, 28, 654.] doi: 10.3866/PKU.WHXB201112232
   

4. [4]

   (4) Zhang, L.;Wang, H.; Chen, Z.;Wong, P. K.; Liu, J. Appl. Catal. B: Environ. 2011, 106, 1.
   

5. [5]

   (5) Zhang, L.; Zhu, Y. Catal. Sci. Technol. 2012, 2, 694. doi: 10.1039/c2cy00411a
   

6. [6]

   (6) Kudo, A.; Hijii, S. Chem. Lett. 1999, 10, 1103.
   

7. [7]

   (7) Shang, M.;Wang,W.; Sun, S.; Zhou, L.; Zhang, L. J. Phys. Chem. C 2008, 112, 10407. doi: 10.1021/jp802115w
   

8. [8]

   (8) Shang, M.;Wang,W.; Xu, H. Cryst. Growth Des. 2009, 9, 991. doi: 10.1021/cg800799a
   

9. [9]

   (9) Shang, M.;Wang,W.; Ren, J.; Sun, S.;Wang, L.; Zhou, L. J. Mater. Chem. 2009, 19, 6213. doi: 10.1039/b907849e
   

10. [10]

    (10) Amano, F.; Nogami, K.; Abe, R.; Ohtani, B. J. Phys. Chem. C 2008, 112, 9320. doi: 10.1021/jp801861r
    

11. [11]

    (11) Zhang, L.;Wang,W.; Chen, Z.; Zhou, L.; Xu, H.; Zhu,W. J. Mater. Chem. 2007, 17, 2526. doi: 10.1039/b616460a
    

12. [12]

    (12) Shi, R.; Huang, G.; Lin, J.; Zhu, Y. J. Phys. Chem. C 2009, 113, 19633. doi: 10.1021/jp906680e
    

13. [13]

    (13) Fu, H.; Zhang, S.; Xu, T.; Zhu, Y.; Chen, J. Environ. Sci. Technol. 2008, 42, 2085. doi: 10.1021/es702495w
    

14. [14]

    (14) Guo, S.; Li, X.;Wang, H.; Dong, F.;Wu, Z. J. Colloid Interface Sci. 2012, 369, 373. doi: 10.1016/j.jcis.2011.12.007
    

15. [15]

    (15) Duan, F.; Zheng, Y.; Chen, M. Appl. Sci. Res. 2011, 257, 1972.
    

16. [16]

    (16) Zhou, L.; Yu, M.; Yang, J.;Wang, Y.; Yu, C. J. Phys. Chem. C 2010, 114, 18812. doi: 10.1021/jp107061p
    

17. [17]

    (17) Song, X. C.; Zheng, Y. F.; Ma, R.; Zhang, Y. Y.; Yin, H. Y. J. Hazard. Mater. 2011, 192, 186.
    

18. [18]

    (18) Zhang, L.; Man, Y.; Zhu, Y. ACS Catal. 2011, 1, 841. doi: 10.1021/cs200155z
    

19. [19]

    (19) Amano, F.; Yamakata, A.; Nogami, K.; Osawa, M.; Ohtani, B. J. Phys. Chem. C 2011, 115, 16598. doi: 10.1021/jp2051257
    

20. [20]

    (20) Fu, H.; Zhang, L.; Yao,W.; Zhu, Y. Appl. Catal. B 2006, 66, 100. doi: 10.1016/j.apcatb.2006.02.022
    

21. [21]

    (21) Zhang, C.; Zhu, Y. Chem. Mater. 2005, 17, 3537. doi: 10.1021/cm0501517
    

22. [22]

    (22) Amano, F.; Nogami, K.; Ohtani, B. J. Phys. Chem. C 2009, 113, 1536. doi: 10.1021/jp808685m
    

23. [23]

    (23) Chakrabarti, S.; Ganguli, D.; Chaudhuri, S. Physica E 2004, 24, 333.
    

24. [24]

    (24) Maczka, M.; Macalik, L.; Hermanowicz, K.; KepiDski, L.; Tomaszewski, P. J. Raman Spectrosc. 2010, 41, 1059.
    

25. [25]

    (25) Graves, P. R.; Hua, G.; Myhra, S.; Thompson, J. G. J. Solid State Chem. 1995, 114, 112. doi: 10.1006/jssc.1995.1017
    

26. [26]

    (26) Bordun, O. M. Inorg. Mater. 1998, 34, 1270.