Tuning the Electronic Structure of Anatase Through Fluorination

A highly fluorinated anatase lattice has been recently reported, providing a new class of materials whose general chemical formula is Ti1−x□xX4xO2−4x (X− = F− or OH−). To characterise the complex structural features of the material and the different F environments, we here apply a computational screening procedure. After deriving a polarisable force—field from DFT simulations, we screen in a step-wise fashion a large number of possible configurations differing in the positioning of the titanium vacancies (□) and of the fluorine atoms. At each step only 10% of the configurations are retained. At the end of the screening procedure, a configuration is selected and simulated using DFT-based molecular dynamics. This allows us to analyse the atomic structure of the material, which is strongly disordered, leading to a strong decrease (by 0.8 eV) of the band gap compared to conventional anatase.

Recently, a novel synthesis technique conducted in our laboratory [24] led to the preparation of a highly fluorinated anatase material in which fluorine replaces oxygen atoms in their lattice sites and the resulting charge deficiency is compensated by the formation of a cationic vacancy every four substitutions.The material obtained has thus the general formula Ti 1−x x X 4x O 2−4x , where the symbol is used to represent a cationic vacancy and X is either a F − ion or an OH − group (the amount of F − may vary depending on the synthesis condition).Elemental analysis and synchrotron diffraction revealed the existence of more than 20% cation vacancies, in fact the stoichiometric formula Ti 0.78 0.22 X 0.88 O 1.12 was assigned to the material.By using 19 F NMR spectroscopy it was also possible to discern three different coordination modes for the F atoms in the material: F − Ti 1 , F − Ti 2 and F − Ti 3 , highlighting the complex structural arrangement present in the material.* Corresponding Author.E-mail: mathieu.salanne@upmc.fr In this Letter, we report the results of a computational study of the fully-fluorinated material (e.g.Ti 0.78 0.22 F 0.88 O 1.12 ) performed in order to better characterize its structural features and the effects induced by the three different environments around the F atoms.The enormous number of possible structural arrangements of the vacancies and of the F atoms in the anatase structure render the problem untreatable directly by ab initio simulations.Therefore we apply a screening procedure on the possible configurations in the spirit of the emerging high-throughput techniques, see e.g.Refs.[25][26][27], by using classical Molecular Dynamics (MD).Several force-fields were previously proposed for pure TiO 2 [28][29][30][31][32][33][34].Here we derive a new polarizable force-field valid for the pure phase as well as for the fluorinated material by extracting the force-field parameter from DFT simulations, via a force and dipole fitting procedure wellestablished in our group (see for example Refs.35 and 36).The details of the force-field are discussed in the Appendix.The parameters of the force-field obtained from DFT simulations are reported in Appendix Table I and II.The classical polarizable force-field includes a damped Born-Mayer-Huggins-Fumi-Tosi term, a Coulombic term and a polarization term described by making use of the induced dipoles.
In order to generate and assess fluorinated samples starting from the pure TiO 2 anatase, we apply a screening procedure, similar for example to what done by Wilmer et al. for metal-organic frameworks (MOFs) [37].At the fixed target composition Ti 0.78 0.22 F 0.88 O 1.12 we consider samples containing: a) only F in the environment F − Ti 2 1 ; b) F in both environments F − Ti 2 1 and F − Ti 3 ; c) F in all possible environments F − Ti 1 2 , F − Ti 2 1 , and F − Ti 3 .We leave the ratio of F in the different environments free to vary at random (albeit with constraints, see below).The starting fluorinated structures are generated from the 4 × 4 × 2 pure anatase TiO 2 structure [38] (Ti 128 O 256 ) leading to a system thus arXiv:1403.5525v1[cond-mat.mtrl-sci]21 Mar 2014 composed: Ti 100 28 F 112 O 144 .These configurations are generated applying the following criteria: a) for the case with only F in the environment F − Ti 2 1 , we erase 28 Ti ions at random with the constraint that they must not be adjacent to a previously selected Ti either at the same height along [001] or on different heights.At each cationic vacancy thus created we assign a random number of localized F N i F between 1 and 6 and then for each vacancy we replace the closest N i F O with F with the constraint that the total number of F must correspond to the stoichiometric number, N F = 112.All F are therefore either found at the shortest −F distance 1.93 Å (at "equatorial" positions in the original TiO 6 octahedron), or at second shortest − F distance 1.98 Å (at "axial" positions in the original TiO 6 octahedron).b) for the case with F in both environments F−Ti 2 1 and F−Ti 3 , we create the vacancies with same procedure described in a) and then randomly substitute 112 O with 112 F. The constraint of having non-adjacent vacancies is sufficient to ensure there are no F with environment F−Ti 1 2 in the system.c) for the case with in all possible environments F − Ti 1 2 , F − Ti 2 1 and F − Ti 3 we erase 28 Ti ions at random with no constraints on the creation of adjacent vacancies and we randomly substitute 112 O with 112 F. We impose that all F and O must be attached to at least one Ti.The configurations not respecting this criterion are discarded.
For the three cases a), b) and c) a screening procedure is then initiated.The protocol of the screening procure is as follows: 1) we perform single-point energy calculations on ∼ 1.5 • 10 5 configurations generated for each case a); b) and c) as described above; we sort the configurations according to their ascending energy and we retain the ∼ 1.5 • 10 4 configurations with the lowest energy for the following step.
2) we perform 0 K geometry optimizations of the atomic positions, keeping the length of the cell vectors fixed; we sort the configuration according to their ascending energy and we retain at maximum the ∼ 1.5•10 3 configurations with the lowest energy for the following step.The configurations crashed during the run are discarded.
3) we perform 0 K cell optimizations allowing the contemporary adjustment of the atomic positions and of the lengths of the cell vectors, while keeping the box angles fixed at α = β = γ = 90 • ; we sort the configuration according to their ascending energy and we retain at maximum the ∼ 1.5 • 10 2 configurations with the lowest energy for the following step.The configurations crashed during the run are discarded.
4) we temper the configurations performing 10 ps N V T runs at finite temperatures from T 1 = 25 K to T 12 = 300 K, every ∆T = 25 K.We control the temperature through the Nosé-Hoover thermostat [39,40] with the time constant τ T = 1 ps.The 15 configurations with the lowest energy at T 12 = 300 K are retained for the following step.
5.1) we equilibrate the remaining samples for 100 ps in the N V T ensemble at T = 300 K, with τ T = 1 ps.
5.2) we equilibrate the remaining samples for 200 ps in the N P T ensemble at T = 300 K (τ T = 2 ps) and at ambient pressure.Pressure is controlled by using the Martyna barostat [41] with τ P = 4 ps.

5.
3) for the remaining samples, we perform a 1 ns production run in the N P T ensemble.The time constants of the thermostat and the barostat are set respectively to τ T = 10 ps and τ P = 10 ps.
At the end of the screening procedure we are left with 15 configurations from group a) i.e. initially containing only F − Ti 2 1 , 8 configurations from group b) i.e. initially containing only F − Ti 2 1 and F − Ti 3 and 10 configurations from group c), i.e. initially containing F in all the three environments.
In Fig. 1 we show the the relative frequency histograms built considering all the potential energies obtained at the end of each screening step from 1) to 4) for all the cases a); b) and c).Together we show the sorted energies of the configurations retained to be run in the following step.We stress here that we consider the ratios of F−Ti 1 , F−Ti 2 and F−Ti 3 as the ones calculated for the starting unrelaxed configuration.We will see in the following that their relative amount can change during the simulations.Looking at Fig. 1 we can immediately observe that the configurations prepared with only F − Ti 2 1 tend always to be the lowest in energy at all steps.This is consistent with the stabilization of the vacancies induced by the presence of the F atoms around them as observed in Ref. [24].We can also observe that upon relaxation of the structures and upon tempering of the configurations the difference in potential energy between the best configurations initially containing only F − Ti 2 1 and the best ones prepared with F in all possible environments diminishes dramatically, see Fig 1(b,d,f,h).The unrelaxed structures at 0 K show a difference in energy from 1.2 to 2.1 eV/f.u.(where we consider 128 formula units, f.u.).After the cell relaxation this ∆E reduces to 45.5-81.7 meV/f.u.Finally the tempered structures, at T = 300 K, show a ∆E 22.4-38.2meV/f.u.This difference is then comparable to the thermal energy k B T , which equals 25.9 meV at T = 300 K.
Therefore we see that the configurations initially prepared with all F environments are stabilized both by allowing the relaxation of the structure and by the effect of temperature.We underline here the importance of the thermal treatment since the lowest energy difference at 0 K after the cell optimization remains larger than 45 meV/f.u.The case initially containing only, F − Ti 2 1 and F − Ti 3 is intermediate between the other two and the difference in energy between these configurations and the ones containing initially only F−Ti 2 1 becomes very small at T = 300 K. Capturing these thermal effects in a generic, entirely ab initio based high-throughput procedure would be impossible.Generally, such studies involve static calculations only since performing molecular dynamics simulations is computationally too expensive (and, very often, useless).Nevertheless, it is interesting to test whether our selected configurations (i.e.those remaining at the end of all the screening procedure) would also have been selected if ab initio static calculations had been performed.To test this, we take their initial structures and calculate their DFT energy.Then we take the same number of random configurations from the pools a); b); c).We find that (i) the configurations given by the classical screening all have a lower DFT energy than the ones taken at random and (ii) the difference in energy between configurations containing only F − Ti 2 1 and the ones containing F − Ti 2 1 and F − Ti 3 are reproduced quantitatively in the DFT calculations while the energy difference with the configuration containing F in all environment is slightly lower in DFT.The results of this validation are shown in Appendix Fig. 5.
Next, we analyze how the initial structural arrangements correlate with the energy of the configurations.The results are shown in Fig. 2. For the configurations containing only F − Ti 2 1 , see Fig. 2(a), we see that lower energies correlate with a higher fraction of F in the "equatorial" configuration.In fact in the anatase structure each Ti (or vacancy) has 4 "equatorial" neighboring O (or F) at a shorter Ti (or )-O (or F) distance and 2 "axial" neighboring O (or F) at a ∼ 2.5% longer Ti (or )-O (or F) distance.As already indicated by the DFT calculations shown in Ref. 24, having the F closest to the vacancy stabilizes the structures.We note that the a priori probabilities of having equatorial F is 86% and the one of having axial F is 14% with the criteria mentioned above for building the pool a).In the same graph we also report the initial energies of the best configurations given at the end by the screening procedure.We observe that they are found in the lowest energy-higher fraction of equatorial F corner.
For the case in which we have F−Ti 2 1 and F−Ti 3 , see Fig. 2(b), lower energies correlate with a higher fraction of F − Ti 2 1 , again consistently with the observed stabilization of the vacancies by the F. In this case the a priori probabilities of having F − Ti 2 1 /F − Ti 3 is 65.6/34.4% according to the criteria used to build the pool b).In this case the final screening configurations are found closer to the average F speciation values rather than at the highest F − Ti 2 1 relative compositions.
Also for the case containing all the three F environments, see Fig. 2(c), lower energies correlates with the percentage of F − Ti 2 1 .Here the a priori probabilities of having F − Ti 1 2 /F − Ti 2 1 /F − Ti 3 is 11/43/46 %.The configurations surviving the screening procedure in this case are even more dispersed in terms of their F − Ti 2 1 content.
As a matter of fact the 19 F NMR analysis on the experimental sample reveals at T = 150 K relative ratios of F − Ti 1 /F − Ti 2 /F − Ti 3 of 13/70/17 % [24], indicating therefore a very special F speciation.This is not inconsistent with our findings.In fact upon relaxing the structures and increasing their temperature we also find a large increase of the number of F − Ti 2 .The F ratios after the dynamic runs is difficult to calculate exactly because the presence of disorder induces a large dependence of the results on the choice of the cutoff distance used for assigning whether a fluorine is coordinated to a titanium atom.Nevertheless we systematically observe  an increase in the number of 2-coordinated F after relaxing and heating.In the following we will consider a quite large F-Ti distance, d F−Ti = 3.0 Å for evaluating the F-Ti 1 /F-Ti 2 /F-Ti 3 speciation for the final configurations (at ambient conditions).This is also justified by the fact that the Ti-F radial distribution function does not decay to 0 between the first and the second peaks (in the cases not containing only F − Ti 2 1 initially, see Fig. 4 in the following).
In order to compare the structural properties of the material as found in the experiments with our simulations, we compare the experimental x-ray structure factor S(k) to the one we calculate from our trajectories using the Ashcroft-Langreth partial structure factors according to the formula: where x α are the relative atomic concentrations and f α (k) are the k-dependent atomic scattering.They are calculated using the analytic approximation: where the coefficients a α,i , b α,i and c α are taken from Ref. [42] for O 2− and from Ref. [43] for Ti % (0.9/87.we see that the classical potential used might have the drawback to exaggerate the percentage of 2-coordinated F atoms with respect to the experiment where the F speciation is found to be 13/70/17 % (albeit this was found at T = 150 K).
Nonetheless we can observe in Fig. 3 that the agreement between the experimental and the calculated structure factor is quite good.The peak slight misalignment is due to an enlargement of the simulation box with the respect to the experiment.In fact one average we found a 2.6% increase of the lattice parameters for pure anatase, to be compared with the 0.5% increase in DFT calculations.It is clear that the intensity of the peaks for the configurations initially prepared with only F − Ti 2 is overestimated but in the other two cases the intensity of the peaks is almost in quantitative agreement with the experiment.The higher level of noise in the simulations results is due to the relatively small size of the box compared to the experiment.The results described for the S(k) validate the study of the structure by our model.Consequently, in Fig. 4 we show the partial radial distribution functions at ambient conditions g αβ (r) for the simulated pure anatase TiO 2 and for the configurations labelled A in Fig. 3.The effect of the disorder introduced by the vacancies is immediately evident looking at the g Ti−Ti (r) and at the g O−O (r).The g Ti−O (r) structure seems to be conserved to a large extent at least for the first two shells, although the effect of disorder is evident when looking for example at the region in between the first two peaks.When comparing the configuration initially having only F − Ti 2 to the other two cases it is evident that the first conserves a greater degree of order with respect to the others.In fact all its g αβ (r) present peaks and shells which are more defined.As we have seen in Fig. 3 however the experimental structure should be more similar to the one of the configurations prepared having initially also F−Ti 1 or F−Ti 3 .For those cases we can observe the similarity of the g Ti−O (r) and g Ti−F (r), confirming the substitutional nature of the F on the O sites also at ambient temperature.However we also note that the Ti − O first neighbor distance become smaller in the fluorinated samples with the first peak of the g Ti−O (r) centered at 1.81 Å, to be compared with 1.94 Å for g Ti−F (r).Connected to this, the first peak of the g O−O (r) is at a larger distance ( 2.9 Å) with respect to the first peak of the g F−F (r) ( 2.6 Å).We also observe that the g O−O (r) is fairly similar to the g F−F (r).Summarizing, we find that fluorine is on average farther from Ti than O and and closer to O than other oxygen ions, while the distance of fluorine from other F ions is found to be similar to the distance between F and O.
In conclusion, in order to characterize fluorinated anatase Ti 0.78 0.22 F 0.88 O 1.12 as produced in experiments [24], we develop a classical polarizable force-field from DFT simulations.With the derived force-field we initiate a screening procedure on three large pools of configurations containing a) only F − Ti 2 1 ; b) F − Ti 2 1 and F−Ti 3 ; c) F−Ti 1 2 , F−Ti 2 1 and F−Ti 3 .We run N P T simulations at ambient conditions on the best configurations left after the screening procedure.The comparison of structure factor calculated from the our simulations to the one measured in experiments reveals a very good agreement of the structural features.We find that upon relaxation and heating the number of 2-coordinated F increases, consistently with experimental results.The study of the partial radial distribution functions shows that the average first neighbor distance is shorter for Ti-O than Ti-F, the reverse is true for O-O and O-F.Our force-field combined to the screening procedure allowed us (i) to select the best configurations starting from a very large pool of possible configurations; (ii) to reproduce the experimental structure to a good extent, (iii) to study details of the partial atomic structure unaccessible in experiments.Furthermore, this procedure appears promising in general for the study of materials presenting complex structural arrangements.
Appendix Figure 5 shows the validation of the screening procedure performed by DFT calculations as described in the text.The derivation of the polarizable classical force field is also described in the Appendix.The parameters of the potential are reported in Appendix Tables 1 and 2.

Classical model and derivation
We describe the interaction potential between the ions by a classical polarizable force-field whose parameters we derive from ab initio DFT simulations.Details on how to extract the parameters of the classical potential from DFT simulations are reported in earlier works [35,36].
Here in the following we give the details of the DFT simulations performed and we illustrate the details of the chosen model for the ions.We then present the parameters derived by fitting the forces and the dipoles obtained by the classical model to the ones obtained in DFT.

DFT calculations
The initial structures are obtained form the pure anatase TiO 2 3 × 3 × 1 structure (36 TiO 2 units, total 108 atoms) creating vacancies and substituting F to O according to the stoichiometric formula Ti 1−x x F 4x O 2−4x .The DFT simulations are performed using the GGA PBE [44] functional to describe the exchange-correlations interactions.We use the QUICK-STEP algorithm [45] for the evaluation of the forces as implemented in CP2K [46].We use the DZVP-MOLOPT-SR-GTH basis set [46,47] and we employ the Goedecker-Teter-Hutter pseudopotentials [48][49][50].The valence orbitals explicitly represented are 3s 2 3p 6 3d 2 4s 2 for Ti atoms, 2s 2 2p 4 for O atoms and 2s 2 2p 5 for F atoms.The plane-wave cutoff is set to 400 Ry.We add dispersive interactions using the DFT-D3 correction [51] and a cutoff radius of 30 Å.We apply isotropic boundary conditions.We make use of the maximally localized Wannier functions (MLWFs) [52] to localize the Kohn and Sham orbitals.

Classical force-field
Here we present the classical polarizable model employed for the description of the atomic species and we show the parameters that we derived from the DFT simulations.The repulsive and dispersive terms of the interactions are taken into account using the the Born-Mayer-Huggins-Fumi-Tosi-damped (BMHFTD) form of the interaction potential: The damping functions are Tang-Toennies functions of the form The parameters extracted from DFT calculations are reported in Table I.
For the Coulombic part of the interaction potential the formal charges for the ionic species are used, −2e for O ions, −e for F ions and +4e for Ti ions.The manybody electrostatic interactions are described by the induced dipoles µ i , obtained at each MD step minimizing the polarization energy q i µ j α g ij (r ij ) − q j µ i α g ji (r ij ) where the Einstein summation convention is assumed, α i is the atomic polarizability and T are the multipole interaction tensors.The damping function g ij (r ij ) is of the Tang-Toennies form The parameters extracted for the polarization part of the potential are reported in Table II.The classical MD simulations are performed using the software CP2K (step 1, 2 and 3 of the screening procedure, see main text) and the in-house simulation software PIMAIM (step 4 and 5 of the screening procedure, see main text), the short-range interactions are cut off at half the norm of the shortest box vector (or less in N P T runs).Isotropic periodic boundary condition are applied.The dipolar Ewald summation method is employed to deal with the electrostatics [53,54].The time step for the integration of the equations of motion in the finite temperature runs is 1 fs.

FIG. 1 .
FIG.1.The panels on the left show the relative frequency histograms of the energies calculated for all the configurations tested at each step, while the panels on the right show the sorted energies of the configurations retained after each step.The energies shown are: (a,b) at 0 K for the starting unrelaxed structures; (c,d) at 0 K after the optimization of the atomic positions; (e,f) at 0 K after the optimization of the atomic positions and cell vector lengths; (g,h) at 300 K after the tempering from 25 to 300 K.All steps are performed for configurations initially containing F − Ti2 1 only (black), F − Ti2 1 and F − Ti3 (red) and F − Ti1 2, F − Ti2 1 and F − Ti3 (blue).

FIG. 2 .
FIG. 2. For the initial configurations at 0 K (a) fraction of "equatorial" F atoms (at the shortest F − distance) vs. energy of the configuration for the case containing F − Ti2 1 only; (b) fraction of F − Ti2 1 vs. energy of the configuration for the case containing F − Ti2 1 and F − Ti3; (c) fraction of F − Ti2 1 vs. energy of the configuration for the case containing F − Ti1 2, F − Ti2 1 and F − Ti3.The points are represented as small brown circles.For each different value of the fraction of "equatorial" F (a) or F − Ti2 1 (b, c) we calculate the mean (black circles) and the standard deviation (black bars) of the corresponding energies.We also report the values assumed at this stage by the configurations run at the final screening step (red squares).

FIG. 3 .
FIG.3.Comparison of the structure factor S(k) at ambient conditions measured in experiments (black solid line) and calculated as described in the text from N P T runs of two configurations given by the screening procedure that we label A (red solid line) and B (green dashed line) for each case: initial F − Ti2 1 only, F − Ti2 1 and F − Ti3 and F − Ti1 2, F − Ti2 1 and F − Ti3.

FIG. 4 .
FIG.4.Partial radial distribution functions g αβ (r) at ambient conditions calculated from N P T simulations of pure TiO2 anatase (black solid line) and from simulations of one of the configurations given by the screening procedure (configuration labelled A of Fig.3) for the cases: initial F − Ti2 1 only (red solid line), F − Ti2 1 and F − Ti3 (blue dashed line) and F − Ti1 2, F − Ti2 1 and F − Ti3 (green dot-dashed line).
FIG.5.For the configurations left at the end of the screening procedure, we take their initial unrelaxed structures and calculate their DFT energy (circles).Those are compared with the DFT energies of the same number of configurations taken at random (squares) from the initial pools of configuration with F − Ti2 1 only (black); F − Ti2 1 and F − Ti3 (red) and F − Ti1 2, F − Ti2 1 and F − Ti3 (blue).

TABLE I .
Parameters of the BMHFTD potential extracted from DFT simulations.

TABLE II .
Parameters of the polarization part of the interaction potential extracted from DFT simulations.