2.1. Computational details

The Quantum Espresso (QE) package ^24,25 was used to perform all the spin-polarized DFT calculations in the generalized-gradient approximation (GGA) with the Perdew-Burke and Ernzerhof (PBE) ²⁶ exchange and correlation functional and the van der Waals dispersion forces (vdW) ²⁷. A plane wave expansion was used to represent the electronic states of the valence electrons and Vanderbilt ultra-soft pseudopotentials were used to represent the interactions of the valence electrons with the ions. The k-point grid was generated through the Monkhorst & Pack method ²⁸. All the crystalline structures were generated using the XCrysDen program ²⁹.

2.2. Convergence study of the simulation parameters

In ideal DFT periodic calculations for crystalline solids, it is desirable to use the largest possible k-points grid and also high values for the kinetic energy cutoff for the wave-functions (E _cut) and the kinetic energy cutoff for the charge density and potential (E _cutrho), with the aim of representing the system in the most realistic way possible. Nevertheless, given the time constraints of simulation research works and the computational limitations of the current numerical processors, a balance should be found between the quality of the results and the efficiency of the calculations. This is the motivation for a parameter convergence study, which seeks to reduce the relative error of the results by calculating the energy differences between two crystalline phases of the solid as a function of a specific parameter.

In this work, the energy differences between the tetragonal P4/nmm and the simple cubic phases of WO₃ were calculated. The E _cut value was started at 5 Ry and increased by 5 Ry at each step up to a maximum value of 120 Ry. E _cutrho values of 8, 10, and 12 times the value of E _cut were tried. The number of k-points for the x axis (k _x ) was varied between 2 and 10, with successive increments of 2, and the number of k-points for the y and z axes (k _y and k _z ) were obtained from k _x , depending on the conventional crystal size.

If the lattice parameters of a unit cell obey the constraints a=b<c where c=2a, then the number of k-points for the x, y, and z axes can be written as k _x , k _x , and k _x /2, respectively. Otherwise, if the parameters obey the constraints a=c<b where b=2a, then the number of k-points can be written as k _x , k _x /2, and k _x , respectively.

The optimal values of the parameters were found by means of an iterative process in which the parameter values that were obtained at the end of the first convergence test were used as the starting values for the second convergence test, and the values obtained from that second test were used as the starting values for the third test, and so on until the parameters reached stable values.

The iterative process used to obtain the optimal values of the parameters is comprised of the following steps:

Step 1. The value of E _cut is varied until convergence is achieved while using the values of E _cutrho and k _x that were obtained in the previous cycle. For the first cycle, the starting values of E _cutrho and k _x are 10 and 2, respectively, and k _y and k _z are defined by following the rules explained above.

Step 2. The value of E _cutrho is varied until convergence is achieved by keeping unchanged the value of k _x that was obtained from the previous cycle and using the value of E _cut obtained from step 1.

Step 3.The value of k _x is varied until convergence is achieved while using the values of E _cut and E _cutrho obtained from steps 1 and 2, respectively.

Step 4. The values of E _cut , E _cutrho , and k _x obtained from steps 1 to 3 are compared with the values obtained at the end of the previous cycle. If the two sets of values are in disagreement, another iterative cycle (steps 1 to 3) is started. If the two sets of values are in agreement, the iterative process is terminated.

A similar parameter optimization method has previously been implemented by one of the authors of the present work ³⁰. For WO₃ (tetragonal P4/nmm), the optimization of the parameters required two iterative cycles and the corresponding results are given in Table 1.

The optimal values of the parameters are as follows: E _cut= 50 Ry, E _cutrho = 10 E _cut y n _k= 8. These values were used in all subsequent calculations for the tetragonal phase of WO₃.

2.3. Surface energy

The surface energy, E _{surf ,} was calculated from the following expression ^31-35:

where A is the area of the slab perpendicular to the z axis, E _slab is the total energy of the slab, E _bulk is the energy of the structure as a bulk, n _slab is the number of formula units of the slab and n _bulk is the number of formula units in the bulk. When the structure is stoichiometric, there is no loss or gain of energy by removing or adding atoms to the slab; in consequence, the summation term ±∑i=1M∆NiEi is equal to zero and Eq. (1) simplifies to the first two terms. If the structure is not stoichiometric, it is necessary to take the summation term into account, where ∆N _i stands for the number of atoms being removed from or added to the slab and E _i is calculated from the following expression:

Where Eix and ∆μix are the total energy and the chemical potential change of molecular species i _x , and x is the number of atoms of element i in species i _x.

3.1. Tetragonal (P4/nmm) WO ₃ bulk

Based on the optimal values of the parameters, a full relaxation of the tetragonal P4/nmm WO₃ unit cell was performed. All atomic positions and lattice parameters of WO₃ were relaxed by using the variable-cell command of the QE package. The starting values of the lattice parameters and atomic coordinates were taken from the study by Locherer et al. ²¹.

The WO₃ structure has three different W-O bonds, designated as W-O1, W-O2 y W-O3. Also, two angles, designated as θ and (, are required to describe the structure’s geometry. The results obtained in this work agree with those reported in the open literature (see Table 2 and Figure 1).

Figure 1 Tetragonal P4/nmm bulk. The oxygen atoms are in red, the tungsten atoms are in gray. 

3.2. Density of states and band gap of the P4/nmm WO ₃ bulk

A density of states (DOS) calculation for the bulk was performed by taking into account the geometry described in Section 3.1. Figure 2 shows that the 2p states of the oxygen atoms are predominant in the valence band and the 5d states of the tungsten atoms are predominant in the conduction band.

Figure 2 Density of states (DOS) of tetragonal P4/nmm WO₃. 

The band gap value obtained in this work is similar to the values reported in other theoretical studies (see Table 3). As far as the authors’ knowledge is concerned, no experimental values have been reported for the band gap of tetragonal P4/nmm WO₃. However, Wang et al.³⁶ suggested that the band gap of the tetragonal phase should be similar to that of the cubic phase because the tetragonal phase is just a deformation of the cubic phase. Indeed, several computational methodologies applied to these phases yield similar results ³⁶. Even so, the band gap of 2.6 eV reported by Hjelm et al. ³⁷ for the cubic phase is quite different from the values listed in Table 3 for the tetragonal phase. Such discrepancies are common in DFT studies for solids and are related to the fundamental aspects of the theory at the physical interpretation of the bands and orbitals ³⁸. A possible way of elucidating such discrepancies is by implementing a specific hybrid functional to emulate the band gap ³⁹. In this work, the widely used functional known as PBE was chosen because it has been already optimized using ultrasoft pseudopotentials, which ensures the reliability of the geometry and energy calculations.

3.3. Tetragonal P4/nmm (001) surfaces

The plane (001) is most frequently mentioned in the work by Salazar ²⁰. The (001) surface is equivalent to the WO₃ bulk with a vacuum space in the z axis.

In Figure 1, there are two protruding oxygen atoms; this could be a problem during the adsorption of electronegative molecules such as water and ethylene (the reactants for the production of ethanol) due to the electrostatic repulsion between negative point charges ^41,42. Also, it is convenient to work with structures without dipoles because these can be troublesome for the calculations ⁴³.

From these considerations, two surface models are proposed. Model I has no protruding oxygen atoms on the surface but it is not stoichiometric with regard to the WO₃ unit. Model II has the same number of oxygen atoms at both sides of the surface and it is stoichiometric with regard to the WO₃ unit. The layers are defined for each model as shown in Figure 3.

Figure 3 Models proposed for the (001) surface of the tetragonal P4/nmm WO₃ a) Model I. Non-stoichiometric and non-dipole structure b) Model II. Stoichiometric and non-dipole structure. 

A proposal similar to model II has been reported for the monoclinic WO₃ (001) surface by taking into account superficial reconstructions of experimental works ⁴⁴.

3.4. Surface energy expressions for the tetragonal P4/nmm WO ₃ (001) surface

Suitable values for both the number of layers and the vacuum distance between the slabs have to be chosen in order to represent the nanostructures with a modest computational cost and achieve an appropriate model of the surface.

By following a procedure similar to the parameter convergence tests for the bulk described in Section 2.2, the optimum number of layers and the vacuum distance between the slabs were found by searching for the values of those parameters that turn out to be associated to consecutive calculations with a surface energy difference less than a prescribed tolerance. In this work, we used a tolerance of 0.5 meV/ Ų, which looks much lower than the one used in similar works ^44,45.

In contrast to the parameter convergence tests for the bulk (in which total energies are involved), the convergence of the surface model is determined from the computation of surface energies.

As model I can be obtained from model II by removing two oxygen atoms (one atom at each end of the slab, see Figure 3), the summation term for model I is not null and the energy for model I should be larger than the energy for model II. Therefore, for model I, Eq.(1) takes the following form:

For a convenient calculation of E ₀, the easiest choice for the molecular species i _x in Eq.(2) is the molecular oxygen (triplet) O₂. In that way, Eq.(2) takes the following form:

Where the chemical potential change ∆μ02 is given by the expression

In Eq.(5), R is the universal gas constant, T is the absolute temperature, a02 is the activity of O₂, and a02o is the standard chemical potential of the oxygen molecule and is given by the following expression ⁴⁶:

Where n stands for the number of oxygen molecules, and GO20, HO20 and SO20 are the standard Gibbs free energy, enthalpy and entropy of formation of the oxygen molecule, respectively, which by definition are set to zero, because the triplet oxygen is the most stable structure relative to the oxygen atom. Besides, as the system is calculated purely from DFT, the absolute temperature is equal to 0 K. In this way, from Eqs. (5) and (6), it follows that ∆μO2=0 .

Taking all of the above into account and the fact that ∆NO=2 because two oxygen atoms are removed when the molecular species is O₂, the surface energy of model I is given by

The summation term in Eq.(1) vanishes for model II because this is stoichiometric. Therefore, for model II, Eq.(1) takes the following form:

3.5. Layer and vacuum convergence for the tetragonal P4/nmm (001) WO ₃ surface

The procedure used to obtain the number of layers and the vacuum distance between the slabs can be summarized as follows. Firstly, the number of layers was converged using an arbitrary vacuum distance (16 Å), i.e., the number of layers was successively increased until the surface energy E _surf did not change more than a prescribed tolerance of 0.5 meV/Ų. By using the resulting number of layers, the vacuum distance was changed until the surface energy did not change more than the same prescribed tolerance. This procedure differs from the one previously used for the parameter convergent tests because it is not iterative. These kinds of calculations are computationally more costly, and the assumption is made that changes will not occur in a second iteration because as shown in Figure 4b, no significant changes occur by varying the vacuum distance.

Figure 4 a) Surface energy as a function of the number of layers and b) surface energy as a function of the vacuum distance between slabs for the tetragonal P4/nmm (001) WO₃ surface. 

Both models of WO₃ can be represented accurately by four layers and at least 10 Å of vacuum. However, in order to avoid the self-interaction between periodic slabs, a final value of 14 Å was chosen for all subsequent calculations.

3.6. Relaxation of the tetragonal P4/nmm WO ₃ (001) surface

A full relaxation of only the atoms inside the slabs was performed for both models by using the values of the number of layers and vacuum distance found previously. This relaxation process is similar to the relaxation of the bulk (vc-relax), except for the avoidance of supercell movements.

As shown in Figure 5, the structural changes of both models due to such a relaxation are quite drastic. The surface model I turned into a tungsten coplanar structure with respect to the z axis (W1-W2, W3-W4, etc.). On the other hand, in the relaxed surface model II each column with a WO end is pulling along its length.

With regard to the bond distances (Table 4), surface model I has symmetric W-O bonds along the z axis, and surface model II simultaneously exhibits very short and long bonds. That model II changes more than model I by relaxation follows from the fact that the sum of the absolute values of the bond length differences before and after relaxation of model II is about three times the corresponding sum for model I.

Figure 5 Views of the tetragonal P4/nmm WO3 (001) surface: (a, b) and (e, f) model I and model II before optimzation, respectively. (c, d) and (g, h) model I and model II after optimization, respectively.  

3.7. Surface energy of the tetragonal P4/nmm WO ₃ (001) surface

The relaxation or geometry optimization of the surface models not only yields a change of structure but also a change of surface energy, which decreases for both models, and the larger decrease occurs for model II.

The decrease in surface energy is related with the stability of the slab and the rearrangement of the superficial atoms, mainly. This result is as expected: more stability in exchange with geometrical distortions.

The relaxed surface model II turns out to be energetically more favorable than the relaxed surface model I because all the tungsten atoms in model I are penta-coordinated (instead of hexa-coordinated as in the bulk), and in consequence, surface model I has a larger energy.

Model II turns out to be more suitable to represent the WO₃ catalyst because of its lower surface energy and correspondingly greater stability. This model has a surface energy comparable to that of the same plane of the monoclinic phase of WO₃ (Table 5). Despite being different crystalline phases, the monoclinic (001) surface is similar to surface model II of the tetragonal phase. To date, no other works have been published to report the energy of a WO₃ surface without protuberant oxygen atoms. This could be due to the additional difficulty involved in the implementation of the mathematical expression, Eq. (1), for the calculation of the surface energy of non-stoichiometric slabs.

3.8. Projected density of states and band gap of the tetragonal P4/nmm WO ₃ (001) surface

Figure 6 shows the projected density of states (PDOS) of the tetragonal WO₃ surfaces. As in the bulk, in surface models I and II the presence of 2p states for the oxygen atoms in the valence band and 5d states for the tungsten atoms in the conduction band is remarkable; however, in the surface models the additional contribution of 5d states of tungsten in the conduction band and 2p states of oxygen in the valence band can be noticed.

Figure 6 a). Projected density of states of the tetragonal surfaces (a) model I and (b) model II.  

In both cases the band gap is practically zero, which indicates a migration of the semi-conductor character toward a semimetal character ⁴⁷. For model I, a shift occurs for the PDOS of electrons with the same energy but different spin.

The PDOS reported in similar studies show a band gap which tends to be small (but not vanishing) for the WO₃ (001) surfaces (Table 6); this feature seems to be peculiar for the tetragonal P4/nmm WO₃⁴⁸.