Structure and Stability of Small TiO2 Nanoparticles

The results from the two techniques are similar, although the methods based on ..... help to remove bad features (for example, by separating overlappi...
15 downloads 0 Views 312KB Size
Download PDF
J. Phys. Chem. B 2005, 109, 15741-15748

15741

Structure and Stability of Small TiO2 Nanoparticles S. Hamad, C. R. A. Catlow, and S. M. Woodley DaVy Faraday Research Laboratory, The Royal Institution of Great Britain, 21 Albemarle Street, London W1S 4BS, United Kingdom

S. Lago and J. A. Mejı´as* Departamento de Ciencias Ambientales. UniVersidad Pablo de OlaVide, Carretera de Utrera, Km 1. 41089 SeVilla, Spain ReceiVed: April 27, 2005; In Final Form: June 13, 2005

The effect of the nanostructure on the photochemistry of TiO2 is an active field of research owing to its applications in photocatalysis and photovoltaics. Despite this interest, little is known of the structure of small particles of this oxide with sizes at the nanometer length scale. Here we present a computational study that locates the global minima in the potential energy surface of TinO2n clusters with n ) 1-15. The search procedure does not refer to any of the known TiO2 polymorphs, and is based on a novel combination of simulated annealing and Monte Carlo basin hopping simulations, together with genetic algorithm techniques, with the energy calculated by means of an interatomic potential. The application of several different methods increases our confidence of having located the global minimum. The stable structures are then refined by means of density functional theory calculations. The results from the two techniques are similar, although the methods based on interatomic potentials are unable to describe some subtle effects. The agreement is especially good for the larger particles, with n ) 9-15. For these sizes the structures are compact, with a preference for a central octahedron and a surrounding layer of 4- and 5-fold coordinated Ti atoms, although there seems to be some energy penalty for particles containing the 5-fold coordinated metal atoms with square base pyramid geometry and dangling TidO bonds. The novel structures reported provide the basis for further computational studies of the effect of nanostructure on adsorption, photochemistry, and nucleation of this material.

1. Introduction Since the earliest publications (see, for example, refs 1-3), the study of the catalytic properties of semiconductor oxides has evolved into a major field of contemporary research, owing to the great potential of these materials for converting light into chemical energy that is usable, for example, for the removal of pollutants and the production of hydrogen and electricity. These applications may contribute soon to solving many environmental and energy related issues.4 There are several materials that have the potential to act as photocatalysts for a wide range of applications including TiO2, WO3, SrTiO3, R-Fe2O3, ZnO, and ZnS. Among them, TiO2 is the most widely used material since it can show high activity while being chemically and biologically inert as well as inexpensive.5 Moreover, it is thoroughly documented that the method of preparation influences strongly the catalytic properties of this material, rendering it completely inactive or, contrarily, highly efficient.6 The control of the crystal polymorph, dopants, surface defects, supported metal clusters, and interaction with a support are topics of intensive research.7,8 Additionally, the control of TiO2 at the nanometer length scale and the study of its effect on photocatalysis and related processes are becoming increasingly topical.9-19 Thus, identification of the structure of the most stable small TiO2 aggregates is needed since it is the starting point for further studies on quantum size effects on photophysics, substrate-particle or particle-particle binding, and charge transfer as well as for understanding the nucleation of larger particles.

Current knowledge on the structure of small TiO2 nanoparticles is however incomplete. It appears that, although the most stable polymorph for the bulk material at room temperature and atmospheric pressure is rutile, TiO2 nanoparticles adopt the anatase structure, most probably due to its lower surface energy. EXAFS, XANES, UV-vis, and Raman studies of TiO2 nanoparticles show the increasing presence of distorted octahedral as well as 5- and 4-fold coordinated Ti atoms as the particle size decreases.20-24 There has also been a computational approach to the structure of small TiO2 particles based on the anatase crystal structure and clusters that fulfill some general rules such as lack of dipole moment and appropriate coordination.25 There are also modeling studies of the morphology of TiO2 nanocrystals26-29 which are, however, based on the stability and surface energies of the different bulk polymorphs. Other calculations have focused on the interaction of water and catechol with small TiO2 particles that also have anatase-like structures.30 Computational studies on the structure and adsorption on rutile and anatase surfaces have been reported as well, but none of them considers other type of structures when the particles become small.31-36 However, it is quite possible that nanoparticles consisting of a few TiO2 units may have structures other than those deriving from the bulk material. Here we focus on finding the most stable structures of small TiO2 particles, with sizes 1 makes the search of GM with a standard (deterministic) minimization algorithm difficult because of the need to start the search from many different geometries, each converging to a local minimum; second, the energy as a function of geometry must be evaluated accurately and inexpensively. Here we put forward a combination of methods that allow us to tackle these problems. For the search of GM we implement several algorithms, namely, a hybrid of simulated annealing-Monte Carlo basin hopping (SA-MCBH), a standard genetic algorithm (GA), and a genetic algorithm hybrid. During the search for the GM the energy of formation is defined using interatomic potentials (IP), after which we employ, and compare with, more accurate density functional theory (DFT) calculations for the more plausible configurations found. Our approach is used for TinO2n with n ) 1-15 particles. We find that, up to the largest size studied here, Ti15O30 with diameter ≈ 1 nm, the structures do not resemble an anataselike structure. Unlike small ZnS clusters45 that have ring or bubble-like structures, the small TiO2 particles are more compact, with Ti atoms having 4- and 5-fold coordination. The larger particles have a central octahedral TiO6 unit enclosed by lower coordination sites. The comparison between the IP approach used here and DFT shows a general good agreement between both approaches, although there are some differences between the structures calculated with the IP and those calculated with DFT for some sizes. The stability of the clusters is predicted correctly by the IP as compared to DFT calculations, although some improvements in the IP for these small clusters are suggested. 2. Method and Computational Details 2.1. Interatomic Potential for GM Search. For the search of the global minima (GM), we have chosen the interatomic potential of Matsui and Akaogi.37 As shown by Oliver et al.,38 this force field is able to reproduce the structure of TiO2 polymorphs. Both tetrahedral and octahedral Ti coordination states exist in these polymorphs. More recently, Bandura and Kubicki39 used this potential as a starting point for the simulation of hydrated TiO2 surfaces. Therefore, the force field is a reasonable choice for the search of local minima in small particles that may have octahedral as well as lower coordination sites. The potential form used by Matsui and Akaogi consists of a Buckingham potential plus electrostatic, charge-charge, interactions. The total energy of the nanoparticle is given by

E)

1

∑ 2 i*j

{

qiqj rij

+ Aije-rij/Fij -

Cij

}

rij6

(1)

The charges are Ti ) 2.169 |e| and O ) -1.098 |e|, and the values of A, F, and C are given in Table 1. 2.2. Locating GM Structures for Clusters. A review of the development and application of global optimization techniques to structural prediction of clusters can be found elsewhere.40-42

Typically a genetic algorithm (GA) or simulated annealing (SA) approach43,44 is implemented to generate a subset of configurations, which hopefully does not contain too many unwanted structures because the energy for each one will need to be computed. Among these structures, one is the desired, or more stable, or at least a configuration that is easily refined to the GM structure. Because an exhaustive, or complete, search of the configurational space is not performed, there is always uncertainty over whether the GM has been located. As n increases, the search space and the difficulty of finding the target configuration increase rapidly and so does the uncertainty of whether the most stable structure found is the GM. Even when a particular method has been found or developed that is able successfully to generate what is typically the desired GM structure, occasionally another low energy structure is mistakenly labeled as the GM structure.45 To increase the confidence in whether we have truly found the most stable structure, or GM, we have implemented three different search algorithms: a SA hybrid, a GA, and a GA hybrid approach. Further details of all three approaches are given below. The energy, as defined in eq 1, is minimized within each of the three methods. For any particular cluster size, the SA hybrid approach requires much more CPU time than our other two approaches. However, the SA hybrid is only used once for each size of cluster to generate the “first proposed” GM structure and, more importantly, the energy of this configuration becomes the target energy that the other two approaches must at least find. Our GA and GA hybrid are then used to check whether all the important areas of the energy landscape have been probed. That is, the latter approaches are continuously implemented with different random starting configurations and optimization parameters until either the same energy structure or a lower energy structure is found. Thus, the target structure or new lower energy structure becomes our predicted GM structure. It is important to note that we report the GM structures rather than a comparison of the three search algorithms. Thus, we do not attempt to tune each method so that a minimum number of candidate structures are evaluated before the GM structure is found. Moreover, as described below, our GA approach can be tailored for a particular target structure, or indeed so that the target structure is excluded from the space search. This is achieved by predefining a set of grid points where ions can be sited during the search. Thus we utilize the strengths of each approach during the different stages of our search. 2.2.1. Approach 1: Hybrid Simulated Annealing-Monte Carlo Basin Hopping (SA-MCBH). One of the first GM methods to be implemented was that of simulated annealing (SA).44 This method mimics the natural process by which a molten system gradually decreases its energy as the temperature decreases, and eventually reaches the GM configurationsthe most stable crystalline phase. It is a direct analogy of the physical process of cooling. The first step, when performing a SA, is to simulate the evolution of the system at very high temperature, using, for example, molecular dynamics (MD). At that temperature, the kinetic energy is sufficiently high to allow the system to overcome energy barriers between local minima very easily. The temperature is then lowered by a small amount, and the evolution of the system is simulated again with MD at the new temperature. The fall in temperature permits a lower potential energy configuration to be obtained and decreases the probability of jumps between local minima. If this annealing process is continued and carried out slowly enough, the final structure is likely to be the GM.

Structure and Stability of TiO2 Nanoparticles

J. Phys. Chem. B, Vol. 109, No. 33, 2005 15743

Figure 1. Schematic representation of one-dimensional potential energy ˜ C, dotted line). hypersurface (EC) and its transformed hypersurface (E

There is another powerful GM technique, Monte Carlo basin hopping (MCBH).46 MCBH is based on a deformation of the potential energy hypersurface. The transformed potential energy is defined by

E ˜ C(X) ) min{EC(X)}

(2)

where X represents the vector of nuclear coordinates and “min” means that a local energy minimization is performed starting from the configuration X. The transformed potential energy surface consists of a series of flat steps; the height of each step is the energy value of the local minimum for that region. See Figure 1. Monte Carlo is then used to search this transformed landscape; random changes to the current configuration are applied and accepted or rejected according to the Metropolis criteria. Note that we will always assume that the temperature parameter in the Metropolis criteria44 is fixed. The advantage of basin hopping is that the barriers have been lowered considerably and, in some cases, transitions are now barrierless. In Figure 1 we can see how transitions between configurations in region A and region B are much easier with the transformed potential energy surface. Any Monte Carlo move that changes the configuration from well A to B is automatically accepted. Based on the Metropolis criteria, moves from B to A are also more probable. As a result, a much broader sampling of the configuration space is achieved, largely improving the efficiency of the Monte Carlo sampling. We combine the two techniques, as described above, for the hybrid SA-MCBH method, which works as follows: First we run a 100 ps MD simulation of the system at high temperature (3000 K), using the DL•POLY code.47 Then the lowest energy configuration of the 3000 K simulation is chosen as the starting configuration for the next search using MCBH. Now the best structure obtained is likely to be the GM; however, there is still a possibility that it is only a local minimum. Therefore, this configuration is used to seed the next MD run at 2990 K. Again the MD will act as a generator of sensible structures, the best of which will be optimized with the MCBH technique. This process is repeated until the final temperature is 1500 K. Below this temperature we observed very little structural change. With a slow fall in the temperature parameter, the MD runs help the MCBH to find the “real” GM. The MCBH runs are performed with the freely available GMIN program,48 which has been used previously to model other oxides.49 The temperature in the MCBH calculations is fixed at 3000 K, and during the search, the step size is varied so that the average acceptance of a Monte Carlo move is 50%. The hybrid SA-MCBH method is very robust. Although requiring a large amount of CPU time, this technique does allow us to perform a thorough search of the possible configurations

Figure 2. Various steps in our GA and hybrid GA shown as a flowchart.

for the (TiO2)n clusters. The configurations obtained with this method are regarded as the target structures for the subsequent GA and hybrid GA methods to find. 2.2.2. Approach 2: Genetic Algorithm (GA). A GA is based on the ideas of Darwin’s theory of evolution. First we assume a fixed environment: in our applications this is the number of ions and the expression for the energy of formation (our model of the interactions). We define any complete set of ionic coordinates for the chosen cluster size to be a “candidate structure”. Within a current population of candidate structures, competition to procreate is simulated. If a structure of one cluster is deemed better than the next, that is, its energy is lower, then statistically, during the simulated competition phase, that candidate is more likely to be chosen as input material when generating the structure of a new candidate cluster. Note that the candidate clusters within the initial population are created using randomized coordinates. Once a new population of candidates is generated, the old population is discarded, the new population becomes the current population, and the next GA cycle begins; see Figure 2. Over the course of several GA cycles a candidate evolves that best fits the environment. In our simulations, as energy is used as the measure of likelihood to procreate, the structure of this candidate cluster will be (if the GA is successful) the GM structure. Our GA is implemented within the GULP code.12,50 Although further details are now presented, a more in-depth explanation and discussion of this procedure can be found in refs 50-53. The GA is used to find the basin on the energy landscape that contains the GM rather than the value of the GM itself; although a well-tuned GA is very efficient at searching different regions of the energy landscape, it is never as efficient as standard optimization techniques when used to relax a structure. Thus, after the GA has finished, the best candidate structure is relaxed in order to obtain the “predicted GM structure”.

15744 J. Phys. Chem. B, Vol. 109, No. 33, 2005 Because any candidate structure within the GM basin (basin containing the GM) is sought, then the GA is only required to find approximate coordinates. Therefore, during the GA phase, we restrict the possible location of ions to a predefined Cartesian grid, where grid points are approximately 0.1 Å apart, which dramatically reduces the energy landscape for the search. The definition of the grid enables a binary representation (the length of which is dependent upon the number of ions and grid points) of the coordinates such that there is a unique one-toone relationship between a configuration of ions and a sequence of binary numbers. During the procreation phase, the GA operators, crossover and mutation, mix and slightly alter these binary sequences, in an way analogous to how DNA is passed on in biological systems, to create new binary sequences and thus new candidate structures that resemble the “parents”. Sometimes the new candidate structures contain the best or worst features of both and therefore are more or less likely to be successful in being selected for procreation in the next GA cycle. For each attempt to find the target or a more stable structure, the GA was initialized 10 times (each starting from a different population of random candidates). Typically each run consisted of 500-1000 GA cycles. For each n, we initially tried one of three predefined grids, each 64 × 64 × 64, but with a different gap between grid points of 0.078, 0.112, and 0.137 Å. If after relaxing the best candidate structure neither the target nor a more stable structure was found, then either a different grid or slight variations of the other GA parameters, as defined in ref 45, were tried. 2.2.3. Approach 3: Genetic Algorithm Hybrid. In a recent review article42 it has generally been found that a GA hybrid is better at finding the GM structure for clusters. Therefore, one option of the GA, implemented with the GULP code, was modified to replicate a similar GA hybrid approach implemented by Roberts and Johnston.41 This included their implementation of the Deaven and Ho54 cut-and-splice versions of the crossover and mutation operators. A hybrid GA is one where each new candidate structure is immediately relaxed (the process of siblings maturing), and the energies of the relaxed candidate structures determine the likelihood of the selection to procreate. Thus the hybrid GA and the MCBH described above are used to search the same landscape, E ˜ C(X). Although the evaluation of each new candidate structure requires more CPU time, the maturing phase does help to remove bad features (for example, by separating overlapping ions) so that there are much fewer poor structures battling to procreate and so a smaller population can be used; in our implementation during the maturing phase we restrict the number of line searches to a maximum of 10, and we use a population between 10 and 20 as well as a reduced number of GA cycles (between 100 and 200). One key procedure when creating a new population is the replacement of identical candidates (in fact, within the current population, any candidate whose energy of formation is within 10-6 eV of another candidate) for new random candidates. These new candidates are then relaxed and checked to make sure a different candidate is indeed generated. As pointed out by Roberts and Johnston,41 any duplicates need to be removed; otherwise the GA struggles as the diversity of the population (a measure of the number of different features that can be mixed during procreation) collapses. In our implementation, as well as the restricted number of line searches, we include the option of moving the ions back onto the defined grid so that binary representation and the corresponding versions of crossover and mutation procedures for generating new candidates can continue

Hamad et al.

Figure 3. Lowest energy and target configurations for TinO2n clusters generated by SA-MCBH method. The energy of formation, E, is in eV/TiO2 unit.

to be employed. We note that, if during the maturing phase, ions within a cluster have moved outside the region covered by the grid, the cluster is displaced so that the center of mass is located near the central grid point. 2.3. Density Functional Theory Calculations. We have also optimized the geometries of the clusters by using different DFT approaches, starting from the stable structures found with the methods described above. There are two main purposes for these calculations: first, we can refine the geometries and also have a more accurate energy of the cluster, and second, the comparison between IP and DFT is a test of the IP used here. In all the DFT calculations, we make use of the hybrid B3LYP exchange-correlation potential55 and different basis sets. First we use an effective core potential (ECP) by Stevens et al. that removes 10 electrons for Ti and 2 electrons for O, combined with a double-ζ basis set.56,57 Then we reoptimize the structure with a double-ζ plus polarization basis set (DZVP) that is specific for DFT.58 For the ECP and the DZVP calculations we also use a Coulomb fitting basis set.58 Finally, we use a more accurate but more costly 6-311G** basis set for a final energy calculation at the DZVP optimized geometry. The comparison among these different basis sets, to estimate their accuracy and to discern whether the ECP choice is accurate enough, will prove useful for future calculations when larger clusters are modeled. All the geometry optimizations have used the NWCHEM59,60 program running in parallel on a cluster. The single point B3LYP/6-311G** calculations employed Gaussian03.61 3. Results and Discussion The target structures, generated by the SA-MCBH method, for the clusters of sizes 1-15 are shown in Figure 3. Many different low energy structures (where currently we only consider the formation energy as defined by eq 1) were also

Structure and Stability of TiO2 Nanoparticles

J. Phys. Chem. B, Vol. 109, No. 33, 2005 15745

Figure 4. Configurations for TinO2n clusters found by GA and/or hybrid GA methods that are more stable, when judged using IP, than that of the respective target structures shown in Figure 3. The energy of formation, E, is in eV/TiO2.

Figure 5. Example higher energy configurations found by GA and/or hybrid GA methods when searching for the GM configuration for (TiO2)n clusters. The energy of formation, E, is in eV/TiO2.

generated by the GA and hybrid GA methods during our search for these target configurations. Three were more stable than their respective target structures and are shown in Figure 4. Less stable configurations, but where the relative energy of formation compared with that of the target configuration is within 0.05 eV/TiO2 unit, are displayed in Figure 5. The energies of formation per TiO2 unit computed using IP for the configurations presented are also shown in Figures 3-5. With the use of the same interatomic potential parameters, the energies of formation for relaxed anatase, brookite, and rutile bulk phases are -39.50, -39.62, and -39.80 eV/TiO2 unit, respectively. Clearly the stability of the clusters increases with size, but is always less stable than that of the bulk phases. In a previous study of (ZnS)n clusters,45 the shape of the GM configurations for the similarly sized clusters were ring-like (regular polygons for n < 6, where rigid ion model used) and then bubble-like. Here, we find that the target clusters are more compact. For example, the ring structure found for n ) 3, shown in Figure 5b, is 0.34 eV higher in energy than the more compact GM configuration for n ) 3, shown in Figure 3c. Upon implementing the hybrid GA approach, three configurationstwo configurations for Ti8O16 (Figure 4b,c) and one for Ti7O14 (Figure 4a)-were found to be more stable than that of the target for the respective cluster sizes. Thus, 13 out of the 15 target structures found using the SA-MCBH approach remained as our predicted GM configurations. However, for n where the SA-MCBH approach was not as successful, the difference in energy between our target and our predicted GM structures is very small (a maximum of 0.022 eV/TiO2). We note that our standard GA approach also readily found the GM configuration for n ) 8 as well as the target structures. The calculated structures are also more stable than those obtained by optimizing small particles with anatase-like struc-

Figure 6. Geometries of TinO2n optimized with B3LYP/DZVP. The optimizations are started from the target structures initially found with SA-MCBH (Figure 5).

Figure 7. Geometries of TinO2n (n ) 7,8) optimized with B3LYP/ DZVP. The optimizations are started from the global minima found with GA or hybrid GA (Figure 4a,b).

tures, such as those studied in ref 30. For example, by optimizing the anatase structures for n ) 16, the smallest particle studied in ref 30, we obtain an energy of -37.018 eV/TiO2 for Ti16O32, which is therefore less stable than our GM clusters with n > 9. The geometries obtained by relaxing the anatase-like particles are less compact than the geometries found using global optimization techniques. This results in less favorable electrostatic interactions as well as bond lengths and angles. For the DFT calculations, we begin by relaxing the GM structures found either with SA-MCBH or GA. Additionally, the Ti7O14 and Ti8O16 structures from SA-MCBH, which according to the IP are not global minima, but differ little from the GM, have been optimized as well. All the DZVP geometries are displayed in Figures 6 and 7, and the total energies of the clusters are summarized in Table 2. The geometries optimized with both ECP plus DZV basis set and with the all-electron plus DZVP basis set are similar except for some small differences in bond lengths and angles. For n ) 1 the DZVP optimization predicts an O-Ti-O angle of 110.6°, whereas ECP and IP predict a linear configuration. We have also optimized the TiO2 molecule with B3LYP/6311G** and MP2/6-311G** methods, which predict angles of 111.41° and 107.90°, respectively. Thus, the correct geometry of the TiO2 molecule seems to be the angular one. The difference between IP, ECP, and DZVP is attributed to the inability of

a The B3LYP/6-311G** energies (last row) are calculated at the B3LYP/DZVP geometries. Most optimizations start from the SA-MCBH coordinates except for n ) 7 and n ) 8, for which the global minima predicted by GA have also been optimized, as specified in the second row.

-89.948 -180.177 -270.404 -360.627 -450.838 -541.0789 -631.241 -631.163 -721.471 -721.421 -811.642 -901.936 -992.084 -1082.292 -1172.564 -1262.783 -1352.958 -999.873 -1999.885 -2999.964 -4000.035 -5000.089 -6000.129 -7000.191 -7000.176 -8000.263 -8000.199 -9000.309 -10000.403 -11000.417 -12000.477 -13000.548 -14000.611 -15000.625 -999.968 -2000.082 -3000.291 -4000.494 -5000.669 -6000.829 -7001.016 -7000.937 -8001.225 -8001.159 -9001.385 -10001.624 -11001.742 -12001.918 -13002.113 -14002.320 -15002.463 ECP DZVP 6-311G**

n ) 12 SA n ) 11 SA n ) 10 SA n)9 SA n)8 GA n)8 SA n)7 GA n)7 SA n)6 SA n)5 SA n)4 SA n)3 SA n)2 SA method n)1 starting geometry SA

TABLE 2: Total Energies of TinO2n Clusters [in hartrees (1 hartree ) 27.212 eV)] for B3LYP/ECP and B3LYP/DZVP Optimized Structuresa

n ) 13 SA

n ) 14 SA

n ) 15 SA

15746 J. Phys. Chem. B, Vol. 109, No. 33, 2005

Hamad et al.

Figure 8. Energy of stabilization, TiO2 f (1/n)TinO2n, in eV, calculated with the different methods used here. IP, interatomic potential; DFT VDZ+ECP, effective core potential and valence double-ζ basis set; DFT DZVP, all-electron double-ζ plus polarization basis set. DFT Geom ) DZVP, E ) 6-311G**, denotes an energy calculation with the 6-311G** basis set for the geometry optimized cluster with the DZVP basis set. All the DFT calculations make use of the B3LYP exchange-correlation potential. The inset is the same graph with a different scale for both axes so that the additional stabilization of Ti10O20 and Ti14O28 can be seen more clearly.

either IP or ECP to account for the double bond character that may occur for an oxygen atom that is bonded to one Ti atom only. This effect becomes less important for larger particles. To gain a more quantitative view of the accuracy of the IP and the ECP approach as compared to the DZVP, we may compare Ti-O distances as in Table 3. There we have included the average distances for each structure as a function of coordination number of Ti and method. The table shows an excellent agreement between distances calculated with IP and the two different DFT methods. IP and ECP tend to overestimate slightly the bond lengths for the low coordination numbers and the smaller particles, but the agreement improves for higher coordination numbers and especially for the larger particles. There are also three structures, Ti5O10, Ti7O14-SA, and Ti8O16SA, for which the 5-fold coordination predicted by IP is not found with either ECP or DZVP. For example, by comparing the structures of Ti5O10 from IP, Figure 3e, and from DZVP, Figure 6e, it can be clearly seen that a 5-fold coordination in IP turns into a 4-fold coordination plus a dangling TidO bond, which is probably due to the preference of the IP for higher coordination numbers, even though that implies straining the structure somewhat, rather than forming partial double bonding. This result is not unexpected given that the potential is derived from TiO2 polymorphs having 4-, 5-, and 6-fold coordination numbers. Nevertheless, for particles with n ) 9-15, the agreement between IP and DFT is good. It is also noteworthy that the DFT GM for Ti7O14 and Ti8O16 are those optimized by starting from the hybrid SA geometries, unlike the IP GM which are those obtained with the GA method. Nevertheless, the energy differences between both types of structures of about 0.05 eV are small. The stability of the particles per TiO2 unit (the energy of the process TiO2 f (1/n)TinO2n) vs n is plotted in Figure 8. The IP agrees reasonably well with the DFT calculations, especially with the energy calculated with the 6-311G** basis set at the DZVP geometry. The DFT energies show that the stability does not decrease monotonically with increasing size, which can be seen more clearly in the inset of Figure 8. There are two sizes,

Structure and Stability of TiO2 Nanoparticles

J. Phys. Chem. B, Vol. 109, No. 33, 2005 15747

TABLE 3: Average Ti-O Distances (in Å) for Different Coordination Numbers (CN) of Ti in TinO2n Particles for Different Methodsa,b CN method 2

1

2

3

4

5

6

7SA

7GA

8SA

8GA

IP 1.728 ECP 1.722 DZVP 1.644 IP 1.823 1.814 ECP 1.813 1.784 DZVP 1.801 1.774 IP 1.878 1.900 1.880 1.890 1.877 1.879 1.876 1.874 ECP 1.894 1.886 1.876 1.877 1.880 1.857 1.867 1.869 DZVP 1.881 1.873 1.865 1.865 1.864 1.857 1.858 IP 1.955 1.939 1.945 1.951 1.953 1.947 ECP 1.932 1.954 DZVP 1.938 1.932 IP 1.988 2.001 ECP 1.998 DZVP 1.997 1.987 1.996

3 4 5 6

9

1.891 1.876 1.874 1.942 1.944

2.021

10

1.863 1.852 1.843 1.947 1.953 1.940

11

1.872 1.877 1.866 1.952 1.941 1.951 1.945 2.002 1.972

12

1.880 1.847 1.866 1.950 1.968 1.942 1.891 1.927 2.007

13

1.867 1.855 1.848 1.933 1.942 1.937 1.921 1.922 1.910

14

1.870 1.860 1.846 1.937 1.956 1.955 1.934 1.937 1.958

15

average

1.874 1.864 1.862 1.946 1.951 1.945 1.944 1.945 1.934

1.728 1.722 1.644 1.818 1.798 1.787 1.878 1.869 1.862 1.946 1.949 1.942 1.946 1.955 1.976

IP ) interatomic potential; ECP ) B3LYP with effective core potential and VDZ basis set; DZVP ) B3LYP with double-ζ plus polarization basis set. b Only global minima are considered, except for n ) 7 and n ) 8, for which the most stable structures for both SA-MCBH and GA are included. For the calculation of the coordination number of Ti, a Ti-O distance cutoff of 2.2 Å has been used. a

Ti10O20 and Ti14O28, that are more stable than their neighbors, which is due to the fact that for these sizes the GM consist of a set of well-connected polyhedra: Ti10O20 can be described as a combination of six tetrahedra and four trigonal bipyramids with only edges and shared apexes, oxygen species, at the surface. However, its neighbors, Ti9O18 and Ti11O22, have less compact structures, each having at least one dangling oxygen atom at the surface. For Ti14O28 there are six tetrahedra, three square base pyramids, three trigonal bipyramids, and two octahedral combined compactly without dangling oxygen species on the outside. Neither Ti13O26 nor Ti15O30 has dangling oxygen species, but to achieve such a compact structure only one central octahedron is allowed, and as many as four and five square base pyramids are formed for Ti13O26 and Ti15O30, respectively. The result is that both clusters are slightly less stable, per TiO2 unit, than Ti14O28. Nonetheless, the stability differences are small, typically a fraction of an electronvolt. Only the more accurate DFT calculations are able to predict these subtle effects. 4. Conclusions We have used a combination of simulated annealing, Monte Carlo basin hopping simulation, and genetic algorithms that provide the global minima of small TinO2n particles (n ) 1-15) with a high level of confidence. These calculations have initially used an interatomic potential that allows a large number of configurations to be explored. The structures have then been refined with more accurate DFT calculations, for which we have also compared several alternatives: effective core potentials, all-electron double- and triple-ζ plus polarization basis sets. The agreement between the IP and DFT is good, particularly for the larger particles (n ) 9-15), for both structure and energy. The calculated global minima consist of compact structures, with titanium atoms reaching high coordination rapidly as n increases. For n ) 11 onward, the particles have at least a central octahedron surrounded by a shell of surface tetrahedra, trigonal bipyramids, and square base pyramids. We have also found that structures with no dangling TidO bonds, a small number of square base pyramids around Ti atoms, and at least one octahedron have a small additional stability, although one needs DFT calculations to capture these subtle effects. The use of these methods opens new possibilities for the study of larger TiO2 particles and of other oxides. Additionally, we report structures that can be used for simulations and modeling

of many interesting properties of this important material, such as quantum size effects on photophysical properties, adsorption on rather small nanostructured TiO2, sintering, and nucleation. Acknowledgment. We would like to thank the EPSRC Network for Crystal Growth and Nucleation. J.A.M. and S.L. thank the Spanish Ministerio de Educacio´n y Ciencia (Project Nos. VEM2003-20574-C03-01 and CTQ2004-00582/BQU). References and Notes (1) Gerisher, H. Photochem. Photobiol. 1972, 16, 243. (2) Fujishima, A.; Honda, K. Nature 1972, 238, 37. (3) Fujishima, A.; Honda, K. Bull. Chem. Soc. Jpn. 1971, 44, 1148. (4) Hoffmann, M. R.; Martin, S. T.; Choi, W.; Bahnemann, D. Chem. ReV. 1995, 95, 69. (5) Volodin, A. M. Catal. Today 2000, 58, 103. (6) Tanaka, K.; Hisanaga, T.; Rivera, A. P. In Photocatalytic Purification and Treatment of Water and Air; Al-Ebaki, H., Ollis, D. F., Eds.; Elsevier: Amsterdam, 1993; pp 169-178. (7) Mejias, J. A.; Jimenez, V. M.; Lassaletta, G.; Fernandez, A.; Espinos, J. P.; Gonzalez-Elipe, A. R. J. Phys. Chem. 1995, 100, 1484. (8) Zhang, H.; Penn, R. L.; Hamers, R. J.; Banfield, J. F. J. Phys. Chem. B 1999, 103, 4656. (9) Zhang, Z.; Wang, C.; Zakaria, R.; Ying, J. Y. J. Phys. Chem. B 1998, 102, 10871. (10) Chun, H.; Yizhong, W.; Hongxiao, T. Appl. Catal. B: EnViron. 2001, 35, 95. (11) Kang, M. J. Mol. Catal. A: Chem. 2003, 197, 173. (12) Du, Y.; Rabani, J. J. Phys. Chem. B 2003, 107, 11970. (13) Aran˜a, J.; et al. J. Mol. Catal. A: Chem. 2003, 197, 157. (14) Zhang, Z.; et al. Langmuir 2003, 19, 8230. (15) Irie, H.; Watanabe, Y.; Hashimoto, K. J. Phys. Chem. B 2003, 107, 5483. (16) Kasuga, T. Langmuir 1998, 14, 3160. (17) Chen, Q.; Zhou, W.; Dou, G.; Peng, L. AdV. Mater. 2002, 14, 1208. (18) Ma, R.; Bando, Y.; Sasaki, T. Chem. Phys. Lett. 2003, 380, 577. (19) Tian, Z. R. J. Am. Chem. Soc. 2003, 125, 12384. (20) Chen, L. X.; Rajh, T.; Wang, Z.; Thurnauer, M. C. J. Phys. Chem. B 1997, 101, 10688. (21) Yeung, K. L.; et al. J. Phys. Chem. B 2002, 106, 4608. (22) Rajh, T.; et al. J. Phys. Chem. B 2002, 106, 10543. (23) Grubert, G.; Stockenhuber, M.; Tkachenko, O. P.; Wark, M. Chem. Mater. 2002, 14, 2458. (24) Choi, H. C.; Jung, Y. M.; Kim, S. B. Vib. Spectrosc. 2005, 37, 33. (25) Persson, P.; Christof, J.; Lunell, S. J. Phys. Chem. B 2003, 107, 3336. (26) Barnard, A. S.; Zapol, P. J. Chem. Phys. 2004, 121, 4276. (27) Barnard, A. S.; Zapol, P. Phys. ReV. 2004, 70, 235403. (28) Barnard, A. S.; Zapol, P. J. Phys. Chem B 2004, 49, 18435. (29) Barnard, A. S.; Zapol, P.; Curtiss, L. A. J. Chem. Theory Comput. 2005, 1, 107. (30) Redfern, P. C.; Zapol, P.; Curtiss, L. A.; Rajh, T.; Thurnauer, M. C. J. Phys. Chem. B 2003, 107, 11419.

15748 J. Phys. Chem. B, Vol. 109, No. 33, 2005 (31) Tilocca, A.; Selloni, A. J. Chem. Phys. 2003, 119, 7445. (32) Gillan, M. J.; Lindan, P. J. D.; Kantorovich, L. N.; et al. Mineral. Magn. 1998, 62, 669. (33) Lindan, P. J. D.; Harrison, N. M.; Gillan, M. J.; White, J. A. Phys. ReV. B 1997, 55, 15919. (34) Calzado, C. J.; San Miguel, M. A.; Sanz, J. F. J. Phys. Chem. B 1999, 103, 480. (35) Miguel, M. A. S.; Calzado, C. J.; Sanz, J. F. Int. J. Quantum Chem. 1998, 70, 351. (36) Giordano, L.; Pacchioni, G.; Bredow, T.; et al. Surf. Sci. 2001, 471, 21. (37) Matsui, M.; Akaogi, M. Mol. Simul. 1991, 6, 239. (38) Oliver, P. M.; Watson, G. W.; Kelsey, E. T.; Parker, S. C. J. Mater. Chem. 1997, 7 (3), 563-568. (39) Bandura, A. V.; Kubicki, J. D. J. Phys. Chem. B 2003, 107, 1107211081. (40) Hartke, B. Angew. Chem., Int. Ed. 2002, 41, 14683, and references within. (41) Roberts, C.; Johnston, R. L. Phys. Chem. Chem. Phys. 2001, 3, 5024. (42) Hartke, B. In Structure and Bonding; Johnston, R. L., Ed.; SpringerVerlag: Heidelberg, 2004; Vol. 110, p 33. (43) Coley, D. A. An introduction to genetic algorithms for scientists and engineers; World Scientific: Singapore, 1999. (44) Kirkpatrick, S.; Gellat, J. C. D.; Vecchi, M. P. Science 1983, 220, 671. (45) Woodley, S. M.; Sokol, A. A.; Catlow, C. R. A. Z. Anorg. Allg. Chem. 2004, 630, 2343. (46) Wales, D. J.; Scheraga, H. A. Science 1999, 285, 1368. (47) Smith, W.; Forester, T. R. J. Mol. Graphics 1996, 14, 136. (48) Wales, D. J. GMIN: A program for finding global minima; http:// www-wales.ch.cam.ac.uk/software.html. (49) Flikkema, E.; Bromley, S. T. J. Phys. Chem. B 2004, 108, 9638. (50) Woodley, S. M.; Battle, P. D.; Gale, J. D.; Catlow, C. R. A. Phys. Chem. Chem. Phys. 1999 1, 2535. (51) Woodley, S. M.; Battle, P. D.; Gale, J. D.; Catlow, C. R. A. Phys. Chem. Chem. Phys. 2004, 6, 1815. (52) Woodley, S. M. Phys. Chem. Chem. Phys. 2004, 6, 1823. (53) Woodley, S. M. In Structure and Bonding; Johnston, R. L., Ed.; Springer-Verlag: Heidelberg, 2004; Vol. 110, p 95. (54) Deaven, D. M.; Ho, K. M. Phys. ReV. Lett. 1995, 75, 288.

Hamad et al. (55) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (56) Stevens, W. J.; Basch, H.; Krauss, M. J. Chem. Phys. 1984, 81, 6026. (57) Stevens, W. J.; Krauss, M.; Basch, H.; Jasien, P. G. Can. J. Chem. 1992, 70, 612. (58) Godbout, Salahub, D. R.; Andzelm, J.; Wimmer, E. Can. J. Chem. 1992, 70, 560. (59) Straatsma, T. P.; Apra`, E.; Windus, T. L.; Bylaska, E. J.; de Jong, W.; Hirata, S.; Valiev, M.; Hackler, M. T.; Pollack, L.; Harrison, R. J.; Dupuis, M.; Smith, D. M. A; Nieplocha, J.; Tipparaju, V.; Krishnan, M.; Auer, A. A.; Brown, E.; Cisneros, G.; Fann, G. I.; Fruchtl, H.; Garza, J.; Hirao, K.; Kendall, R.; Nichols, J.; Tsemekhman, K.; Wolinski, K.; Anchell, J.; Bernholdt, D.; Borowski, P.; Clark, T.; Clerc, D.; Dachsel, H.; Deegan, M.; Dyall, K.; Elwood, D.; Glendening, E.; Gutowski, M.; Hess, A.; Jaffe, J.; Johnson, B.; Ju, J.; Kobayashi, R.; Kutteh, R.; Lin, Z.; Littlefield, R.; Long, X.; Meng, B.; Nakajima, T.; Niu, S.; Rosing, M.; Sandrone, G.; Stave, M.; Taylor, H.; Thomas, G.; van Lenthe, J.; Wong, A.; Zhang, Z.; NWChem, A Computational Chemistry Package for Parallel Computers, version 4.6; Pacific Northwest National Laboratory: Richland, WA, 2004. (60) Kendall, R. A.; Apra`, E.; Bernholdt, D. E.; Bylaska, E. J.; Dupuis, M.; Fann, G. I.; Harrison, R. J.; Ju, J.; Nichols, J. A.; Nieplocha, J.; Straatsma, T. P.; Windus, T. L.; Wong, A. T. High Performance Computational Chemistry: an Overview of NWChem a Distributed Parallel Application. Comput. Phys. Commun. 2000, 128, 260-283. (61) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision B.04; Gaussian, Inc.: Pittsburgh, PA, 2004.