1. Introduction
The control of crystal morphology is a complex and difficult process, which depends on the crystal internal structures and different synthetic methods, such as, reaction times, impurities incorporation, surfactants etc^{1}. The morphological importance of nanoparticles in the search of materials with high efficiency properties in optical, magnetic, catalytic, chemical, and other physical properties became increasingly relevant due to the amount of recent scientific reports^{2}^{,}^{3}.
On the other hand, computational simulation applied to materials sciences and related to researches of nanostructured systems has been more and more used aiming to predict and confirm experimental results. This tool is useful in collaborative projects between theoretical and experimental researchers.
Recently, our research group has applied theoretical methods to study the effects of surface stability on the morphological transformation of different metals and metal oxides (Ag, TiO_{2} (anatase), BaZrO_{3}, and αAg_{2}WO_{4})^{1}^{,}^{4}^{,}^{5}. By means of this method, morphologies of crystalline materials can be used as a very useful tool for the design and knowledge of synthesis of new materials. This study has a simple way of overcoming the experimental variables in order to understand the relationship between properties and structure.
In a logical sense, highenergy crystal surfaces are usually lost during crystal growth to a total surface energy minimization. Therefore, the selective exposure of high energy facets at the surface of micro and nanocrystallites is an important and challenging research topic.
In particular, the variety of phase, morphology and properties of the TiO_{2} makes this system an interesting target for this approach. TiO_{2} is a multifunctional material, which has been used in a wide variety of applications in many fields, such as in ceramics, cosmetics, medicine, food, and catalysts^{6}. TiO_{2} has three main crystalline polymorphs: rutile, anatase, and brookite^{7}. Rutile is a stable phase, whereas anatase and brookite are metastable that can be transformed into rutile when annealed^{8}.
Even though anatase TiO_{2} is considered more active than the other two main TiO_{2} polymorphs, especially when TiO_{2} is employed as a catalyst and photocatalyst^{7}^{,}^{9}, the rutile polymorph exhibits some superior physical properties, such as enhanced lightscattering properties on account of its higher refractive index, which is beneficial from the perspective of effective light harvesting. Rutile TiO_{2} is a semiconductor due to its wide bandgap, around 3.0 eV and it has a known potential to the application as a photocatalyst material^{10}.
According to Park et al.^{11}, the dyesensitized solar cells based on rutile TiO_{2} exhibit photovoltaic characteristics compared to those of conventional anatase TiO_{2}based solar cells. Furthermore, the orientation of the rutile particles may improve the charge transport of the Dyesensitized solar cell^{12}. Facts like these are important since rutile is potentially cheaper to produce.
Besides the different characteristic due to the TiO_{2} polymorphism, the anisotropy properties of this system are also frequently reported^{13}. These properties are facetdependent and, therefore, the exposed crystal facets play an important role in determining their magnitudes. Based in this argument, the synthesis of morphology with exposed facets of superior properties need an elaborate design.
Some studies reported show the importance of the influence of atomic arrangements of the (110), (100), and (001) TiO_{2} (rutile) surfaces on catalytic and photocatalytic^{7}^{,}^{14}^{}^{16}. Therefore, the modulation of rutile TiO_{2} together with its stability and low cost can be of great use in applications.
In this sense, the main objective of this paper is to show how the values of surface energies, obtained from computational simulations, can be used to obtain some possible morphologies under thermodynamic equilibrium conditions for the rutile TiO_{2} phase and also to better understand the mechanism of possible morphological transformation. This method provides an approach with both predictive and explanatory capabilities. The calculated diagrams relate the crystal growth conditions with the observed morphologies in an attempt to rationalize the morphologies obtained under different experimental conditions.
2. Computational Method and Models Systems
All computational simulations were performed by periodic density functional theory (DFT) in conjunction with Becke’s three parameter hybrid nonlocal exchange functional^{17} combined with the Lee−Yang−Parr gradientcorrected correlation functional^{18} (B3LYP) using the CRYSTAL14 software^{19}^{,}^{20}.This code uses a Gaussiantype basis set to represent crystalline orbitals as a linear combination of Bloch functions defined in terms of local functions (atomic orbitals). The titanium and oxygen atoms were described by 8651(3d)G^{21} and 8411d1G1^{22}, respectively. The functional B3LYP was chosen due to its potential for reproduce with accuracy the electronic properties of semiconductor materials^{23}, and the lattice parameters in our case gave the best agreement if compared with experimental results^{24}. Calculations with B97H, HISS, HSE06, HSESOL, LCwPBE, M062x, PBE0 and wB97x functional were conducted. However, even though all resulted in the same order of surface stability, the structural parameters were not well described.
The level of accuracy in evaluating the infinite Coulomb and HF exchange series is controlled by five parameters, T_{i}=1, 2, 3, 4 and 5, such that twoelectron contributions are neglected when the overlap between atomic functions is below 10^{Ti}. For our calculations T_{i} has been set to 8, 8, 8, 8 and 18, and the shrinking factor (PackMonkhorst and Gilat net) was set to 6.
At the beginning of the computational procedure, the optimizations of the lattice parameters and internal coordinates of rutile TiO_{2} were conducted to minimize the total energy of the structure using experimental parameters. The rutile phase belongs to the space group P42/mnm (136) with a bravais lattice (a = 4.594 Å, c = 2.959 Å and internal coordinate u = 0.305)^{24}. The optimized parameters were a = 4.571 Å, c = 2.987 Å, u = 0.305, which were in good agreement with the experimental values (variation of 0.023 and 0.028Å for a and c parameters, respectively).
From these optimized parameters the low index (100), (001), (101), (110), (111) surfaces were modeled and bidimensional slabs were built, and a new optimization of those surfaces was made as a function of the fractional coordinates. All surfaces models used here have the bottom and top planes equivalent in symmetry.
The termination, i.e., first and last layer of the surface model is the main problem of the simulation. This happens since not all terminations are stable and may not exist or quickly undergo a reconstruction. Some exposure of atoms can generate a nonzero dipole moment perpendicular to the surface, and such surfaces, when ionic are generally unstable.
The description of the selected slab models are described as follows: i) For the (100) surface, the first two atomic layers exhibit, O bridge atoms linked to two Ti atoms of the second layer, and Ti atoms exhibits coordination five, respectively. ii) The (101) surface exhibits two kinds of atoms, first layer of O atoms with coordination two and Ti atoms with coordination five. iii) The (110) surface has Ti atoms with coordination six and five, and two kinds of oxygen atoms, bridge atom, linked to two inplane Ti atoms, and inplane O atom linked to Ti atoms. In the (001) surface, there is Ti and O atoms with coordination four and two, respectively. iv) The last surface, (111) has one O atom with coordination one. The surface models and the Mulliken charges are presented in the supplementary material (see Figure S1 and Table S1).
It is important to notice that in the (110), (101), (010) and (111) surfaces the oxygen atoms are in the first layer and consequently are the ones mostly exposed. These surfaces can be considered as oxygenated surfaces.
From the thermodynamic point of view, the equilibrium shape of a crystal is determined by the free energies of various facets and can be calculated by classic Wulff construction that minimizes the total surface free energy at a fixed volume^{25}. The Wulff theorem provides a simple relationship between the surface energy, E_{surf}, of the (hkl) and has been successfully used in materials science to obtain the shapes of materials, including nanomaterials^{26}^{}^{28}. The E_{surf} is defined as the total energy per repeating cell of the slab (E_{slab}), minus the total energy of the perfect crystal per molecular unit (E_{bulk}) multiplied by the number of molecular units of the surface (n), divided by the surface area per repeating cell of the two sides of the slab, as follows:
The order of stability of these surfaces depends on the direction and the exposure of its atoms. The experimental determination of the surface energies and their stability order is not a trivial procedure. The main used technique to deduce their stability may be measuring the surface area of the single crystals by microscopy techniques. During the growth process of single crystal under thermodynamic control, facets with smaller surface energy are preferably exposed. For the construction of the Wulff crystals, it was used the VESTA^{29} program that enables the construction of the planes with the distances of interest of the crystal center. This distance is directly proportional to the energy modulation.
3. Results and discussion
The surface relaxation needs to be carefully evaluated when it is intended to develop a mechanism about the morphology changes. The surface energies to the TiO_{2} without and with the surface atoms relaxation,
Surface 



(1 1 0)  1.70  0.96 
(0 1 0)  1.32  1.01 
(1 0 1)  2.76  1.41 
(0 0 1)  3.33  1.90 
As can be seen in Table 1, an important observation must be stated, the full relaxation of the atoms significantly decreases the relaxed surface energy for (110), (101) surfaces whereas the (010) is slightly stabilized. Therefore, from the optimized results, the calculated surface stability order is (110) < (010) < (101) < (001)< (111). This stability order is in accordance with other theoretical studies^{24}.
From these surface energy values and the Wulff construction method (see Figure 1), it is possible to obtain the morphologies of
The significant difference observed in both morphologies is mainly due to the right relaxation of the (110) surface. It is noticed that the majority of the experimental particles presents a different shape comparing with the
It is known that morphologies can undergo changes because of its different synthesis environments. Modifications in the synthetic route, presence of impurities and adsorption of surfactants can act in different ways in each of the surfaces modifying the rate of growth in the respective directions. Therefore, the simulation of a real system which considers experimental factures can be very complex due to free variables that influence the surfaces, mostly unknown by experimentalists. Several methodologies can be found in the literature about theoretical simulation to explain surface influences and the modification at the equilibrium stage^{34}^{,}^{35}. Our group has used a method to evaluate the routes based on the modulation of the relationship between the surface stabilities of different faces and the areas of those parts which are exposed in the final shape^{1}^{,}^{4}^{,}^{5}^{,}^{36}.
Assuming the morphology obtained in vacuum is the starting point, it was possible to create a set of the morphological routing changes as a result of surface energy variation that takes into account the (001), (010)/(100) ratio, (101)/(011) ratio, (110) and (111) surfaces. The map of some morphologies routes of rutile TiO_{2} is showed in Figure 2. Finetuning of the desired morphologies can be achieved by controlling the values of the E_{surf} of the different surfaces.
It is worth emphasizing that there are a large number of combinations referring to all morphology possibilities in this space group, however we mainly focus in the reported experimental morphologies^{30}^{}^{33}. Experimental nanoparticles usually present morphologies found on the right side of the map showing that the surfaces (111) and (110) are more commonly stable in comparison with other surfaces. In spite of the simplicity, the explanatory power of this method has helped in the evaluation of morphological changes with the variation of the synthesis environment besides the connection between some system characteristics and the exposed surfaces^{4}^{,}^{37}. There are several situations that the morphology map can be useful. An interesting example of its usefulness is to explain the changes observed by Lai et al. (see scheme 1 of article)^{30}. In this paper, there are morphology modification by the addition of Sodium Fluoride in order to control exposed surfaces resulting in the decrease of the (110) surface in the equilibrium. In order to better study the reported modification, it is possible to perform an inverse Wulff construction choosing the exact experimental surfaces sizes to obtain the surface energies ratios. The exact morphology of the first particle obtained by Lai is presented in Figure 3a (also in Figure 2).
The first morphology observed on the left of Figure 3a has a similar size of the one obtained by Lai et al. and can be achieved by a ratio of 0.60 Jm^{2} to (110) and 1.20 Jm^{2} to (111). Therefore, it is possible to evaluate the experimental route by the stabilization of the (111) surface in relation to the (110) at a critical point at 0.40 Jm^{2} to (111) when the (110) surface is not thermodynamically stable in the equilibrium. Above the ratio value the (110) “surface” does not coexist with the (111) exposed plane. It must be clear that the assumption is about the ratio, consequently the destabilization of the (110) surface also results in the same change. The experimental particles of the cited article also show the apparition of a trait of (001) surface with a great amount of NaF addition also showing its stability. However, the figure does not allow a reliable measurement of the (001) surface area to an exact study of the ratio. Notwithstanding the above, the critical point of the ratio to the (001) apparition in relation to the (111) surface is 0.80 ((001)/(111)).Based on this last observation, Figure 3b shows the morphology transformation to the similar structure shown on the right side of Figure 3a. This structure can be seen as an octahedral structure (two inverted pyramid). The observation of this structure through experimental methods can lead to the conclusion that the formers of this morphology are the same plane. But, as seen in Figure 3b, the two pyramid structures are formed by (101) surfaces with 0.27 Jm^{2}. Once again, these observations can lead to the conclusion that distinct synthesis route can produce nanoparticles with the same morphology, but formed by distinct planes with different properties.
The above example elucidates the assumption of the main interaction of the additive with the (110) surface. However, in similar situations, the modifications are not so explicit adding even more value to this method. Besides that, even that some surface properties are well defined; the stabilization of the surface causes difference in some characteristics, mainly in the optical and electrical properties.
The method used presents a set of advantages as mentioned above. Therefore it consists in a way of how to circumvent the large number of variables of the real experiment that influences the final morphology, and also the massive amount of computational processing time that would be necessary for all calculations in every step of the map. However, these simplifications generate limitations when it is necessary to conduct more detailed electronic analyses in other structures of the map instead of the initial one. This is because the method is a classic approach after the initial structure obtained by quantum calculations to supply information about the relative surface energy and the variation needed to a morphology modification.
4. Conclusions
In this paper, rutile phase TiO_{2} was investigated by computational simulations based on DFT/B3LYP calculations and Wulff construction model on the mechanism of the morphology transformation. The low index, (100), (001), (101), (110), (111) surfaces were modeled and the calculated surface energy follows the sequence: (110)<(010)<(101)<(001)<(111). With these surface energies as the starting point, the theoretical ideal morphology and a map of some morphologies routes were constructed. The present strategy provides a way to study the morphology change routes caused by factors as different precursors, pH, capping agent and impurities. The theoretical results were used to previously evaluate reported experimental particles morphologies and the shape modification caused by differences on the synthesis methods.
The morphology map can be used as guide to support experimental works in order to analyze microscopy results and to be related with final properties of the material. Therefore this method can provide important information about the surface structure and other physical and chemical related applications.