Probing the nature of Lewis acid sites on oxide surfaces with ³¹P(CH₃)₃ NMR: a theoretical analysis

Abstract

The characterization of catalytic oxide surfaces is often done by studying the properties of adsorbed probe molecules. The ³¹P NMR chemical shift of adsorbed trimethylphosphine, P(CH₃)₃ or TMP, has been used to identify the presence of different facets in oxide nanocrystals and to study the acid–base properties of the adsorption sites. The NMR studies are often complemented by DFT calculations to provide additional information on TMP adsorption mode, bond strength, etc. So far, however, no systematic study has been undertaken in order to compare on the same footing the chemical shifts and the adsorption properties of TMP on different oxide surfaces. In this work we report the results of DFT+D (D = dispersion) calculations on the adsorption of TMP on the following oxide surfaces: anatase TiO₂(101) and (001), rutile TiO₂(110), tetragonal ZrO₂(101), stepped ZrO₂(134) and (145) surfaces, rutile SnO₂(110), (101) and (100), wurtzite ZnO(101̄0), and cubic CeO₂(111) and (110). Beside the stoichiometric surfaces, also reduced oxides have been considered creating O vacancies in various sites. TMP has been adsorbed on top of variously coordinated Lewis acid cation sites, with the aim to identify, also with the support of machine learning algorithms, trends or patterns that can help to correlate the ³¹P chemical shift with physico-chemical properties of the oxide surfaces such as adsorption energy, Bader charges, cation–P distance, work function, etc. Some simple correlation can be found within the same oxide between the ³¹P chemical shift and the adsorption energy, while when the full set of data is considered the only correlation found is with the net charge on the TMP molecule, a descriptor of the acid strength of the adsorption site.

---

Introduction

Nuclear magnetic resonance (NMR) is a powerful and sensitive technique that can be used to study catalytic sites and functional groups anchored on a surface.¹ To this end one can directly investigate the NMR properties of the ions constituting the material, metal cations² and oxygen anions,^3,4 or use simple adsorbed molecules to probe the nature of the surface,⁵ in particular the acid or basic character of the adsorption sites. In this context, trimethylphosphine, P(CH₃)₃ or TMP, has been used to investigate oxide surfaces and characterize their acidity by studying the chemical shift of the ³¹P nucleus of adsorbed TMP.⁶ Probably one of the first studies in this field is that of Lunsford and co-workers who in 1984 studied the acidic properties of a zeolite.⁷ Since then, the method has been adopted for other oxides. More recently, Peng and co-workers⁸ have monitored the nature of various facets of anatase TiO₂ (a-TiO₂) nanoparticles using ³¹P NMR of adsorbed TMP as a surface probe. They showed that ³¹P NMR is capable of differentiating the facets of the titania nanoparticles. A peak at −36 ppm and a shoulder at −29 ppm were assigned to the interaction of TMP with Ti⁴⁺ cations of the (101) and (001) facets, respectively. They also found that the TMP adsorption energy, as obtained from DFT calculations, correlates with the NMR chemical shift.⁸ Another oxide that has been investigated with this approach is zirconia. Isolated TMP molecules, and TMP coadsorbed with CO₂, were used to study the basic properties of ZrO₂ using ³¹P NMR.⁹ The ³¹P chemical shift of TMP adsorbed on zirconia has been discussed also in other studies, finding values in the range −28 to −50 ppm.^10–13

SnO₂ nanosheets and nano-shuttles were investigated using TMP probe molecules and ³¹P NMR to study the facet-dependent acidity.¹⁴ This work showed that the (001), (101), (110), and (100) facets can be differentiated by ³¹P NMR chemical shifts of adsorbed TMP. Also in this case, a linear correlation has been reported with the calculated adsorption energies of TMP on the various surfaces and the measured ³¹P chemical shift.¹⁴

Even ZnO plates, rods and powder samples have been investigated by means of ³¹P NMR of adsorbed TMP combined with DFT calculations.¹⁵ DFT was also used to study the interaction energy and the effect of oxygen vacancies on both nonpolar (101̄0) and polar (0002) surfaces of ZnO nanoparticles.¹⁵ As for other oxides, a linear relationship has been found between the adsorption energy of TMP and the calculated ³¹P NMR chemical shifts.

Finally, ³¹P NMR of TMP probe molecules has been used to investigate various facets of CeO₂ octahedral, rod, and cube samples, as well as the concentrations of the corresponding sites.¹⁶ The dominant ³¹P NMR signal of adsorbed TMP is observed at −33, −47.5, and −58 ppm for octahedral, rod, and cube structures with dominant exposed (111), (110), and (100) facets, respectively. The adsorption energy of TMP on the regular CeO₂(111), (110), and (100) surfaces and on the hydroxylated (111) and (100) surfaces was calculated at the DFT level and plotted against the ³¹P NMR chemical shifts. Notice however that these correlations are always found within the same acid center.¹⁷ The question whether similar correlations can be established for different oxides or not, so as to provide a scale of acidity of the surface from the analysis of the ³¹P chemical shift, remains open.

This brief summary shows the wide use of the ³¹P NMR technique to study oxide surfaces and the different morphologies of oxide nanostructures using adsorption of TMP probe molecules. In several of these studies the experimental results are complemented by specifically designed DFT calculations, showing a good potential of this approach to assign specific adsorption sites to a given NMR signal. However, so far, no systematic study of different oxide surfaces using the same computational approach has been reported. In the present study we compare, using the same DFT approach, the NMR chemical shifts of ³¹P for TMP molecules adsorbed on anatase and rutile TiO₂, tetragonal ZrO₂, rutile SnO₂, cubic CeO₂ and wurtzite ZnO. For some oxides different facets have been considered, and beside regular surfaces, in a few cases we also studied stepped surfaces with low-coordinated ions. Coverage effects on the chemical shifts have also been addressed. Finally, oxygen vacancies have been created on the surface of these oxides, and the consequences on the TMP adsorption and on the ³¹P NMR chemical shifts have been discussed.

Computational method

Density functional theory (DFT) calculations have been performed using the Vienna Ab Initio Simulation Package (VASP 5.4.4).^18–20 For the exchange–correlation functional we used the Perdew, Burke and Ernzerhof (PBE) formulation,²¹ and to correct the self-interaction error, the PBE+U^22,23 approach has been adopted using the following set of Hubbard U parameters: U = 3 eV for the 3d states of Ti,²⁴U = 4.7 eV for the 3d states of Zn,²⁵U = 4 eV for the 4d states of Zr,²⁶U = 3.5 eV for the 4d states of Sn,²⁷U = 4 eV for the Ce 4f states.²⁷ Dispersion was included by means of the Grimme D3 approach.^28,29 To describe the effect of the core electrons we used the projector augmented wave (PAW) method.^30,31 The valence electrons explicitly considered are H (1s), C (2s, 2p), O (2s, 2p), P (3s, 3p), Ti (3s, 3p, 3d, 4s), Zn (3d, 4s), Zr (4s, 4p, 4d, 5s), Sn (4d, 5s, 5p), and Ce (5s, 6s, 5p, 4f, 5d). The energy cutoff for plane waves was set to 400 eV, and the optimizations were performed using the conjugate gradient scheme until the change in total energy between successive steps was less than 10⁻⁵ eV.

We considered a-TiO₂(101) surfaces with [3 × 1] and [2 × 1] supercells with 5 layers of Ti and 10 layers of O (5 TiO₂ trilayers, Ti₆₀O₁₂₀ and Ti₄₀O₈₀ formula, respectively). Another titania facet, a-TiO₂(001), was built with a [2 × 2] supercell with 6 layers of Ti and 12 layers of O (6 TiO₂ trilayers, Ti₂₄O₄₈ formula). [3 × 2] and [2 × 1] supercells models were used for tetragonal ZrO₂(101) with 5 layers of Zr and 10 layers of O (5 ZrO₂ trilayers, Zr₆₀O₁₂₀ and Zr₂₀O₄₀ formula, respectively). The zirconia stepped surfaces, ZrO₂(134) and ZrO₂(145), were described as in our previous work.³² We described the SnO₂(110) surface with [3 × 2] and [2 × 1] supercells with 5 layers of Sn and 10 layers of O (5 SnO₂ trilayers, Sn₆₀O₁₂₀ and Sn₂₀O₄₀ formula, respectively). The same supercells were considered for the rutile r-TiO₂(110) surface. SnO₂(101) and (100) surfaces with [2 × 2] supercells with 5 and 7 layers of Sn and 10 and 14 layers of O, respectively (5 and 7 SnO₂ trilayers, Sn₄₀O₈₀ and Sn₂₈O₅₆ formula, respectively) were considered. For CeO₂(111) [2 × 2] and [1 × 1] supercells were used with 5 layers of Ce and 10 layers of O (5 CeO₂ layers, Ce₈₀O₁₆₀ and Ce₂₀O₄₀ formula, respectively). The [2 × 1] and [1 × 1] supercells of CeO₂(110) were built with 9 atomic layers of Ce and O (9 CeO₂ layers, Ce₃₆O₇₂ and Ce₁₈O₃₆ formula, respectively). Finally, we studied two wurtzite ZnO(101̄0) models with [3 × 2] and [2 × 2] supercells with 14 atomic layers of Zn and O (7 ZnO layers, Zn₈₄O₈₄ and Zn₅₆O₅₆ formula, respectively). In all cases the slabs were separated by more than 15 Å of vacuum.

We optimized the bulk lattice parameters of each oxide using a kinetic energy cutoff of 600 eV and a 8 × 8 × 8 Monkhorst–Pack k-point grid. All atomic positions were fully optimized with convergence obtained when the ionic forces were smaller than |0.01| eV Å⁻¹. Atomic charges have been determined using the Quantum Theory of Atoms in Molecules (QTAIM) or Bader analysis.³³

The work function was computed as the energy of the vacuum level (determined applying a dipole correction to the unit cell) with respect to Fermi level (in this case the top of the valence band of the oxide).

The adsorption energies of TMP have been calculated as:

E_ads = E(TMP/MO_x) − E(TMP)(g) − E(MO_x) (1)

where TMP = P(CH₃)₃ and MO_x refers to the oxide surface. The formation energy of an O vacancy (E_f) is calculated according to the following equation:

E_f = E(MO_x−1) + ½E(O₂)(g) − E(MO_x) (2)

Notice that at the PBE level the binding energy of the O₂ molecule is overestimated, and this may reflect in some uncertainty in the computed O vacancy formation energies. However, the scope of this work is not to provide accurate values for this property but rather to assess the possibility to detect the presence of vacancies via³¹P NMR.

The isotropic chemical shift (δ_iso) is computed as δ_iso = σ_ref − σ_calc, where σ_ref is the shielding value of the reference and σ_calc is the chemical shielding obtained in VASP. VASP calculates the chemical “shifts” for non-metallic crystalline systems using the linear response method of Yates, Pickard, and Mauri.^34,35 For ³¹P NMR, on the basis of experimental reports, we consider H₃PO₄ and (NH₄)H₂PO₄ compounds as references since these are adopted in experimental studies. The VASP chemical shieldings for these compounds are 294 ppm and 283 ppm, respectively. With our setup the computed chemical shift for a free, gas-phase TMP molecule is of −60 ppm to be compared with an experimental value of −62 ppm, Table 1.³⁶

Table 1 Work function (Φ, eV) of the clean surfaces, isomers of adsorbed P(CH₃)₃, adsorption energy (E_ads, eV), Bader charge (q, |e|), bond distance of the P atom from the surface cation M (R_M–P, Å), ³¹P NMR chemical shielding (σ_calc, ppm), isotropic chemical shifts (δ_iso, ppm) and experimental chemical shifts (δ_iso, ppm) of adsorbed P(CH₃)₃

		Φ	Isomer	E ads	Q			R M–P	σ calc	δ iso	δ exp	Ref.
					P	M	P(CH3)3
a Reference: H3PO4. b Reference: (NH4)H2PO4.
Free P(CH3)3		—	—	—	1.20	—	0.00	—	354	−60^a	−62^a	34
a-TiO2(101)	[3 × 1]	6.98	1	−1.49	1.27	2.16	0.21	2.644	310	−16^a	−36^a	8
	[2 × 1]	6.93	2	−1.48	1.23	2.26	0.19	2.631	310	−16^a
			1	−1.46	1.27	2.26	0.19	2.649	313	−19^a
a-TiO2(001)	[2 × 2]	5.93	2	−1.29	1.20	2.23	0.15	2.602	319	−25^a	−29^a
		5.93	1	−1.26	1.29	2.22	0.16	2.613	320	−26^a
r-TiO2(110)	[3 × 2]	6.93	2	−1.76	1.23	2.25	0.24	2.659	298	−4^a	—	—
	[2 × 1]	6.89	2	−1.05	1.21	2.25	0.17	2.700	319	−25^a
ZrO2(101)	[3 × 2]	6.42	2	−1.33	1.20	2.68	0.13	2.781	323	−29^a	−28 to −50^a (average = −40)	9–13
		6.42	1	−1.32	1.23	2.68	0.13	2.796	323	−29^a
	[2 × 1]	6.48	2	−1.30	1.25	2.68	0.10	2.789	328	−34^a
ZrO2(134)	[1 × 1]	5.87	2	−1.41	1.18	2.63	0.10	2.812	332	−38^a
ZrO2(145)	[1 × 1]	6.31	1	−1.43	1.22	2.65	0.12	2.809	333	−39^a
		6.31	2	−1.44	1.21	2.63	0.10	2.799	332	−38^a
SnO2(110)	[3 × 2]	7.15	1	−2.46	1.37	2.19	0.41	2.617	284	−1^b	−23^b	14
	[2 × 1]	6.97	1	−2.13	1.35	2.25	0.30	2.636	302	−19^b
SnO2(101)	[2 × 2]	6.79	1	−2.11	1.33	2.22	0.37	2.585	277	6^b	−14^b
SnO2(100)	[2 × 2]	7.69	1	−1.57	1.31	2.23	0.27	2.698	302	−19^b	−29^b
CeO2(111)	[2 × 2]	5.91	2	−0.95	1.31	2.32	0.14	3.039	314	−20^a	−33^a	16
	[1 × 1]	5.91	2	−0.95	1.31	2.32	0.10	3.042	322	−28^a
		5.91	1	−0.93	1.29	2.33	0.10	3.079	323	−29^a
CeO2(110)	[2 × 1]	4.16	1	−0.85	1.25	2.30	0.07	3.148	340	−46^a	−47.5^a
	[1 × 1]	4.22	1	−0.84	1.21	2.30	0.05	3.141	344	−50^a
ZnO(101̄0)	[3 × 2]	5.82	1	−1.68	1.37	1.07	0.27	2.339	335	−41^a	−43^a	15
	[2 × 2]	5.87	1	−1.63	1.33	1.03	0.24	2.347	337	−43^a
		5.87	2	−1.60	1.37	1.05	0.24	2.354	340	−46^a

---

A search of correlations between some of the computed properties has been attempted using the classical approach to plot two variables. Beside this, a multiple linear regression method (MLR) has been used to fit (X₁, X₂, X₃,…, y) with linear relationships to search for non-obvious correlations. Eqn (3) and (4) are multiple linear regression for one sample (X_n, y_n) and multiple samples (X, y), respectively. In eqn (4), X is a N × (d + 1) matrix, where N is the number of data points and d is the number of parameters, so each row of X corresponds to one sample data [1, X_n1, X_n2, …, X_nd]. W is a column vector [b, W₁, W₂,…, W_d]^T, where b is a scalar, representing the linear offset. The basic idea of multiple linear regression is to reduce the error between the fitted results y_fit and the original set of values y.

y = XW^T (4)

Results and discussion

1. TMP adsorption on stoichiometric oxide surfaces

We start the discussion from TMP adsorption on anatase and rutile TiO₂, Fig. 1(a–c) and Table 1. TMP always binds preferentially with a P atom pointing towards a surface Ti⁴⁺ cation with binding energies that, at low coverage, go from −1.76 eV on r-TiO₂(110) to −1.49 eV on a-TiO₂(101) and −1.29 eV on a-TiO₂(001). Notice that slightly different isomers can form, with very similar stabilities, Table 1. This depends on the rotation of the methyl groups with respect to the surface. At higher coverage we found a significant change in binding only for the r-TiO₂(110) surface where E_ads becomes −1.05 eV. In all cases the distance of the TMP molecule from the surface Ti atom is similar, around 2.6 Å, and increases slightly for higher coverages, Table 1. The computed ³¹P chemical shifts, referred to H₃PO₄, show a dependence both on the surface and on the coverage. On a-TiO₂(101) at low coverage δ_iso is −16 ppm and it becomes −19 ppm at higher coverage. Experimentally, a chemical shift of −36 ppm has been assigned to TMP adsorbed on this surface.⁸ A more pronounced chemical shift, −25 ppm, is computed for the a-TiO₂(001) surface, to be compared with the value of −29 ppm reported experimentally.⁸ On this surface virtually no change is found by increasing the coverage.

Fig. 1 Most stable isomers of adsorbed TMP on the different oxides and facets considered in this work.

If we plot the computed ³¹P chemical shift versus the adsorption energies we found a linear behavior with good correlation (R² = 0.97, Fig. 2); of course, one has to consider that this correlation is obtained with only five values; furthermore, the changes in chemical shift are rather small (from −16 to −26 ppm). Nevertheless, a trend is present and shows more negative values for weaker TMP–surface interaction.

Fig. 2 Plot of computed ³¹P chemical shift versus adsorption energy for individual oxide surfaces: (a) CeO₂, (b) a-TiO₂, (c) ZrO₂.

On the rutile structure, r-TiO₂(110), the chemical shift shows a pronounced dependence on the coverage, and goes from −4 ppm at low coverage to −25 ppm at high coverage, Table 1. No experimental data seem to exist for this surface. When we try to add the two points computed for r-TiO₂(110) to the plot of Fig. 2 the regression deteriorates and R² becomes 0.85. This seems to indicate that a correlation between adsorption energy and chemical shift exists, but only within a given polymorph.

We move now to the zirconia surface, Fig. 1d. Here we considered the most stable (101) surface of tetragonal zirconia, t-ZrO₂(101), and we studied the dependence of ³¹P chemical shift of adsorbed TMP on both coverage and coordination of the cation sites. To this end we considered the stepped ZrO₂(134) and ZrO₂(145) surfaces,³⁷Fig. 1e, f and Table 1. On t-ZrO₂(101) TMP binds to Zr_7c⁴⁺ ions with E_ads = −1.33 eV, a value that does not change with coverage, Table 1. On the low-coordinated Zr_6c sites of the stepped surfaces the adsorption is only slightly higher, −1.41 eV or −1.44 eV. On all these systems the distance of TMP from the surface Zr ion is similar and is around 2.80 Å, Table 1. The chemical shifts of ³¹P go from −29 ppm on t-ZrO₂(101) at low coverage, to −34 ppm by increasing the coverage. So, the density of the TMP molecules weakly affects both the adsorption energy and the chemical shift. More pronounced is the effect of coordination of the Zr ion, as the chemical shift becomes −38/−39 ppm on the two stepped zirconia surfaces. Experimentally, a range of −28 to −50 ppm has been reported for TMP on this surface,^9–13 suggesting a possible role of low-coordinated sites and a relatively high coverage of TMP molecules in the measured samples.

In Fig. 2 we report the correlation between the computed ³¹P chemical shifts and adsorption energies (6 values). The correlation is not as good as found for the a-TiO₂ system (R² = 0.70, Fig. 2). On the other hand, also in this case a trend is found but the slope of the curve is opposite to the titania case, since here more negative values of the chemical shift correspond to more strongly bound TMP molecules, Fig. 2.

The next oxide analyzed is SnO₂. Here we considered the rutile structure and the (110), (101) and (100) surfaces, Fig. 1g–i. Notice that for this oxide the reference of the chemical shifts is not the H₃PO₄ molecule as in other cases, but the (NH₄)H₂PO₄ ammonium compound in order to compare our results with existing measured data¹⁴ (in order to compare with H₃PO₄ one has to add +11 ppm to the values of Table 1). On the r-SnO₂(110) surface TMP is strongly bound with E_ads = −2.46 eV, a value that decreases to −2.13 eV for higher coverage, Table 1. On the other surfaces TMP is less strongly bound and exhibits E_ads = −2.11 eV on the (101) surface and −1.57 eV on the (100) facet. These values are somewhat larger than those reported by Zhang et al.¹⁴ because no dispersion was included in their study. Furthermore, their calculations refer to high coverage adsorption, and we have shown that the adsorption energy is sensitive to the coverage (see above). A great variability of the adsorption energy is found as a function of the surface structure. This is partly reflected in the Sn–P distance that shows a larger variation than in previous cases: it goes from a minimum of 2.585 Å on SnO₂(101) to 2.698 Å on the (100) surface, Table 1.

The variability in adsorption energy and Sn–P distance is only partly reflected in the NMR chemical shifts. On r-SnO₂(110) the ³¹P chemical shifts (referred to (NH₄)H₂PO₄) is of −1 ppm at low coverage. This becomes −19 ppm at high coverage, a value close to the experimental measurement of −23 ppm.¹⁴ On the (101) and (100) surfaces, and for low coverage, the chemical shifts are 6 and −19 ppm, respectively, to be compared with experimental values of −14 and −29 ppm attributed to TMP adsorbed on these surfaces. Thus, the results show a significant dependence both on facet exposed and density of TMP molecules on the surface. In general, the agreement with the measured values is better when higher coverage of TMP is considered. Due to the small number of points, we did not attempt to correlate chemical shift and adsorption energy.

Two surfaces of CeO₂ have been considered, the (111) and the (110), Fig. 1j and k, and Table 1. Of all the oxides considered, CeO₂ shows the smallest binding with TMP and the longest cation–TMP distance: the binding energies, close to those reported by Tan et al.,¹⁶ go from −0.84 to −0.95 eV depending on the facet and on the coverage, while the distance goes from a minimum of 3.039 Å on CeO₂(111) (low coverage) to a maximum of 3.148 Å on CeO₂(110) (low coverage). These changes result in different ³¹P chemical shifts of TMP. On CeO₂(111) we compute a value of −20 ppm at low coverage (reference H₃PO₄) and −29 ppm at high coverage (showing once more a dependence on the coverage). These values are in excellent agreement with a peak in the NMR spectra at −33 ppm observed experimentally.¹⁶ On the CeO₂(110) surface the chemical shift is considerably larger in absolute value, −46 ppm at low coverage and −50 ppm at high coverage, again in close agreement with the experiment that shows a peak at −47.5 ppm attributed to TMP adsorbed on the (110) CeO₂ surface.¹⁶ A nice correlation is found when we plot the chemical shifts for the five cases considered versus the strength of the TMP bond to CeO₂ (R² = 0.95) (Fig. 2). Here a weaker bonding corresponds to a more negative chemical shift, the same trend observed for TiO₂ but opposite to what found for ZrO₂.

Finally, only the (101̄0) surface of ZnO has been considered, Fig. 1l. Notice that the main difference in this system compared to the other oxides is the formal charge of the cation, which is +II while it is +IV in the rest of systems considered. We used two supercells that allow us to study the role of coverage. This is relatively minor on the adsorption energy, which goes from −1.68 eV to −1.60 eV going from the [3 × 2] to the [2 × 2] supercell, and modest also on the chemical shift, as this goes from −41 ppm (low coverage) to −46 ppm (high coverage). These values are in close agreement with the reported experimental shift of −43 ppm for TMP on ZnO, Table 1 (reference H₃PO₄).¹⁵ The adsorption energies are about 50% larger than those reported by Peng et al.¹⁵ due to the inclusion of dispersion in our work. The Zn–P distance is the shortest among the cation–P distances of all oxides considered, about 2.35 Å, Table 1. This, however, is not a factor contributing the final value of the chemical shift, as we will show below.

2. TMP adsorption on reduced oxide surfaces

In the previous section we have discussed the adsorption of TMP on the regular, non-defective surfaces of TiO₂, ZrO₂, SnO₂, CeO₂ and ZnO. With the exception of zirconia,³⁸ all the other oxides are classified as reducible oxides characterized by a band gap of 2–4 eV, typical of wide gap semiconductors, and by the relatively easy loss of oxygen under reactive conditions with formation of oxygen vacancies, V_O.³⁹ For these oxides, again with the exception of bulk zirconia, the presence of oxygen vacancies is rather common and often it is difficult to prepare non defective surfaces. In reducible oxides the removal of O atoms results in excess electrons that are transferred to the bottom of the conduction band. This can lead to charge localization with polaron formation.⁴⁰ The most studied cases in this respect are those of TiO_2−x and CeO_2−x oxides where the excess of charge associated to V_O centers localizes on Ti 3d or Ce 4f states. Even ZrO₂, when produced in form of nanoparticles or nanofilms, increases considerably its reducibility and can lose oxygen quite easily, as recently shown both theoretically,^41,42 and experimentally.⁴³

In this respect, it is interesting to study the possibility to use the ³¹P NMR signal of adsorbed TMP to address the effect of the presence of oxygen vacancies on the surface of the oxide. So far, the role of vacancies on the ³¹P chemical shift of TMP has not been studied in a systematic way. We removed an O atom from various positions of the most stable surfaces (see Fig. 1 where the O_nc atoms removed are indicated) and we computed the NMR chemical shift of the TMP probe molecule adsorbed on the same cation shown in the figure, which is next neighbor or in close proximity of the defect center. Notice that it is not our aim to study the relative stability of O vacancies in these oxides (for instance in anatase TiO₂ O vacancies form in the sub-surface, not on the surface), nor the level of electron localization that follows the O removal. Here we only want to see how the ³¹P NMR properties are depending on the presence of these defects. The analysis has been done for the larger supercells in order to avoid spurious effects due to an excessive number of O vacancies. The results, Table 2, are directly compared with those of the same system in absence of vacancies.

Table 2 Formation energy of an oxygen vacancy (E_f, eV), isomers of adsorbed P(CH₃)₃, adsorption energy (E_ads, eV), Bader charges (q, |e|), bond distance of P atom from the surface cation (R_M–P, Å), ³¹P NMR chemical shielding (σ_calc, ppm), isotropic chemical shift (δ_iso, ppm) of the adsorbed P(CH₃)₃. The reference of δ_iso is the chemical shielding of the H₃PO₄ molecule

		E f	Isomer	E ads	q			R M–P	σ calc	δ iso
					P	M	P(CH3)3
a-TiO2(101)	[3 × 1]	—	1	−1.49	1.27	2.16	0.21	2.644	310	−16
	[3 × 1]_VO2c	4.67	1	−1.66	1.21	2.17	0.24	2.567	299	−5
	[3 × 1]_VO3c	5.68	1	−1.60	1.23	2.00	0.16	2.629	317	−23
r-TiO2(110)	[3 × 2]	—	2	−1.76	1.23	2.25	0.24	2.659	298	−4
	[3 × 2]_VO2c	3.05	2	−1.64	1.18	2.24	0.20	2.644	277	17
	[3 × 2]_VO3c	4.90	2	−1.51	1.23	2.20	0.21	2.700	273	21
ZrO2(101)	[3 × 2]	—	2	−1.33	1.20	2.68	0.13	2.781	323	−29
	[3 × 2]_VO	5.97	2	−1.35	1.16	2.39	−0.02	2.745	321	−27
SnO2(110)	[3 × 2]	—	1	−2.46	1.37	2.19	0.41	2.617	284	10
	[3 × 2]_VO2c	2.32	2	−2.11	1.28	2.24	0.33	2.605	304	−10
	[3 × 2]_VO3c	3.47	1	−2.23	1.27	1.97	0.42	2.631	273	21
CeO2(111)	[2 × 2]	—	2	−0.95	1.31	2.32	0.14	3.039	314	−20
	[2 × 2]_VO	2.85	2	−1.08	1.15	2.30	0.07	3.131	333	−41
ZnO(101̄0)	[3 × 2]	—	1	−1.68	1.37	1.07	0.27	2.339	333	−41
	[3 × 2]_VO	3.09	1	−1.47	1.27	0.70	0.24	2.396	334	−40

---

On a-TiO₂(101) both O_2c and O_3c atoms from the top layer have been removed, Fig. 1a. The vacancy formation energy, E_f = 4.67 eV, obviously indicates preferential removal of O_2c with respect to O_3c, E_f = 5.68 eV, Table 2. The adsorption energy of TMP to the adjacent Ti⁴⁺ site is similar to the corresponding regular surface: defective TiO₂E_ads = −1.66 eV (near V_O2c) or −1.60 eV (near V_O3c), respectively, to be compared with the regular surface E_ads = −1.49 eV. The Ti–P distance is affected only for the case of V_O2c, where it becomes 2.567 Å (2.644 Å on the non-defective surface). When the O_2c atom is removed, the ³¹P signal of TMP shifts by 11 ppm, and goes from −16 ppm on the regular surface to −5 ppm on the reduced one. An opposite effect is found when the O_3c atom is removed, as here the ³¹P chemical shift is −23 ppm, with an additional shift of −7 ppm. Notice however that V_O3c centers hardly form due to the large difference in formation energy of the two defects, see Table 2.

On the rutile polymorph, r-TiO₂(110), Fig. 1c, the removal of O has a lower cost than on anatase TiO₂, in particular for the O_2c atoms (E_f = 3.05 eV) compared to the O_3c ones (E_f = 4.90 eV). When a V_O2c is created, the adsorption energy of TMP to an adjacent Ti⁴⁺ site, E_ads = −1.64 eV, is smaller in absolute value than on the regular surface, E_ads = −1.76 eV (the effect is even larger for O_3c, Table 2). No major change in the Ti–P distance is observed, Table 2. Interestingly, the ³¹P chemical shift goes from δ_iso = −4 ppm on regular r-TiO₂(110) to δ_iso = 17 ppm on the reduced surface, r-TiO_2−x(110), with a positive shift of 21 ppm, the same direction found for the anatase polymorph but more pronounced. A similar situation is found for a V_O3c defect since also in this case δ_iso goes from −4 to +21 ppm.

Removing an O_3c atom from the (101) surface of ZrO₂, Fig. 1d, has a cost that, not surprisingly, is almost twice that of r-TiO₂(110) (E_f = 5.97 eV, Table 2). The high band gap of zirconia results in a different electronic structure of the defect, with the excess of charge mostly localized inside the vacancy³⁸ instead of being transferred to the low-coordinated Zr ion. This fact results in virtually no change in the properties of adsorbed TMP: E_ads (−1.33 eV on the regular surface) becomes −1.35 eV, the Zr–P distance from 2.781 Å becomes 2.745 Å, and δ_iso from −29 ppm becomes −27 ppm. Even the charge distribution remains unchanged, Table 2. It seems that on ZrO₂ the NMR signal is insensitive to the presence of O vacancies.

The SnO₂(110) surface has the same structure as r-TiO₂(110) but the removal of oxygen, Fig. 1g, has a much lower cost. E_f(V_O2c) is in fact 2.32 eV, and also E_f(V_O3c), 3.47 eV, is relatively low. When a TMP molecule is adsorbed near a vacancy the adsorption strength is reduced, −2.11 eV versus −2.46 eV on the regular surface (a similar effect is found for O_3c, Table 2). Virtually no change is observed on the Sn–P distance, Table 2. On the other hand, the ³¹P NMR chemical shift is sensitive to the presence of the vacancy, and δ_iso goes from 10 ppm on stoichiometric SnO₂ to −10 ppm on the reduced one, Table 2, with a shift of −20 ppm (here the shifts are given with respect to the H₃PO₄ molecule, as for the other systems). The effect is of similar magnitude to that computed for rutile TiO₂, but of opposite sign! When the O_3c atoms are removed, a hypothetical situation given the large difference in formation energy, δ_iso becomes 21 ppm, with a positive shift of +11 ppm due to the presence of the vacancy.

The next oxide studied is CeO₂(111). Here the formation of an O_3c vacancy, Fig. 1j, has a cost of 2.85 eV, comparable to that of r-TiO₂ and ZnO, Table 2. The adsorption energy of TMP on the Ce ion is only moderately affected and goes from −0.95 eV on stoichiometric CeO₂ to −1.08 eV on reduced CeO₂, Table 2. The Ce–P distance shows an opposite behavior and becomes slightly larger (3.131 Å) for the reduced oxide where the bonding is stronger, compared to the regular surface (3.039 Å) which exhibits a weaker interaction. The presence of the V_O3c center results in a significant change in the chemical shift of ³¹P that goes from −20 to −40 ppm. This is the same change in chemical shift observed for SnO₂(110), but is associated to a reinforcement of the adsorption energy, while in SnO₂ the additional shift of −20 ppm on the reduced surface is accompanied by a decrease of the TMP bond strength.

The last oxide considered is ZnO(101̄0), where the formation of an O_3c vacancy, Fig. 1l, costs 3.09 eV, Table 2. On the reduced surface TMP binds less strongly than on the regular surface, the same trend observed on r-TiO₂(110) and on SnO₂(110) (from −1.68 eV to −1.47 eV), but the ³¹P chemical shift remains unchanged, −41 ppm versus −40 ppm, the same effect found on zirconia. The practically identical chemical shift for regular and reduced ZnO surfaces is accompanied by slightly different Zn–P distances, 2.339 Å (regular) versus 2.396 Å (reduced), showing that a change of 0.06 Å in the distance does not result in a significant change in the ³¹P chemical shift.

These results clearly show that it is not straightforward to find a simple direct correlation between the ³¹P chemical shift and the nature of the oxide surface, stoichiometric or reduced. Despite the fact that the TMP molecule is adsorbed near the defect, the analysis of the six oxides shows a scattered behavior when an O atom is removed. On TiO₂, both anatase and rutile, when the surface is reduced there is an additional positive shift in δ_iso of about 10–20 ppm. On SnO₂ and CeO₂ the adsorption of TMP on the reduced surfaces results in an additional chemical shift of about −20 ppm, opposite to the TiO₂ case. Finally, for ZnO and ZrO₂ there is virtually no difference in the chemical shift of TMP adsorbed on the regular or on the reduced surface. Furthermore, if we add the computed adsorption energies for the reduced oxides to the correlation curves obtained for the stoichiometric surface, Fig. 2, we observe a deterioration of R². While disappointing, these results reflect the complexity of the origin of the ³¹P chemical shift in these systems.

3. Trends and correlations

In the previous sections we have presented the raw data related to the properties of TMP molecules adsorbed on the cation sites of various oxide surfaces. A first question to address is the robustness of the computed data set. To address this point, in Fig. 3a and b we have plotted the computed versus the measured ³¹P NMR chemical shifts for all the systems and sites considered. In Fig. 3a we included the δ_iso values obtained for both the stoichiometric and reduced (defective) surfaces. Of course, the number of computed points is larger than that of the experimental data, as these cannot easily distinguish the presence of steps, defects, low-coordinated sites. Fig. 3a shows a general trend but a moderate correlation with a R² factor of 0.74. In Fig. 3b we report the same analysis, but we averaged the computed values for each surface facet (regular or defective) and we obtain an overall good linear correlation with R² = 0.90. The use of an average for each facet of the computed values for stoichiometric/reduced oxides is motivated by the fact that NMR spectra provide an integrated measure of oxide nanoparticles that may contain defects and irregularities. This shows that the computational approach followed reproduces the experimental trend for the various oxides, and that the computed δ_iso values, despite the small data set, can be used to try to find non-obvious correlations with the help of machine learning algorithms (see below).

Fig. 3 Correlation between ³¹P NMR experimental and calculated chemical shifts by considering (a) δ_iso values of adsorbed P(CH₃)₃ on all surfaces (regular and defected) and (b) the average of δ_iso for each facet. Reference δ_iso: H₃PO₄ (294 ppm).

Before to attempt the use of a machine learning approach, we have plotted in Fig. S1 of the ESI† the computed ³¹P chemical shifts versus some of the fundamental properties reported in Tables 1 and 2: adsorption energy, E_ads, charge on P atom as obtained from a Bader analysis, q(P), charge on the surface cation where the TMP molecule is adsorbed, q(M), net charge on the whole TMP molecule, q(P(CH₃)₃), cation–P distance, R_M–P, Table 1. All computed values are included in these plots. The plot of the ³¹P chemical shift versus the adsorption energy, E_ads, (Fig. S1a, ESI†) shows a general trend towards more positive values for stronger adsorption bonds, but with a very poor correlation (R² = 0.41). No correlation at all exists with the charge on the P atom of TMP, Fig. S1b (ESI†), or with the charge of the metal cation, Fig. S1c (ESI†). On the contrary, some trend emerges when the Bader charge on the TMP molecule is considered, Fig. S1d (ESI†). Here all the points are along a straight line with the exception of the values associated to ZnO. Even the plot of the M–P distance versus³¹P chemical shift, Fig. S1e (ESI†), shows a complete absence of correlation. As we mentioned already, ZnO is characterized by a bivalent cation, at variance with the +IV nature of the Lewis acid sites in TiO₂, ZrO₂, SnO₂, CeO₂. Therefore, the same properties have plotted removing ZnO from the analysis. While there is a general improvement in the distribution of the data, Fig. S2a–d (ESI†), still no clear trend emerges, with the expect of the correlation between Bader charge on P(CH₃)₃ and ³¹P chemical shift, Fig. 4. Here a linear correlation is found, with a R² = 0.92. One can assume that a more positive charge on TMP corresponds to a stronger acid character of the surface (larger charge donation from the probe molecule). This is associated to a more positive value of the chemical shift. In particular, the smallest q(P(CH₃)₃), 0.05 |e|, is found in correspondence of CeO₂(110) with large negative chemical shifts, around −50 ppm. Notice that a free P(CH₃)₃ molecule with q(P(CH₃)₃) = 0 |e| has δ_iso = −60 ppm, Table 1. Intermediate values of q(P(CH₃)₃), about 0.20 |e|, are found for TiO₂ and correspond to chemical shifts of the order of −20 ppm. Finally, large positive charges, q(P(CH₃)₃) ≈ 0.4 |e|, are found in correspondence of stoichiometric or reduced SnO₂, where the ³¹P chemical shift is from about −10 to +20 ppm. According to this classification, the order of Lewis acid strength of the MO₂ oxides considered should be SnO₂ > TiO₂ > ZrO₂ > CeO₂. When we take the ³¹P chemical shift as a measure of the surface acidity, then ZnO should have a behavior intermediate between that of CeO₂ and ZrO₂ (δ_iso ≈ −40 ppm).

Fig. 4 Plot of computed ³¹P chemical shift versus the Bader charge on P(CH₃)₃ adsorbed on various M⁴⁺ sites of a-TiO₂, r-TiO₂, ZrO₂, SnO₂, and CeO₂. We also include in the plot the free P(CH₃)₃ molecule, with chemical shift −60 ppm and q(P(CH₃)₃) = 0.

Before concluding this section, we briefly discuss the use of machine learning algorithms to find non-obvious correlations. Of course, one could use other regressors, like the random forest, and it would be useful to increment the data set. However, using the methodology outlined in the computational section, we found two general descriptors that provide a sufficiently good correlation with some computed properties. In the first descriptor, the chemical shift δ_iso can be expressed as a combination of three variables, the adsorption energy of TMP, E_ads, the Bader charge on TMP, q(P(CH₃)₃), and the charge on the surface cations, q(M), δ_iso = aE_ads + bq(P(CH₃)₃) + cq(M). This descriptor correlates nicely with the computed ³¹P NMR chemical shifts, and the regression has a R² = 0.90, Fig. S3a (ESI†). Another descriptor is also function of three variables, the work function of the oxide Φ, defined as the position of the top of the valence band with respect to the vacuum level, and the same Bader charges that enter in the definition of φ₁. The corresponding equation is therefore: δ_iso = aΦ + bq(P(CH₃)₃) + cq(M). This second descriptors also shows a linear correlation with the computed δ_iso values, with an R² factor only slightly smaller than in the previous case (0.86), Fig. S3b (for the coefficients, see ESI†).

Of course, one should mention that these descriptors, while interesting, depend on quantities such as the Bader charges that can only be derived from a calculation and that are not physical observables. Therefore, it is difficult to predict the chemical shift of another system without performing a DFT calculations to obtain the Bader charges.

Conclusions

Using a DFT+D approach we have studied the adsorption properties of P(CH₃)₃ molecules on various oxide surfaces (a-TiO₂, r-TiO₂, ZrO₂, SnO₂, CeO₂ and ZnO) with the aim to evaluate the possible use of this probe molecule combined with ³¹P NMR to titrate Lewis acid sites. To this end we have considered various facets, different coverages, and both stoichiometric and reduced surfaces. In general, the computed ³¹P chemical shifts correlate nicely with the experimental values reported in the literature (Fig. 3). However, finding a universal behavior that allows one to define a scale of acid strength based on ³¹P chemical shifts turns out to be not straightforward. A linear correlation exists between adsorption energy of P(CH₃)₃ and ³¹P chemical shifts only within a given oxide polymorph (Fig. 2), in agreement with what previously found. However, if one plots the ³¹P chemical shift versus the adsorption energy for different oxides no correlation is found (Fig. 2).

When one considers a homogenous set of oxides where the cations have all the same +IV oxidation state, a-TiO₂, r-TiO₂, ZrO₂, SnO₂, and CeO₂, a linear correlation is found between the ³¹P chemical shift and the Bader charge of the P(CH₃)₃ molecule (Fig. 4), showing that this is a good descriptor of the acid strength of the surface site. In fact, a larger positive charge on the probe molecule corresponds to a stronger acid site and to a more positive chemical shift. However, Bader charges are not physical observables. Furthermore, they need to be computed, limiting the predicting power of this descriptor to establish scales of acid strength. Nevertheless, based on this result one can propose the following trend in Lewis acidity of the four oxides: SnO₂ > TiO₂ > ZrO₂ > CeO₂.

The results presented show a clear dependence of the ³¹P chemical shift on the coverage of the adsorbed probe molecules, an aspect that is rarely discussed in the experiment. On the contrary, no clear trend or effect has been found when the oxide surface is reduced by removing O atoms and forming surface O vacancies.

Finally, the use of machine learning approaches resulted in descriptors that are functions of several variables but, one more, some of these properties are related to the nature of the adsorbed probe molecules and need to be computed before to be included in the descriptor, thus limiting the potential use to predict the behavior of other systems not considered in this work.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the Italian Ministry of University and Research (MIUR) through the PRIN Project 20179337R7 and by the Australian Research Council under the Discovery Project Grant DP 200100313.

References

1. A. Marchetti, J. Chen, Z. Pang, S. Li, D. Ling, F. Deng and X. Kong, Adv. Mater., 2017, 29, 1605895.
2. H. Eckert and I. E. Wachs, J. Phys. Chem., 1989, 93, 6796–6805.
3. F. Maleki and G. Pacchioni, J. Phys. Chem. C, 2019, 123, 21629–21638.
4. F. Maleki and G. Pacchioni, J. Chem. Phys., 2019, 151, 224705.
5. E. Lima, L. C. de Ménorval, D. Tichit, M. Laspéras, P. Graffin and F. Fajula, J. Phys. Chem. B, 2003, 107, 4070–4073.
6. A. Zheng, S. B. Liu and F. Deng, Chem. Rev., 2017, 117, 12475–12531.
7. W. P. Rothwell, W. X. Siten and J. H. Lunsford, J. Am. Chem. Soc., 1984, 106, 2453–2455.
8. Y. K. Peng, Y. Hu, H. L. Chou, Y. Fu, I. F. Teixeira, L. Zhang, H. He and S. C. E. Tsang, Nat. Commun., 2017, 8, 675.
9. L. Zhang, B. Yue, X. Chen and H. He, J. Phys. Chem. C, 2018, 122, 24094–24102.
10. H. Yu, H. Fang, H. Zhang, B. Li and F. Deng, Catal. Commun., 2009, 10, 920–924.
11. T. Riemer and H. Knozinger, J. Phys. Chem., 1996, 100, 6739–6742.
12. J. Xu, A. Zheng, J. Yang, Y. Su, J. Wang, D. Zeng, M. Zhang, C. Ye and F. Deng, J. Phys. Chem. B, 2006, 110, 10662–10671.
13. K. Shimizu, T. N. Venkatraman and W. Song, Appl. Catal., A, 2002, 224, 77–87.
14. W. Zhang, Z. Lin, H. Li, F. Wanga, Y. Wena, M. Xua, Y. Wanga, X. Kea, X. Xiac, J. Chen and L. Peng, RSC Adv., 2021, 11, 25004–25009.
15. Y. K. Peng, L. Ye, J. Qu, L. Zhang, Y. Fu, I. F. Teixeira, I. J. McPherson, H. He and S. C. E. Tsang, J. Am. Chem. Soc., 2016, 138, 2225–2234.
16. Z. Tan, G. Li, H. L. Chou, Y. Li, X. Yi, A. H. Mahadi, A. Zheng, S. C. E. Tsang and Y. K. Peng, ACS Catal., 2020, 10, 4003–4011.
17. Y. Chu, Z. Yu, A. Zheng, H. Fang, H. Zhang, S. J. Huang, S. B. Liu and F. Deng, J. Phys. Chem. C, 2011, 115, 7660–7667.
18. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
19. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
20. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251.
21. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
22. V. I. Anisimov, J. Zaanen and O. K. Andersen, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 44, 943.
23. S. Dudarev, G. Botton, S. Savrasov, C. Humphreys and A. Sutton, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505.
24. F. Maleki and G. Pacchioni, J. Phys.: Condens. Matter, 2021, 33, 494001.
25. S. Tosoni and G. Pacchioni, J. Phys. Chem. C, 2020, 124, 20960–20973.
26. F. Maleki and G. Pacchioni, ACS Catal., 2021, 11, 554–567.
27. F. Maleki and G. Pacchioni, Surf. Sci., 2022, 718, 122009.
28. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
29. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
30. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
31. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
32. F. Maleki and G. Pacchioni, Top. Catal., 2020, 63, 1717–1730.
33. R. Bader, Chem. Rev., 1991, 91, 893–928.
34. J. Pickard and F. Mauri, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 245101.
35. R. Yates, C. J. Pickard and F. Mauri, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 024401.
36. L. Baltusis, J. S. Frye and G. E. Maciel, J. Am. Chem. Soc., 1987, 40(109), 40–46.
37. S. Tosoni, H.-Y. Chen and G. Pacchioni, ChemPhysChem, 2015, 16, 3642–3651.
38. C. Gionco, M. C. Paganini, E. Giamello, R. Burgess, C. Di Valentin and G. Pacchioni, Chem. Mater., 2013, 25, 2243–2253.
39. Z. Helali, A. Jedidi, O. A. Syzgantseva, M. Calatayud and C. Minot, Theor. Chem. Acc., 2017, 136, 1–16.
40. C. Franchini, M. Reticcioli, M. Setvin and U. Diebold, Nat. Rev. Mater., 2021, 6, 560–586.
41. A. Ruiz Puigdollers, S. Tosoni and G. Pacchioni, J. Phys. Chem. C, 2016, 120, 15329–15337.
42. A. Ruiz Puigdollers and G. Pacchioni, Nanoscale, 2017, 9, 6866–6876.
43. M. A. Rahman, S. Rout, J. P. Thomas, D. McGillivray and K. T. Leung, J. Am. Chem. Soc., 2016, 138, 11896–11906.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2cp03306b

---