Computational characterization of the structural and mechanical properties of nanoporous titania

Nanoporous titania is one of the most commonly used biomaterials with good biocompatibility and mechanical strength. Understanding to the influence of pore structures on their performances is crucial for the design and preparation of titania-based materials. Two kinds of structural models for nanoporous titania were constructed and used to investigate the effect of pore size and/or porosity on their mechanical properties by using molecular dynamic simulations with the Matsui–Akaogi potentials. The porous structures were relaxed and their elastic constants were computed and used to evaluated their bulk, shear and Young's moduli. Overlap effect in small pores, pore size and porosity have considerable influence on computed elastic moduli. Compared to bulk rutile TiO2, reduced mechanical moduli were predicted. Simulations on uniaxial tensile tests revealed an anisotropic stress–strain relationship and a brittle-to-ductile transition for structures with large porosities. Fracture failure was predicted for all the studied porous structures. The maximum stress decreases with pore size and porosity, while the corresponding strain decreases with pore size, but increases with porosity.


Introduction
Porous materials can be classied into three types by their pore sizes: microporous (<2 nm), mesoporous (2-50 nm) and macroporous materials (>50 nm). Nanoporous materials are dened as those porous materials with pore diameters of 1-100 nm. 1 Nanoporous materials have shown great potential in catalysis, hydrogen storage, drug delivery and biological medicine, etc. [2][3][4][5] because of their high surface-to-volume ratios, moderate and tunable pore sizes, pore shapes, porosity and framework compositions. Being used as biomedical biomaterials for tissue engineering and drug delivery, some porous materials with good biocompatibility and bioreductivity have been characterized, such as porous titanium and porous calcium phosphates. In fact, implant three-dimensional substitute materials with porous structures are widely used in the treatment of bone trauma in clinics because of their good osteoreductivity that is essential for bone regeneration. [6][7][8][9] The pore structures provide ideal anchors for biomolecules/cell adhesion, proliferation, differentiation and formation of new organs. On the other hand, the porous materials by careful design may possess mechanic strengths that are sufficient to support damaged organs or tissues. 10,11 Typically, most biological implant materials for bone repair have a Young's modulus higher than those of hard tissues. Stress occlusion can impede bone remodeling and healing, resulting in increased bone porosity. One of the major problems associated with implants in orthopedic surgery is the mismatch in Young's modulus between natural bone and implants. One way to alleviate this problem is to reduce the Young's modulus of the materials by introducing holes, [12][13][14] thereby minimizing damage to tissues near the implant and ultimately extending device lifetime. In addition, bone/implant xation is achieved by the mutual intersection between the bone and the porous implant matrix. 15 Among various bone gra substitutes from natural and synthetic materials, titanium dioxide (TiO 2 ) is considered to be a highly suitable biocompatible material for bone anchoring implants due to its good biocompatibility, osteoconductivity and osteoinductivity. [16][17][18][19] In fact, TiO 2 exhibits little or no toxicity in both in vitro and in vivo experiments. 20 Extensive studies from their preparation, characterization, through applications have been conducted for various TiO 2 -based materials. Fattakhova-Rohlng et al. 21 surveyed the synthesis of various porous TiO 2 nanostructures including membranes, spheres, bers, etc. Sabetrasekh et al. 22 reported that porous TiO 2 particles exhibit better activity in cell proliferation, but less toxicity than commercial bone gras. Natix®, Straumann® BoneCeramic and Bio-Oss®. Haugen et al. 23 developed a highinterconnectivity and high-porosity TiO 2 that favors osseous integration, union and regeneration. Santos et al. 24 explored the mechanical properties of biomedical titanium dioxide lms prepared by anodizing on commercial pure Ti (cp-Ti), including nanohardness, elastic modulus and brittleness. The modied lms have relatively high nano-hardness and small elastic modulus in comparison with the cp-Ti. The reduction in modulus of elasticity was attributed to the porosity and inherent roughness of the oxide surface.
A number of computational studies have been conducted to explore the structures and properties of rutile titania and to help with the design of rutile-based materials. The computations focused mainly on the electronic, thermodynamic and mechanical properties of bulk or nano-sized TiO 2 . For example, Mashreghi 25 calculated the thermal expansion coefficients of rutile titania nanoparticles under 300-1000 K by means of molecular dynamics (MD) simulations with the Buckingham potentials. Kim et al. 26 reproduced the lattice structures, thermal expansion coefficients and bulk modulus of rutile, brookite and anatase by using molecular dynamics simulations with the Morse interactions. Chung and Buessem 27 used Voigt-Reuss-Hill (VRH) approximation to calculate the elastic modulus of isotropic polycrystalline TiO 2 based on the anisotropy single crystal elastic constant. The results are in good agreement with the measured polycrystalline elastic modulus. Mahmood et al. 28 performed a rst-principles calculation of rutile TiO 2 using the general gradient approximation (GGA) proposed by Perdew-Burke-Ernzerhof (PBE), and found that the single crystal elastic constant (C ij ), bulk modulus (B), Young's modulus (Y), shear modulus (G), linear bulk modulus (B a ) and the corresponding single crystal elastic constant (S ij ) alternately increase or decrease with pressure. However, most of those studies focused on non-porous crystalloid rutile TiO 2 . In the research of mesoporous titania, Solveyra et al. 29 performed molecular dynamics simulations of water conned in mesoporous TiO 2 -rutile pores at different water contents, obtaining water density and diffusion coefficient as a function of distance from the surface. Similarly, Velasco et al. 30 performed a detailed analysis of the hydrodynamics in the pore space of titanium dioxide by a combination of experimental nuclear magnetic resonance NMR data and MD simulation. Gautam et al. 31 used MD simulation to study the structural and kinetic properties of propane restricted in mesoporous TiO 2 .
Compared to tremendous studies on crystallized titanium dioxides, fewer researches have been conducted to their porous samples. Precise design and control over the pore structures is crucial for modulating the properties of porous TiO 2 . Computational studies provide an effective way to establish the correlation between pore structures and performances. In this work we designed various computational models for porous titania with different pore sizes and porosities, and studied their structural and mechanical properties. Our purpose is to gure out the inuence of pore structures on the properties of porous titania, which would be helpful for the design and preparation of titania-based materials.

Methodology
Titania has three phases in nature, rutile, anatase, and brookite. 32 Rutile TiO 2 is thermodynamically stable in its bulk form with a large crystallite size at normal pressure and temperature up to its melting point 1830 K. 33 Rutile belongs to P4 2 /mnm space group in which a titanium atom sits at the center of the unit cell and is surrounded by an octahedron formed by six oxygen atoms. The centered Ti atom has six coordinates, while each O atom is coordinated by three Ti atoms. Starting with the rutile TiO 2 structure, two series of computational models for porous titania were designed and constructed by excavating one cylinder along the z-axis direction. In Model I, all of the rutile TiO 2 porous structures were established at a xed porosity of 8.1% and the cylinders are 1.3, 2.8, 3.4 and 5.1 nm in diameter, labeled as IA, IB, IC and ID, respectively. This model was used to evaluate the effect of pore size on the structures and properties of porous titania. In Model II, cylinders with a xed diameter of 2.8 nm were excavated in the supercell. The porosity (from 8.1% to 22.6%) was controlled by changing the number of unit cells in the supercell. This model was used to evaluate the effect of porosity on the structures and properties of porous titania, labeled as IIA, IIB, IIC, IID and IIE, respectively. In both models, the Ti/O ratio is kept to 1 : 2 (ref. 29) and the wall thickness between adjacent pores are larger than 1 nm, as shown in Fig. 1. The three-dimensional structures of Model I and Model II were also given in Fig. S1. † The structural parameters of the constructed models are presented in Table  S1. † Arising from the extensive research interest, several kinds of force elds have been proposed for titanium dioxide. [34][35][36][37][38][39][40][41][42][43][44][45][46] The Matsui-Akaogi (MA) force elds 46 are commonly used to describe the TiO 2 systems. Other potentials used for the TiO 2 systems include modied-MA 42,43 and MS-Q. 44,45 Using molecular dynamics simulations with the MA force elds, Matsui and Akaogi 46 calculated the structure and physical properties of four kinds of TiO 2 . Naicker et al. 36 and Fuertes et al. 38 simulated the surface energy and optical properties of TiO 2 nanoparticles using MD simulations with the MA potential. The MA potential was also used to simulate the surface structure of TiO 2 (ref. 35) and the deformation of anatase-type nanocrystals under uniaxial stretching and compression at room temperature. 37 In addition, the MS-Q potential 44,45 also produced accurate descriptions for the structures and properties of TiO 2 crystal. [39][40][41] The MA potential is dened as where U ij (r), C ij and r ij are the interaction energy, van der Waals (vdW) coefficient and distance between atoms i and j, respectively. r and q are ionic-pair dependent length parameter and atomic charge, respectively. 3 0 is the dielectric constant in vacuum. Alderman and Skinner et al. 43 further proposed a modied-MA potential as where s is an interaction-dependent length parameter. The MS-Q potential model is dened as where A ij , r 0 , and B ij are the parameters of Morse potential. r 0 is the equilibrium interatomic separation. The above mentioned three potentials were employed in this work. The parameters are listed in Table S2-S4 in ESI. † Molecular dynamics simulations were carried out using the LAMMPS code. 47 The time integral of the Newton's equation of motion was performed using the velocity Verlet algorithm, and the time step was set to 1.0 fs. The Nosé-Hoover thermostat 48,49 was used to control the temperature. The periodic boundary conditions were applied in the three directions. All the structures were rstly relaxed with the conjugate gradient (CG) algorithm, 50 and then the constant-pressure constanttemperature (NPT) ensemble was employed for simulations at given temperatures. Particle-particle particle-mesh (PPPM) method was applied to minimize errors in long range terms coulombic. The cut-off of electrostatic interaction and van der Waals interaction used in all the simulation was set to 10Å. Each simulation was carried out for 1.2 ns. The variations in potential energy with simulation time of all the structures are shown in Fig. S2-S4 in ESI. † All the systems reach their equilibrium at about 200-500 ps. The simulations continue and the data of last 400 ps were used for analysis.
The elastic modulus of rutile TiO 2 are computed with the elastic constants, which are obtained using the energy-strain theory: and where V 0 and E(V 0 ,0) represent the volume and energy of the unstrained structure. m, d and C ij denote the stress tensor, strain tensor and elastic constants respectively. Six independent elastic constants, C 11 , C 12 , C 23 , C 33 , C 44 and C 66 , are dened for the tetragonal bulk and mecroporous rutile TiO 2 structures. 51 The maximum strain was set to 1% in the evaluation of elastic constants. Based on the constants, bulk modulus (K), shear modulus (G) and Young's modulus (E) and Poisson ratio (h) were then computed with the Voigt-Reuss-Hill method. [52][53][54] The calculations of elastic constants were carried out for the equilibrated structures at 300 K. The details of K, G, E, h calculations are given in the ESI. † The volume thermal expansion coefficient where V 0 and T 0 represent the initial reference volume and temperature. V T represents the volume at temperature T. The b variations of titanium dioxide were explored at temperatures between 300 K and 1000 K with a temperature separation of 100 K.

Results and discussion
First, we examined the reliability of the potentials by examining the structures and properties of bulk rutile TiO 2 . A 9 Â 8 Â 18 supercell with 7776 atoms was established. The simulated lattice parameters are given in Table 1, together with the measured and calculated results in previous studies. 55,56 The three potentials produce comparable results, and the MA and MS-Q results match well with the measurements. 56 The computed six elastic constants of rutile TiO 2 are given in Table 2. Previous calculations 33,57 produced similar results with experiments. 58 Our simulations with the MA and the MS-Q potentials are in good agreement with other studies. The modied-MA results deviate considerably from the The MS-Q potential was thus selected to calculate the volume thermal expansion coefficients of the porous structures. Fig. S5 and S6, † which have similar shapes with Fig. 1, in ESI, † show the pore structures with different pore sizes and porosities in rutile TiO 2 aer relaxation with the MA force eld. Compared to Fig. 1, the atoms on the surface of the pores are rearranged to some extent. Both the surface oxygen and titanium atoms basically retain their interatomic connections with only small displacements around their original locations. The locations of inner atoms far away from the holes are almost unchanged. Table 3 lists the structural parameters of the pores in Model I and II. Compared to the unrelaxed structures (Table S1 †), in both models the total volumes of the cubic boxes decrease aer relaxation under the isobaric-isothermal ensemble. Accordingly, the pore sizes become smaller, but the changes are less than 4%. Atoms have a tendency to move toward the vacuum holes, resulting to the reduction of pore size and the distortion of pore edges. In both models, moreover, the surface areas of the pores and porosity also have a little decrease. Table 3 also presents the specic surface area (g) of all the porous models. g decreases with pore size in Model I, and increases with porosity in Model II. Fig. 3 shows the volume thermal expansion coefficient (b) of the porous rutile TiO 2 structures. In both models, b increases with temperature, which is in agreement with many other materials. At some temperatures, b increases with pore size and porosity, but, in general, their changes are rather small. Moreover, the b values are comparable with that of bulk rutile TiO 2 .
Six independent elastic constants, C 11 , C 12 , C 23 , C 33 , C 44 and C 66 were computed for the porous structures with the MA force elds, which are given in Table 4. The mechanical stability of these structures was examined using the following conditions: 61,62 C 11 > 0, C 33 > 0, C 44 > 0, C 66 > 0, C 11 -C 12 > 0, C 11 + C 33 À 2C 23 > 0, and 2C 11 + C 33 +2C 12 + 4C 23 > 0. These conditions are satised by the computed results of all the structures. The results of bulk rutile TiO 2 are also given for comparison. For Model I, the elastic constants of R ¼ 1.3 nm (IA) have rather large changes compared to the corresponding bulk ones. Because of the serious overlap effect of surface atoms in the small pore, IA exhibits different features from the other models. The overlap effect reduces in the models from IB, IC to ID. The computed elastic constants are comparable for IB and IC, and change to some extent for ID. Moreover, changes in C 11 , C 12 and C 66 are more remarkable than in C 23 , C 33 and C 44 . For Model II, the elastic constants decrease with increasing porosity. Similarly, the decrease in C 11 , C 12 and C 66 are more remarkable than in C 33 , C 23 and C 44 . Therefore, both pore size and porosity have inuence on the magnitude of elastic constants, especially in C 11 , C 12 and C 66 .
Bulk modulus (K) and shear modulus (G) are used to measure the rigidity of a material resisting to the volume deformation and shape deformation under external forces. Young's modulus (E), a proportion of stress and strain, is customarily used as a measure of rigidity of a crystalline solid. Fig. 4 and 5 display the variation of these moduli with pore diameter and porosity, respectively. The moduli of bulk rutile TiO 2 are also given in the gures. At the same pore diameter or porosity, an order of E > K > G is always noted. Those three moduli are smaller than the corresponding bulk values. Except IA, the three moduli decrease with pore size and porosity, which again illustrates the overlap effect in small nanopores. The reduced moduli in porous rutile TiO 2     are in agreement with the observations for the rigidity reduction of other porous materials. For example, Li et al. 63 veried that the elastic moduli of porous nickel decrease as the pore size increases. The reduction was attributed to the increased load and moment on the pore walls. It is interesting to note that the decrease of E and G with porosity (Fig. 5) is approximately between linear and quadratic functions. Kováčik et al. 64,65 proposed two models of permeability, which are given as where E and G are the effective Young's modulus and shear modulus of porous material, respectively with the porosity of P. E 0 and G 0 are the Young's moduli and shear moduli of bulk material. P c is the porosity at which E or G becomes zero. f E and f G are the characteristic exponents of the porous samples. In general, the characteristic indices f E and f G are 1.10-1.70 and 1.00-1.40, complying with the computed trend shown in Fig. 5. The ductility/brittleness variations have been observed in other porous materials. Meille et al. 66 reported that porous alumina, which has a xed pore size, becomes brittler when its porosity increases. The variation was attributed to the thinner walls in the structures with larger porosities, which are more difficult to resist against shear loads. Poisson's ratio (h) represents the ratio of transverse strain to axial strain when a material is deformed along the direction of load. It is an elastic constant reecting transverse deformation of a material. The material with a h indicates that the lateral deformation amount is larger than the longitudinal deformation before the plastic deformation occurs. On the contrary, the lateral deformation is smaller than the longitudinal deformation amount. The value of Poisson's ratio (h) is oen used to estimate the ductility/ brittleness with a critical point of 0.25 or lower for a brittle material and 0.33 or higher for a ductile material. 67,68 Bulk rutile TiO 2 is brittle, but little is known about the brittleness changes when the pores are introduced. The h values are between 0.26 and 0.31 for the structures in Model I and between 0.28 and 0.22 for the structures in Model II. However, one notes that the h values were computed by assuming the structures are isotropic, but the porous rutile TiO 2 are anisotropic. Therefore, the estimation with isotropic h values is not feasible for predicting the ductility/brittleness of the porous structures.
Uniaxial tensile measures the mechanical responses of a material subjected to a uniaxial load in an external environment. Molecular dynamics simulations have been conducted by many authors 69-72 to investigate the elastic limit, tensile strength, yield point, yield strength and other tensile properties of crystals, polymers and ceramics. The stress-strain behavior of porous rutile TiO 2 was simulated in this work by performing uniaxial tensile tests on the model structures. The simulations were carried out on the x, y and z directions of rutile TiO 2 , which are perpendicular to the (001), (100) and (010) facets, respectively. In each stress-strain relationship, only one direction is stretched, while pressure in the other two vertical directions was set to zero. The strain rate of stretching was set to 1.0 Â 10 10 s À1 and temperature was kept at 298 K. The stress-strain curves were obtained for all the structures in Model I and Model II, as shown in Fig. 6. The structures were stretched from 0 to 20% under an external load, and the corresponding stresses are s xx, s yy and s zz in the three directions, respectively. The uctuations in the two directions perpendicular to the tensile direction are rather small and are not presented.
For the structures in Model I, their stress-strain curves follow the same behaviors: elastic at small strains, and then failure at great strains. However, the curves vary with pore size and are different from those of bulk rutile TiO 2 . The porous structures have smaller elastic moduli than the bulk and their failure points are smaller than those of the bulk in the three directions. In the x and y directions, the structure with a large pore has a slightly smaller elastic modulus than the structure with a small pore. The maximum stresses s xx of Model IA-ID in the x direction are 17.4, 16.5, 15.4 and 12.7 GPa, the strains corresponding to the maximum stresses are 15.2, 13.6, 12.5 and 10.4%, respectively. In the y direction, the maximum stresses are 17.5, 16.2, 15.7 and 12.6 GPa and their corresponding strains are 16.0, 13.5, 12.9 and 10.3%, respectively. The maximum stresses s xx and s yy , as well as their corresponding strains decrease with pore size, indicating that the structures with larger pore sizes at the same porosity have lower tensile strengths. Therefore, the structure with a large pore tends to be fractured at a small load in the x and y directions. In the z direction, the maximum s zz values are larger those in the other two directions and have few changes with pore size. Their corresponding strains also have rather small changes, 12.7%, 12.5%, 12.7% and 12.6%, respectively. In the direction parallel to the cylindrical pore, therefore, the mechanical behaviors have smaller changes than the other two directions.
One notes that the curves for Model IA in the x and y directions uctuate when the loads are close to the maximum stresses, which indicates a plastic deformation at a strong stress. Such plastic deformation was not noted in the other three structures in Mode I. IA has the smallest pore size and thinnest wall among the four structures when their porosity is xed. The thin wall favors a plastic deformation for the rutile TiO 2 . However, the plastic deformation changes into a brittle fracture when the strain increases further. Therefore, the stretching failure in all the structures in Model I is characterized by their brittle fracture, a sudden drop in stress at a given strain.
For the structures in Model II whose pore sizes are xed and porosities are varying and larger than those of the structures in Model I, their stress-strain curves reect elastic responses at small strains. However, the curves deviate from the elastic responses at larger strains. The stresses uctuate with strain, which represents a plastic deformation at this stage. The deformation is more remarkable for the structures with large porosities, resulting from their thinner wall thickness. Brittle fractures occur when the stains increase further. Over the studied stain range, elastic and inelastic, the moduli of the porous structures are smaller than those of the bulk. The stressstrain relationship varies with porosity. In the x direction, the maximum stress reaches 16.73, 14.87, 13.9, 12.09 and 12.02 GPa at a strain of 13.68%, 14.2%, 14.78%, 14.58% and 15.54%, respectively. The maximum stress decreases with porosity. The structure with a larger porosity has a lower maximum stress, which, however, corresponds to a larger strain. This is in contrast to the variations found in Model I. The structures in Model II have larger porosities than those on Model I. When the external load is larger than 10%, the stress varies in an inelastic pattern with the strain. The deviation from elastic response is more remarkable for the structures with larger porosities. The ve structures are in the order of IIA > IIB > IIC > IID > IIE. The brittle-to-ductile transition results from the rearrangement of atoms around the pores, leading to a small stress/stain over a certain range of strains. Therefore, IIE has the smallest maximum stress but the largest failure strain. Similar variations were noted in the stress-strain curves in the y direction, the maximum stress and the failure strain vary with porosity, but their maximum stresses are lower than those in the x direction. The brittle-to-ductile transition was also noted. In the z direction, the variations with porosity become less remarkable. The ve structures have almost the same stress-strain relationship at small strains. Their maximum stresses are higher than those in the other two directions, and are slightly higher for the structures with larger porosities. The maximum stresses correspond to almost the same strain values. The brittle-to-ductile transition is not clear in this direction, but brittle fractures are noted for all the ve structures.

Conclusion
Two structural models were constructed to mimic the structures of porous rutile TiO 2 . In the rst model, the four structures have different pore sizes of 1.3, 2.8, 3.4 and 5.1 nm in diameter, respectively, but have the same porosity. In the second model, the ve structures have different porosities from 8.1%, 12.7%, 15.9%, 20.5% and 22.6%, respectively, but have the same pore size. The inuence of pore size and porosity on the structures and mechanical properties of porous rutile TiO 2 were then investigated using MD simulations. The MA potential was selected to predict the structures and elastic properties, and the MS-Q potential was used for computing the volume expansion coefficients.
In the relaxed porous structures, the surface atoms basically retain their interatomic connections with only small displacements around their original locations. The locations of inner atoms far away from the holes are almost unchanged. Their volume thermal expansion coefficients are comparable with that of bulk rutile TiO 2 , and increase with pore size and porosity at the same temperature, but the increments are rather small. The computed elastic constants vary with pore size and porosity, especially for C 11 , C 12 and C 66 . The structure with the smallest pore size of 1.3 nm exhibit exceptional variations in its elastic constants because of the overlap effect in its small pore. The elastic modulus (E, K and G) were evaluated based on the computed elastic constants. The three moduli are smaller than the corresponding bulk values, decreasing with pore size and porosity.
The uniaxial tensile tests were performed for all the structures in the two porous models with MD simulations. The stress-strain curves were obtained in three directions. The stress varies in an anisotropic way. In the z direction parallel to the cylindrical pore the stress behaves close to that of the bulk. In the other two directions, the maximum stress decreases with pore size and porosity. At small strains, elastic response was noted for all the structures. When the strain increases, inelastic response was noted for the structures with large porosities. The deviation from elastic response increases with porosity, and brittle-to-ductile transition was noted. When the strain increases further, brittle fracture where a sudden drop in stress occurs was noted for all the structures. The brittle failure occurs at smaller strains for the structures with larger pores in Model I, but at larger strains for the structures with large porosities, resulting from the brittle-to-ductile transitions. Our simulations show that the mechanical behaviors of porous structures depend closely with pore structures, including their size, density, etc. Pore shape and their connectivity may also have great inuence on their mechanical properties, which will be investigated in our future work.

Conflicts of interest
There are no conicts to declare.