Development of novel interatomic potentials for simulation of rutile TiO2

Physica B: Condensed Matter 574 (2019) 311657

Contents lists available at ScienceDirect

Physica B: Condensed Matter journal homepage: www.elsevier.com/locate/physb

Development of novel interatomic potentials for simulation of rutile TiO2 ∗

T

Yu Tian, Jun Ding, Xia Huang , Kun Song, Shi-qing Lu, Hao-ran Zheng College of Mechanical Engineering, Chongqing University of Technology, Banan, 400054, China

A R T I C LE I N FO

A B S T R A C T

Keywords: TiO2 crystals Interatomic pairwise potentials Molecular dynamic simulation Mechanical properties

A novel interatomic potential for molecular dynamic (MD) simulation of rutile TiO2 crystals was developed based on Buckingham pairwise potentials. MD simulation of rutile TiO2 was conducted using the developed potentials to obtain various physical parameters (e.g., lattice constant, bulk modulus, shear modulus, and elastic modulus) of rutile TiO2 crystals. Results based on the developed potentials were then compared with the experimental results and first-principles calculations. The results indicated that the developed interatomic potentials can accurately calculate various physical parameters (e.g., lattice constant, bulk modulus, shear modulus, and elastic modulus) of rutile TiO2 crystals. Calculated based on the developed potentials, the melting point of rutile TiO2 was highly consistent with the experimental values which demonstrates the validity of the developed potentials.

1. Introduction Under ambient pressure, TiO2 crystals exhibit three phases, including rutile, anatase, and brookite phases. Owing to its high catalytic activity and good stability, rutile TiO2 has shown potential as a photocatalytic material [1]. However, rutile TiO2 is only active under UV light with a wavelength below 386 nm because of its relatively broad bandgap. UV light constitutes less than 5% of sunlight, hence the extremely low utilization rate of sunlight by rutile TiO2. Wang et al. [2] modified TiO2 by doping with metallic ions such as Fe3+, Co2+, Cu2+, and Ni2+ to introduce defect energy levels to its bandgap. The response of TiO2 to visible light was enhanced, and its photocatalytic performance was improved. Alamgir et al. [3] prepared 5% Co-doped TiO2 nanostructures by using the sol-gel method. Raman spectroscopy showed that the bandgap was 3.33 eV, and the increased bandgap led to enhanced photocatalytic performance. Villabona et al. [4] modified TiO2 by doping with rare earth ions to alter the crystallinity of TiO2. Electron-hole traps were generated to separate electron–hole pairs and extend the life of photogenerated charges, thus improving the photocatalytic performance of the doped systems. Zhu et al. [5] proposed the photoelectrochemical measurements of black TiO2 with the nucleus–amorphous shell structure. The results demonstrated significantly increased the carrier density and improved the visible light absorption of black TiO2 in photocatalysis relative to those of conventional white TiO2. Despite the significantly enhanced photocatalytic performances of TiO2, the aforementioned methods require accurate measurement of the

∗

radius, valence, and doping concentration of the doping ions as these physical parameters significantly affect the catalytic performance of doped systems. The optimized doping concentration in experiments is extremely challenging to achieve; and inappropriate doping concentrations may cause degradation of the photocatalytic performance. Thus, systematic studies on the stability, physical, and chemical properties of TiO2 are crucial to enhance its photocatalytic performance. Experiments can effectively obtain several macroscopic physical properties of TiO2 but cannot fully elucidate the microscopic mechanism behind these properties. Molecular dynamic (MD) simulation is an effective method to investigate the microscopic mechanisms of materials at the atomic scale. One of the key factors affecting the accuracy and validity of MD simulation results is the selection of potentials [6]. Ti–O and O–O interatomic potentials were developed in the study of defect pyrochlore oxide structures [7]. However, all of these interatomic potentials were based on the shell model and not suitable for displacement cascade simulations owing to their instability under non-equilibrium conditions [8]. Therefore, the development of interatomic potentials that are not based on the shell model is highly important for the intensive study of TiO2 properties. Zhang et al. [9] established a Ti–O fitting database and obtained potential parameters by using a parallel genetic algorithm; in addition, they developed binary modified embedded-atom method potentials of Ti–O on the basis of spline interpolation. Meanwhile, Ti–O interactions at the dilution limit of O were simulated, and the transitional potential barrier of gap levels and gap sites of O and O in stacking faults of Ti prisms. Regardless, this potential function is mainly used for MD simulation of the diffusion of pure Ti or

Corresponding author. . E-mail address: [email protected] (X. Huang).

https://doi.org/10.1016/j.physb.2019.08.034 Received 1 July 2019; Received in revised form 8 August 2019; Accepted 20 August 2019 Available online 20 August 2019 0921-4526/ © 2019 Elsevier B.V. All rights reserved.

Physica B: Condensed Matter 574 (2019) 311657

Y. Tian, et al.

O near Ti defects. Hammond et al. [10] proposed Buckingham potentials applicable for the Fe–Y–Ti–O system and investigated the aggregation of Y–Ti–O solute atoms in the Fe matrix by using the Monte Carlo method. In the current study, Y, Ti, and O atoms were regarded as charge-neutral substances as they focused on the embedment of the Y–Ti–O system into the Fe matrix. Thus, the long-range Coulomb force did not have to be calculated, and the effects of ionic charges on calculations of speed were neglected in the development of the potential function. Bandura et al. [11] proposed a study of structures of water molecules and counterions near the TiO2/water interface by using the first-principles calculation method. The Buckingham potential parameters of Ti–O, Ti–Ti, O–O, and Ti–O (H2) were obtained by fitting in the Ti–O–H system. Although it can describe Ti–O–H interactions on the TiO2 surface in different coordination environments, this potential is mainly applicable for MD simulation of the water/TiO2 interface. In summary, interatomic potentials of Ti - Ti, Ti–O, and O–O were involved in previous studies; however, interatomic potentials proposed in previous reports cannot be directly applied in MD simulation of rutile TiO2 crystals because of their differences in the cutoff radius, number, and type of atoms. Meanwhile, no studies have been conducted on interatomic potentials for MD simulation of rutile TiO2. Therefore, we developed a novel interatomic potential based on Buckingham pairwise potentials by combining the Coulomb potentials of rutile TiO2. We then calculated several key physical parameters of rutile TiO2 crystals. The results based on the proposed potentials were compared with the experimental values and those obtained by first-principles calculation to verify the accuracy of the proposed interatomic potentials.

Table 1 Parameters for TiO2 pairwise potentials. TiO2

Ti–Ti Ti–O O–O Charge(C) Cut off, rc(Å)

In this study, pairwise interatomic potentials of Ti–Ti, Ti–O, and O–O were developed using the Buckingham pairwise potentials [12] by combining the Coulomb potentials of rutile TiO2 crystals. Eq. (1) shows the potential function used. The Buckingham potential is the most commonly used potential for describing atomic interactions in ionic materials. It mainly contains the repulsive term and the attractive term. The former degrades exponentially as the radial distance r increases, and the latter depends on r6. The attractive term and the repulsive term are mainly attributed to dipole–dipole interactions and Pauli repulsion, respectively. Once atoms become close to each other, the repulsive term dominates, and the corresponding atom shells overlap. The repulsive term and the attractive term are key factors in obtaining the physical properties of materials at equilibrium. Moreover, interatomic distances tend to be extremely small under high pressure, resulting in significantly high requirements for the repulsive term. If the overlapping effect of ion shells are weak, the trend of repulsive interactions can be considered consistent with the exponential decay of the electron density function. Thus, the repulsive term in the Buckingham potentials can effectively simulate the repulsions between the electron nucleus and the shell [13]: (1)

UijCoulomb

(2)

= qi qj / rij

UijBuckingham = Aij exp(−

rij ρij

) − Cij / rij6

ρ(Å)

C(eV Å6)

9440561.00 20473.35 8745.77 qTi 1.29 15

0.154 0.194 0.234 qO −0.65

27.99 12.05 29.18

experimental values of the lattice constant, bulk modulus, shear modulus, and elastic constant of rutile TiO2 were used as input values for fitting potential parameters. Essentially, fitting empirical potentials is a process of identifying the minimum of the target function, which is the sum of squares of differences between the input and the calculated parameters. In GULP, the target function was minimized using the Newton–Raphson algorithm. The Buckingham potential model (Eq. (3)) includes three independent parameters—Aij, ρij, and Cij. Table 1 summarizes the fitted potential parameters. Once the Buckingham potentials were obtained by fitting, MD simulations were conducted using GULP to calculate key parameters (e.g., bulk modulus, shear modulus, and elastic constant) of TiO2. The entire simulation process was under the NPT ensemble, with the temperature and pressure of the system controlled by the Nose–Hoover thermostat [15]. The model used was a 2 × 2 × 3 TiO2 supercell containing 72 atoms. The simulation duration was 60 ps, with a time step of 1 fs and relaxation time of 30 ps. To verify the validity of the Buckingham potentials, we employed first-principles calculations because of the absence of some key experimental data of TiO2. The firstprinciples calculations of the lattice constant and the elastic constant of rutile TiO2 crystals at 0 K were achieved using the CASTEP software [16], and the interactions of ions and valence electrons were described using ultrasoft pseudopotentials. The cutoff energy of 390 eV was determined based on convergence. The exchange–correlation function was described by GGA–PBE functional analysis [17]. The unit cell used to calculate the elastic constant and the lattice constant exhibited a tetragonal structure containing six atoms. Integration into the Brillouin zone was achieved using the 5 × 5 × 8 Monkhorst–Pack special K-point method, and the overall energy difference at self-consistent convergence was set to 10-6 eV/atom.

2. Theoretical method and potential function model

Uij = UijCoulomb + UijBuckingham

A(eV)

3. Results and discussion 3.1. Lattice constant As shown in Fig. 1, the rutile TiO2 crystals are tetragonal with a space group number of P42/mnm. Ti4+ ions occupy the end points and body center of the cubic primary cell, while O2− ions occupy the centers of orthogonal triangles with Ti4+ ions as vertices. The lattice parameters at different temperatures are obtained to verify the potential function. The lattice parameter at 0 K cannot be measured experimentally; thus, we calculated the theoretical value of the lattice parameter at 0 K by using the DFT method (see Table 2). As shown in Table 2, the lattice parameter at 300 K (GULP column), calculated using the developed potentials, is highly consistent with that measured experimentally [19]. Meanwhile, the lattice parameter at 0 K (GULP column), calculated using the developed potentials, is highly consistent with that calculated using the DFT method. Therefore, the proposed interatomic potentials exhibit good accuracy. Phonon dispersion curves are significant to verify the accuracy of empirical potentials. Essentially, phonon dispersion curves reflect interatomic bonding in solids, mechanical stability, and thermodynamic properties. Empirical potentials with reasonable accuracy tend to repeat

(3)

where Aij, ρij, and Cij are the Buckingham potential parameters; qi and qj are ionic charges; and rij is the distance between atom i and atom j. The pseudopotential electronic configuration of Ti and O are 3s23p63d24s2 and 2s22p4, respectively. The Born effective charges (BECs) calculated using density functional theory (DFT) were regarded as ionic charges, and the BECs of Ti and O in rutile TiO2 were +1.29 e and −0.65 e, respectively. The cutoff radius was 15 Å. The Buckingham potential parameters were fitted using the GULP software [14], and MD simulations were executed accordingly. The 2

Physica B: Condensed Matter 574 (2019) 311657

Y. Tian, et al.

while the other 15 phonon dispersion curves are optical modes. Acoustic modes reflect the vibrations of the center of mass, and optical modes reflect the relative vibration of particles in the primary cell. Moreover, lattice waves in crystals are equivalent to harmonic oscillators under harmonic approximation, and the energy of lattice waves is equivalent to that of the corresponding harmonic oscillator. The space group number of rutile TiO2 crystals was P42/mnm, and the optical phonon at Point Γ in the Brillouin zone can be described using the following irreducible equation [20]:

Γopt = A1g + A2g + A2u + 2B1u + B1g + B2g + Eg + 3Eu where g and u refer to Raman and IR active, respectively, and E refers to the degenerate mode. A2u and Eu are polar and can split to longitudinaloriented (LO) and transverse-oriented (TO) optical phonons with different frequencies owing to the relevance of their macroscopic electric fields to LO phonons. Short-range interatomic forces cause anisotropy, and the TO or LO modes of A2u and Eu vary in frequency have different frequencies at Point Γ. In LO lattice vibrations with symmetry of A2u and Eu, atoms shifted to the c axis in directions parallel and perpendicular to it, respectively, both of which were activated by IR light. A2u (LO) and Eu (TO) are effective for light polarizations parallel and perpendicular to the c axis, respectively. If the IR deactivation mode is nonpolar, B2g, B1g, and A1g are in Raman active modes. A2g (non-polar) and B1u are the IR mode and the Raman inactive mode, respectively. Table 3 summarizes the phonon spectrum frequencies at Point Γ calculated using the proposed method and the experimental values [21] reported elsewhere. As observed, these values are highly consistent. Table 4 summarizes the optical- and acoustic-mode phonon frequencies at the high-symmetry points X, M, Z, A, R in the Brillouin zone and experimental values. As observed, optical-mode phonon frequencies at points M and Z are consistent with those obtained experimentally, whereas optical- and acoustic-mode phonon frequencies at point X are slightly higher than the experimental values. This finding may be attributed to the fact that potential parameters obtained in this study were based on the static structural properties of materials, and vibrational frequencies calculated based on potential parameters deviated from experimental values.

Fig. 1. Structure of the primary cell of rutile TiO2 crystals. Table 2 Lattice parameters of TiO2. Lattice parameters(Å)

GULP

DFT

Experiment

TiO2 a(0K)

4.52

-

b(0K)

4.52

c(0K)

3.10

a(300K) b(300K) c(300K)

4.53 4.53 3.09

4.66 4.55 [18] 4.66 4.55 [18] 2.97 2.92 [18] -

4.59 [19] 4.59 [19] 3.00 [19]

R

4

25

15 1

Frequency(THz)

20

3.2. Elastic constant and mechanical parameters

11 5 7 14 2

15

10

13

10

8 9

5

3

0

[ , , ]

0.5

[ , ,0]

0.0

[1/2, ,0]

[ ,0,0]

0.5 0.0/0.5 wave vector

The elastic properties of solid materials reflect their strains caused by variations in external conditions [24]. The strength of atomic bonds, the equation of state, and the phonon spectra significantly affect the elastic properties of solid materials. The elastic properties are reflected by the elastic constant, which determines the capacity of crystals to external forces. The elastic constant determines the responses of crystals to external excitations, which are described by bulk modulus (B),

[0,0, ]

0.0

Table 3 Comparison of computed mode frequencies (in THz) at Point Γ with experimental data for rutile TiO2.

[0, ,1/2]

0.5/0.0 0.5

Fig. 2. Phonon dispersion curves of rutile TiO2.

phonon dispersion curves obtained experimentally. Fig. 2 illustrates the phonon dispersion curves of rutile TiO2 crystals at 0 GPa obtained by MD simulation based on the developed potential function. The TiO2 crystal theoretically has 18 dispersion relations as the primary cell of TiO2 contains six atoms. Thus, the TiO2 crystal is expected to have 18 phonon dispersion curves. This conclusion is consistent with the results of this study with regard to the number of phonon dispersion curves. The 18 phonon dispersion curves in Fig. 2 consist of three acoustic modes and 15 optical modes. The lowest three curves at Point Γ are acoustic modes because they are within the acoustic frequency range, 3

Mode

This work

Experiment [21,22]

Experiment [23]

A1g(1) A2g(2) B1g(3) B2g(4) Eg(5) B(1) 1u(6) B(2) 1u(7) A2u(TO)(8) E1 u(TO)(9) E2 u(TO)(10) E3 u(TO)(11) A2u(LO)(12) E1 u(LO)(13) E2 u(LO)(14) E3 u(LO)(15)

20.72 15.05 3.70 25.84 18.78 – 17.11 9.37 7.32 14.30 20.14 – 13.46 15.01 24.87

18.36 – 4.29 24.78 13.41 – – 5.01 5.49 11.64 15.00 24.33 11.19 13.74 24.18

18.30 – 4.25 24.72 13.34 3.39 12.18 5.18 5.66 – 14.81 – 11.23 12.85 25.24

Physica B: Condensed Matter 574 (2019) 311657

Y. Tian, et al.

Table 4 Calculated (cal.) and experimental (exp.) [21] mode frequencies (THz) at highsymmetry points X, M, Z, A, R in the Brillouin zone. X calc.

X exp.

Optical modes 4.974 3.084 – 8.071 7.918 9.033 9.972 9.130 14.337 11.756 17.971 13.528 20.381 – 25.320 25.38 Acoustic modes 4.974 3.084 5.898 5.821

M calc.

M exp.

Z calc.

Z exp.

A calc.

R calc.

– – 8.030 8.838 – 14.412 21.126 23.688

3.058 7.811 9.398 9.442 – 13.552 – 23.26

6.883 9.266 9.896 13.738 13.742 19.311 20.011 22.943

– 9.040 9.720 12.312 12.397 – – –

7.711 8.964 12.530 13.319 17.184 18.541 20.013 23.635

8.970 9.224 9.717 15.162 16.468 18.038 20.376 23.517

2.989 2.989

2.936 3.058

– –

3.697 12.397

3.525 3.525

5.496 –

Table 6 Comparison of the mechanical properties of TiO2 in the current study and in the other results.

Shear modulus, G(GPa)

GULP(this work)

DFT(this work)

Other results

92.57

122.99

a)

69.05 [25] 113.98 [26] a) 223.01 [25] b) 230 ± 20 [27] b) 212 [28] a) 187.76 [25] – – a) 0.364 [25] a) 0.251 [29] a) 8722 [25] b)

Bulk modulus, B(GPa)

225.16

210.33

Young modulus, E(GPa)

105.74(Ex) 105.74(Ey) 370.46(Ez) 0.316

182.93(Ex) 182.93(Ey) 375.80(Ez) 0.366

Longitudinal sound velocity, vp(m/s)

9539.94

9440.50

Transverse sound velocity, vs(m/s)

5253.59

5411.41

Debye temperature(K)

795.46

815.16

Poisson Ratios(σ)

a)

shear modulus (G), Young's modulus (E), the Poisson's ratio (σ), and sound speed (V). Therefore, nine components of the elastic constant of TiO2 crystals were calculated by MD simulation based on the developed potential function and the first-principles calculation method based on the density functional theory (DFT) (see Table 5). Some elastic constants in the 6 × 6 elastic constant matrix were equivalent to each other (e.g., C12 = C22, C13 = C23, C44 = C55), which can be attributed to the symmetry of rutile TiO2 crystals. Therefore, the elastic constants of rutile TiO2 were reduced from nine to six, including C11, C12, C13, C33, C44, and C66. Moreover, elastic constants calculated based on the developed interatomic potentials were highly consistent with those calculated using the first-principles method. Meanwhile, the elastic constants calculated based on the developed interatomic potentials and the first-principles method were highly consistent with the values reported in previous studies [25], verifying the validity and good accuracy of the developed potential function. As shown in Table 6, the mechanical properties of rutile TiO2, including bulk modulus (B), shear modulus (G), Young's modulus (E), the Poisson's ratio (σ), and sound speed (V), can be obtained using the Voigt–Reuss–Hill approximation theory [32]. The coefficients obtained using the Voigt approximation theory and the Reuss approximation theory are the upper and lower limits of the elastic constant, respectively; the average of these two values is defined as the theoretical elastic constant. According to the Voigt approximation theory and the Reuss approximation theory, the bulk modulus and shear modulus of tetragonal rutile TiO2 crystals can be calculated using the following:

GV =

1 (C11 + C33 − 2C13 + 6C44 + 3C66) 15

(4)

BV =

1 (C11 + 3C12 + 4C13 + C33) 9

(5)

15 GR = 4S11 − 8S13 + 4S33 + 6S44 + 3S66 BR =

a)

a) b)

(8)

B = (BV + BR)/2

(9)

5163.98 [30] 815 [35] 778 [31]

The theoretical results (DFT). The experiment data.

using Eqs. (4)–(9). The bulk modulus (225.16 GPa) calculated based on the developed interatomic potential was highly consistent with that calculated using the first-principles method (210.33 GPa) and that reported by Hu et al. (223.01 GPa) [25]. The longitudinal sound velocity (vp) and transverse sound velocity (vs) of rutile TiO2 crystals can be calculated based on their shear modulus and bulk modulus [33]:

vp = (

B + 4G /3 1/2 ) ρ

(10)

vs = (

G 1/2 ) ρ

(11)

where ρ refers to the density of rutile TiO2, B refers to the bulk modulus, and G refers to the shear modulus. Table 6 shows the longitudinal sound velocity (vp) and transverse sound velocity (vs) obtained using Eqs. (10) and (11). Owing to the inconsistency of shear modulus, the transverse wave (5253.59 m/s) and longitudinal wave (9539.94 m/s) speeds calculated based on the developed interatomic potentials exceeded 4083 and 8722 m/s, respectively [25]. Subsequently, the Debye temperature (θD) of rutile TiO2 crystals can be calculated using the Debye model [34]:

θD =

h 3n NA ρ 1/3 [ ( )] vm k 4π M

(12)

where h is the Planck constant, k is the Boltzmann constant, NA is the Avogadro constant, M is the molecular mass, and n is the number of atoms in the molecule. The average sound velocity can be determined using the following:

(7)

G = (GV + GR)/2

a)

9896.29 [30] 4083 [25]

b)

(6)

1 S11 + 3S12 + 4S13 + S33

a)

1 2 1 vm = [ ( 3 + 3 )]−1/3 3 vt vl

(13)

As shown in Table 6, the Debye temperature calculated using the developed interatomic potentials (795.46 K) is consistent with that

Table 6 presents the bulk modulus and shear modulus calculated

Table 5 Elastic constants Cij (in GPa) of TiO2 calculated by MD simulation (using GULP) and DFT calculation and compared with the results of Hu [25].

GULP DFT Hu's [25]

C11

C12

C13

C22

C23

C33

C44

C55

C66

271.84 275.50 289.51

221.82 148.34 183.95

149.69 155.39 165.60

271.84 275.50 –

149.69 155.39 –

448.88 485.27 463.13

97.30 133.83 111.40

97.30 133.83 –

222.23 198.80 229.45

4

Physica B: Condensed Matter 574 (2019) 311657

Y. Tian, et al. 9

820

(a)

1.8

7 1.5

800

790

RDF

6

RDF

Unit cell volume(Å3)

810

300 K 2000 K

(b)

8

5

1.2

0.9

4

0.6

780

0.3

3

8.00

8.25

8.50

8.75

9.00

9.25

9.50

9.75

10.00

r(Å)

770

2

1900K

760 500

1000

1500

1

2000K 2000

2500

0

3000

1

2

3

4

5

6

7

8

9

10

r(Å)

Temperature(K)

Fig. 3. (a) Supercell volume vs temperature; (b) Radial distribution function of TiO2 at 300 and 2000 K.

calculated using the first-principles method (815 K) [35], indicating that the developed interatomic potentials can effectively describe the elastic constant and mechanical properties of rutile TiO2.

MD simulation were consistent with the experimental values, demonstrating that the developed interatomic potentials are applicable in the study of the mechanical properties of rutile TiO2 crystals. (2) Phonon dispersion curves of rutile TiO2 crystals were calculated based on the developed interatomic potentials. The phonon spectrum frequencies at Point Γ were highly consistent with the experimental values and those calculated using the first-principles method. Moreover, the melting point was verified by MD simulation based on the developed interatomic potentials. The melting point of rutile TiO2 was predicted to be 1950 K by using the heating-until-melting method. This value was very close to the experimental value, demonstrating the good validity of the developed interatomic potentials. (3) The semi-empirical method used for fitting of potentials is limited by several intrinsic limitations and inaccuracy. To enhance the accuracy of the potential parameters, future studies are intended to include an in-depth consideration of the effects of various factors in obtaining potential parameters with improved accuracy. This study serves as a reference to our future studies and other efforts in the development of interatomic potentials.

3.3. Verification of melting point Melting point verification is essential in evaluating the accuracy of interatomic potentials. Once the external temperature exceeds the melting point, both solid and liquid phases exist in crystals and transform to each other until dynamic equilibrium is reached. Before melting, crystal structures remain regular, and atoms vibrate near the lattice only. Once the external temperature exceeds the melting point, atoms in crystals move randomly, and the ordered lattice structures are damaged. Crystal volumes are exposed to significant changes before and after melting. Thus, the melting point of rutile TiO2 was obtained using the heating-until-melting method in this work [36]. A rutile TiO2 supercell (2 × 2 × 3) was first established, and system temperature was adjusted under the NPT ensemble by MD simulation (GULP) to reach equilibrium. The relaxation time was 1 ns, and the time step was 0.05 ps. As shown in Fig. 3(a), the temperature corresponding to the discontinuity point in the supercell volume vs. temperature curve during melting is identified as the melting point of TiO2. The volume of the rutile TiO2 supercell abruptly changes at 1900–2000 K; thus, the melting point of rutile TiO2 ranges from 1900 K to 2000 K. However, the precise melting point cannot be determined based on the information in Fig. 3 but rather, on the radial distribution functions of rutile TiO2 crystals at different temperatures. Fig. 3(b) illustrates the radial distribution functions of rutile TiO2 at different temperatures. A radial distribution function can describe the structure, characterize the degree of structural disorder, and relate the atomic structure to the macroscopic characteristics. As observed in Fig. 3(b), peaks at 300 K appear narrow and high, demonstrating that crystal structures are regular and rutile TiO2 is in the solid state. As the temperature increases, peaks tend to be broad and low, demonstrating a disordered crystal structure attributable to ordered thermal vibrations. All peaks, except for the first one, disappeared gradually at 2000 K. Therefore, the melting point of rutile TiO2 is 1950 K, which is close to the experimental value (2143 K) [37].

Conflicts of interest The authors declare that there is no conflict of interest regarding the publication of this paper. Acknowledgments This work is financially supported by NSAF (Grant no. U1530140) and by Chongqing Science and Technology Commission (Nos.CSTC2017JCYJAX0357 and CSTC2016JCYJA0517), by Science and Technology Project Affiliated to the Education Department of Chongqing Municipality (No. KJ1709224), and by Foundation of National Key Defense Laboratory of Computational Physics (No. 6142A0501020217). References [1] A.V. Bandura, J.D. Kubicki, Derivation of force field parameters for TiO2-H2O systems from ab initio calculations [J], J. Phys. Chem. B 107 (40) (2003) 11072–11081. [2] J. Wang, Y.F. Zhao, T. Wang, et al., Photonic, and photocatalytic behavior of TiO2, mediated by Fe, CO, Ni, N doping and co-doping [J], Physical B Physics of Condensed Matter 478 (2015) 6–11. [3] Alamgir, W. Khan, S. Ahmad, et al., Formation of self-assembled spherical-flower like nanostructures of cobalt doped anatase TiO2 and its optical band-gap [J], Mater. Lett. 133 (2014) 28–31. [4] E.G. Villabona-Leal, J.P. López-Neira, J.A. Pedraza-Avella, et al., Screening of factors influencing the photocatalytic activity of TiO 2 :Ln (Ln = La, Ce, Pr, Nd, Sm, Eu

4. Conclusions (1) Novel interatomic potentials for MD simulation of rutile TiO2 crystals were developed based on the Buckingham potential model. MD simulations based on the developed interatomic potentials were conducted to calculate the physical parameters (e.g., lattice parameter, elastic constant, shear modulus, bulk modulus, wave speed, and the Debye temperature) of rutile TiO2. The values obtained by 5

Physica B: Condensed Matter 574 (2019) 311657

Y. Tian, et al.

2, and Mn F 2[J], Phys. Rev. 154 (2) (1967) 522–526. [22] D.M. Eagles, Polar modes of lattice vibration and polaron coupling constants in rutile (TiO2) [J], J. Phys. Chem. Solids 25 (11) (1964) 1243–1251. [23] J.G. Traylor, H.G. Smith, R.M. Nicklow, et al., Lattice dynamics of rutile [J], Phys. Rev. B Condens. Matter 3 (10) (1971) 3457–3472. [24] I.A. Amenzade, Theory of Elasticity [M], Mir Publishers, 1979. [25] C. Hu, Z. Zeng, C. Kong, et al., High pressure structural instability and thermal properties of rutile TiO2 from first-principles [J], Chin. J. Chem. Phys. 27 (1) (2014) 69–74. [26] D.H. Chung, The voigt-reuss-hill (vrh) approximation and the elastic moduli of polycrystalline ZnO, TiO2 (rutile), and α-Al2O3 [J], J. Appl. Phys. 39 (6) (1968) 2777–2782. [27] L. Gerward, J.S. Olsen, Post-rutile high-pressure phases in TiO2, J. Appl. Crystallogr. (30) (1997) 259–264. [28] N. Yu, Keisuke Nakayama, Eiichi Takahashi, et al., P-V-Tequation of state of stishovite to the mantle transition zone conditions [J], Phys. Chem. Miner. 31 (10) (2005) 660–670. [29] T. Mahmood, C. Cao, W.S. Khan, et al., Electronic, elastic, optical properties of rutile TiO2 under pressure: a DFT study [J], Phys. B Condens. Matter 407 (6) (2012) 958–965. [30] E. Shojaee, M.R. Mohammadizadeh, First-principles elastic and thermal properties of TiO2: a phonon approach [J], J. Phy. Con. Matt. Ins. Phy. J 22 (1) (2010) 15401-0. [31] T.R. Sandin, P.H. Keesom, Specific heat and paramagnetic susceptibility of stoichiometric and reduced rutile (TiO2) from 0.3 to 20 K [J], Phys. Rev. 177 (3) (1969) 1370–1383. [32] D.H. Chung, W.R. Buessem, The voigt‐reuss‐hill approximation and elastic moduli of polycrystalline MgO, CaF2, β-ZnS, ZnSe, and CdTe [J], J. Appl. Phys. 38 (6) (1967) 2535-0. [33] E. Schreiber, O.L. Anderson, Elastic constants and their measurement [J], J. Appl. Mech. 42 (3) (1973) 747–748. [34] O.L. Anderson, A simplified method for calculating the Debye temperature from elastic constants [J], J. Phys. Chem. Solids 24 (7) (1963) 909–917. [35] R. Sikora, Ab initio study of phonons in the rutile structure of TiO2 [J], J. Phys. Chem. Solids 66 (6) (2005) 1069–1073. [36] S.N. Luo, T.J. Ahrens, T. Çagin, et al., Maximum superheating and undercooling: systematics, molecular dynamics simulations, and dynamic experiments [J], Phys. Rev. B 68 (13) (2003) 399–404. [37] Y. Li, T. Ishigaki, Thermodynamic analysis of nucleation of anatase and rutile from TiO2 melt [J], J. Cryst. Growth 242 (3) (2002) 511–516.

and Gd) in the degradation of dyes [J], Comput. Mater. Sci. 107 (6) (2015) 48–53. [5] Guilian Zhu, Hao Yin, Chongyin Yang, et al., Black titania for superior photocatalytic hydrogen production and photoelectrochemical water splitting [J], ChemCatChem 7 (17) (2015). [6] W.F. Van Gunsteren, X. Daura, N. Hansen, et al., Validation of molecular simulation: an overview of issues [J], Angew. Chem. 57 (4) (2017) 884–902. [7] P.J. Wilde, C.R.A. Catlow, Defects and diffusion in pyrochlore structured oxides [J]. Solid State Ionics, Diffusion & Reactions 112 (3–4) (1998) 173–183. [8] D.S.D. Gunn, N.L. Allan, H. Foxhall, et al., Novel potentials for modelling defect formation and oxygen vacancy migration in Gd2Ti2O7 and Gd2Zr2O7 pyrochlores [J], J. Mater. Chem. 22 (11) (2012) 4675–4680. [9] P. Zhang, D.R. Trinkle, A modified embedded atom method potential for interstitial oxygen in titanium [J], Comput. Mater. Sci. 124 (2016) 204–210. [10] K.D. Hammond, H.J. Voigt, L.A. Marus, et al., Simple pair-wise interactions for hybrid Monte Carlo-molecular dynamics simulations of titania/yttria-doped iron. [J], J. Phys. Condens. Matter 25 (5) (2013) 055402. [11] A.V. Bandura, J.D. Kubicki, Derivation of force field parameters for TiO2-H2O systems from ab initio calculations [J], J. Phys. Chem. B 107 (40) (2003) 11072–11081. [12] R.A. Buckingham, The classical equation of state of gaseous helium, neon and argon [J], Royal Society of London Proceedings 168 (933) (1938) 264–283. [13] S. Yip, Handbook of materials modeling [J], Mater. Today 8 (11) (2005) 345–348. [14] J.D. Gale, A.L. Rohl, The general utility lattice program (GULP) [J], 29 (5) (2003) 291–341. [15] Simone Melchionna, Giovanni Ciccotti, B.L. Holian, Hoover NPT dynamics for systems varying in shape and size [J], Mol. Phys. 78 (3) (1993) 533–544. [16] V. Milman, B. Winkler, Ab initio modeling in crystallography ☆ [J], Int. J. Inorg. Mater. 1 (3–4) (1999) 273–279. [17] J.P. Perdew, K. Burke, M. Ernzerhof, ERRATA: generalized gradient approximation made simple [J], Phys. Rev. Lett. 77 (18) (1998) 3865–3868. [18] B. Montanari, N.M. Harrison, Lattice dynamics of TiO2 rutile: influence of gradient corrections in density functional calculations [J], Chem. Phys. Lett. 364 (5–6) (2002) 528–534. [19] J.K. Burdett, T. Hughbanks, G.J. Miller, et al., Structural-electronic relationships in inorganic solids: powder neutron diffraction studies of the rutile and anatase polymorphs of titanium dioxide at 15 and 295 K [J], ChemInform 18 (37) (1987) 3639–3646. [20] J.G. Traylor, H.G. Smith, R.M. Nicklow, et al., Lattice dynamics of rutile [J], Phys. Rev. B Condens. Matter 3 (10) (1971) 3457–3472. [21] S.P.S. Porto, P.A. Fleury, T.C. Damen, Raman spectra of Ti O 2, Mg F 2, Zn F 2, Fe F

6