Molecular dynamics simulation of thermodynamic and thermal transport properties of strontium titanate with improved potential parameters

Computational Materials Science 60 (2012) 123–129

Contents lists available at SciVerse ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Molecular dynamics simulation of thermodynamic and thermal transport properties of strontium titanate with improved potential parameters Wen Fong Goh ⇑, Tiem Leong Yoon, Sohail Aziz Khan School of Physics, Universiti Sains Malaysia, 11800 Penang, Malaysia

a r t i c l e

i n f o

Article history: Received 13 October 2011 Received in revised form 12 March 2012 Accepted 14 March 2012 Available online 7 April 2012 Keywords: Molecular dynamics simulation Strontium titanate Thermal expansion Isothermal compressibility Heat capacity Thermal conductivity

a b s t r a c t A molecular dynamics simulation has been performed to investigate the thermal expansivity, isothermal compressibility, heat capacity and thermal conductivity of strontium titanate. The potential model captures the ionic and covalent characteristics of strontium titanate well. The parameters of the model were derived by fitting against the experimental lattice parameters. With these fitted parameters, we then evaluated the variations of lattice, thermal expansion coefficient, isothermal compressibility, heat capacity and thermal conductivity as a function of temperature from room temperature up to 2000 K, and pressure from ambient pressure up to 20.3 GPa. The thermal conductivity calculations were performed using non-equilibrium molecular dynamics method, and corrections for finite size effect were made. The simulation results are in good agreement with the experimental data and the theory. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Strontium titanate SrTiO3 is one of the most common ceramic materials. It is paraelectric above 105 K [1], and has an ideal perovskite structure with lattice constant of 3.9051 Å [2] and space  at room temperature. This material finds applicagroup of Pm3m tions among other thing in electronics and electrical insulations. Various researches had shown that by introducing oxide defects in strontium titanate or doped with alkaline earth metals or transition metals, its electrical, thermal and thermoelectric properties can be modified. For example, niobium-doped strontium titanate can be electrically conductive [3], lanthanum-doped strontium titanate is a good thermoelectric material [4,5] and cerium-doped strontium titanate can be used for fuel cells applications [6]. Knowledge of the reduced structural properties of the material is essentials. Many of the applications of the material operate at high temperatures. Performing experimental researches are costly or difficult due to practical limitations, especially at very high temperature and high pressure. As feasible alternative, computer simulation is routinely carried out to understand the thermo-physical properties and thermal stability of the material under investigation. Computer simulation also allows us to investigate the physical properties of a material at extreme temperature and pressure without actually carrying out the experiments in the lab. In this paper, we present a molecular dynamics simulation on strontium titanate to investigate its thermodynamic and thermal ⇑ Corresponding author. E-mail address: [email protected] (W.F. Goh). 0927-0256/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2012.03.027

transport properties, i.e. isobaric thermal expansivity, isothermal compressibility, heat capacity and thermal conductivity. Comparison is then made against experimental findings to justify the correctness of the numerical results. Deriving an accurate and reliable force field or interacting potential is crucial in molecular dynamics simulation. To our knowledge, there exist two well known sets of potential parameters for strontium titanate, i.e. published by Katsumata et al. [7] and Seetawan et al. [8]. However, (i) the set of potential parameters published by Seetawan et al. [8] has an unphysical attraction between Sr and O at close range, which causes the simulation to be unstable (number of particles in the simulation are not preserved) and leads to failure of reproducing the experimental lattice parameters at high temperatures, (ii) the fact that we are able to reproduce Katsumata et al. [7] lattice parameters using their potential parameters even to elevated temperatures, strengthen the case, that Seetawan et al. [8] reproduction at high temperatures has been over looked and is wrongly presented. This is obvious since using his [8] potential parameters we were not able to reproduce his result. None of these inadvertences have been pointed out. Due to these reasons, we are motivated to address the issues related to the potential parameters of strontium titanate. 2. Simulation method 2.1. Potential model Deriving an accurate and reliable force field is always a preliminary step in molecular dynamics simulation. In the present

124

W.F. Goh et al. / Computational Materials Science 60 (2012) 123–129

Table 1 Potential parameters of the potential model. The values in the table are taken from Katsumata et al. [7] except ci, which are obtained in the present work. The last two rows show the parameters for the Morse-type potential.

ri (Å)

qi (Å)

Sr Ti O

2.0 2.2 1.4

1.846 1.385 1.850

0.13 0.10 0.18

0 0 25.04

Dij (kcal mol1)

xij (Å1)

r0ij (Å)

26.0

2.0

Ti–O

Kcal mol

qi (e)

Present work

200 0

repulsion

2

4

Seetawan et al.'s potential 8

6

8

10

r

200

USr

O

Ion

  1 3 1 ci Kcal2 Å mol 2

400

1.6

400 600

unphysical attraction at close range

800

work, we used pairwise interaction potential of the Ida-type repulsion [9] and Morse-type covalent bonding:

qi qj Uðrij Þ ¼ þ f0 ðqij Þ exp r ij

Fig. 1. Seetawan et al.’s [8] potential for Sr–O pair showing unphysical attraction at close range. Our improved potential avoids the unphysical attraction.

!

rij  rij cij  6 þ Dij fexp½2xij ðrij qij r ij

 r0ij Þ  2 exp½xij ðrij  r0ij Þg

ð1Þ

where i and j refer to the type of ion in a unit cell: i, j 2 {Sr, Ti, O}. The first term represents Coulombic interaction, where qi and qj are the effective charges of the ith and jth intermolecular ions, and rij is the distance between them. The second term represents short range repulsion, where force constant, f0 = 1.0 kcal/Å mol1, rij = ri + rj and qij = qi + qj. The parameters ri and qi cater for the ionic radii and softness respectively. The third term typifies Van der Waals attraction, where cij = cicj. The last term represents Morse-type asymmetric stretching and compression of covalent bond for the intermolecular ions, where Dij is the depth of the Morse potential well, xij is the potential steepness and r0ij is the equilibrium distance for ions i and j. In our calculations, the parameters in the Morse-type potential assume non-zero values only for the Ti-O pair. The potential parameters used in the present work are shown in Table 1. The potential parameterization is discussed in detail in Section 3.1. 2.2. Simulation details Simulations were performed on a cell of size corresponds to 8  8  8 unit cell, which consists of 2560 ions. Initial positions and velocities were generated according to experimental crystal structure and Maxwell velocity distribution. Verlet algorithm with a time step of 1.0 fs was used to solve the equation of motion. Simulations in canonical ensemble (i.e. NVT-constant) and isothermal-isobaric ensemble (i.e. NPT-constant) were performed using Nose–Hoover algorithm with a time constant of 0.1 ps and 1.0 ps for thermostat and barostat respectively. A global cut-off distance of 11.0 Å was used with long-range Van der Waals tail correction being added. Ewald summation with a precision of 104 was used for long-range Coulombic interaction. Periodic boundary condition was imposed. Our simulations were performed using the opensource package Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [10]. All systems were equilibrated for 20,000 time steps. After equilibration, the results were extracted by performing time averaging of 200,000 time steps. 3. Results and discussion 3.1. Potential parameter analysis Deriving an accurate and reliable force field is important in molecular dynamics simulation, because the quality of the force field is the single most significant factor that determines the correctness of the simulation results. Any force field must reproduce the experimental results, such as lattice parameters. For strontium

titanate, there exist two well known sets of potential parameters, i.e. published by Katsumata et al. [7] and Seetawan et al. [8]. However, using the set of potential parameters published by Seetawan et al. [8], we were unable to reproduce their lattice parameters at temperatures above 1200 K, though we were able to reproduce their lattice parameters as indicated in their paper at lower temperatures. We then carefully analyze the set of potential parameters published by Seetawan et al. [8]. Our analysis shows that the interacting potential of Sr–O in their potential parameters [8] has an unphysical attraction at close range (Fig. 1), which causes the simulation to be unstable and results in particle lost. This unphysical closed range attraction is the reason for the failure of the potential parameters published by Seetawan et al. [8], i.e. to being unable to preserve the number of particles in the simulation. This inadvertence has not been pointed out. Our further analysis shows that, at lower temperatures, this effect is not be observable, because other components of the potential parameters keep the ions far apart that the chances of getting the atoms falling into the potential well are negligible and makes the unphysical close-range attraction to have no effect. Conversely, at higher temperatures, the particles have higher energy that the chances of getting Sr and O particles overlapping are higher. It is not clear how Seetawan et al. [8] were able to obtain the value of lattice parameters at temperature beyond 1200 K, while at the same time having an unphysical interaction in their potential. The fact that we are able to reproduce Katsumata et al.1 [7] lattice parameters using their potential parameters even to elevated temperatures strengthen the case that Seetawan et al. [8] reproduction at high temperatures has been over looked and is wrongly presented. Furthermore, according to Katsumata et al. [7], their potential parameters are not yet fully optimized to reproduce the physical properties other than lattice parameters, thermal expansion and compressibility, e.g. heat capacity and thermal conductivity. This implies that their potential parameters need to be improved in order to calculate these two properties. Motivated by such inadvertences, we decided to parameterize some of the potential parameters. We adopted the values of the parameters qi, ri, qi, DTiO, xTiO, r0TiO of Katsumata et al. [7] but ignore their Van der Waals attraction parameters ci. The Van der Waals attraction parameter for anion-anion interaction, i.e. cOO, is obtained by fitting our numerical results against the experimental lattice parameters and assuming cSr and cTi both equal to zero. The latter assumption is based on the consideration that 1 Initially when using the set of potential parameters published by Katsumata et al. [7], we failed to reproduce the lattice parameters as indicated in their paper. Our analysis shows that this is so because their Van der Waals attraction parameters ci (i = Sr, Ti, O) are too large. After that we contacted them and they mentioned that the values are a factor of 1000 short and should be revised to ci  1000. Now we are able to reproduce their lattice parameters.

125

W.F. Goh et al. / Computational Materials Science 60 (2012) 123–129 Table 2 Structural parameters of strontium titanate at 298 K and 1 atm. a, b, c are the lattice constants, A, B, C are the angles and . is the density. Present work

Experiment [2]

a (Å) b (Å) c (Å) A (°) B (°) C (°) . (kg/m3) Molar volume (105 m3 mol1)

3.90516 3.90509 3.90516 89.9992 90.0002 90.0005 5116.09 3.58641

3.9051 3.9051 3.9051 90 90 90 5116.7 3.5863

Experiment 11

Lattice parameter

Structural parameters

Present work

3.98

Seetawan et al. 8

3.96

3.94

3.92

3.90 500

Linear thermal expansion coefficient a can be derived from the variation of lattice parameter with temperature (Fig. 2). Both the experimental [12,15–17] and our result for linear thermal expansion coefficients between 298 K and 2000 K are shown in Fig. 7. The simulated and experimental linear thermal expansion coefficient at room temperature and ambient pressure are 1.010  105 K1 (present work), 1.077  105 K1 [12], 0.863  105K1 [15] and 0.940  105 K1 [16] respectively. The experimental values of the linear thermal expansion coefficient found in the literature are inconsistent. Ligny et al. [12] deemed that it remains almost constant from room temperature up to 1800 K. Our calculation suggests a diminutively growth from 1.010  105 K1 at room temperature to 1.376  105 K1 at 2000 K. 3.4. Isothermal compressibility Isothermal compressibility bT can be derived from the variation of volume with pressure. The isothermal compressibility versus pressure up to 20.3 GPa at room temperature is shown in Fig. 8. There is only one single measured experimental value of

2000

3.85

Molar volume, V 10 5 m3 mole

Present work

3.80

Experiment 12 Katsumata et al. 7

3.75 3.70 3.65 3.60 3.55

500

1000

1500

2000

Temperature, T K Fig. 3. Variation of molar volume as a function of temperature.

3.60

Molar volume, V 10 5 m3 mole

3.3. Thermal expansion

1500

Fig. 2. Variation of lattice parameter as a function of temperature.

3.2. Crystal structure and structural parameters Based on the potential model (Eq. (1)) and the fitted parameters in Table 1, the physical properties of strontium titanate have been computed. The simulated and experimental [2] structural parameters of strontium titanate at room temperature and ambient pressure are shown in Table 2 for comparison. The simulation results match the experimental values very well. Figs. 2 and 3 show the variations of simulated lattice parameter and simulated molar volume respectively, as a function of temperature, together with the experimental data [11,12] and other simulation results [7,8]. The lattice parameter increases monotonically from 3.9051 Å at room temperature to 3.9825 Å at 2000 K, whereas the molar volume increases monotonically from 3.586  105 m3 mol1 at room temperature to 3.804  105 m3 mol1 at 2000 K. We also evaluated the variation of molar volume with respect to pressure at fixed temperature. The variation of molar volume as a function of pressure at 298 K, 381 K and 467 K from our simulation and experiments [13,14] are shown in Figs. 4–6. At 298 K, the molar volume decreases from 3.586  105 m3 mol1 at ambient pressure to 3.247  105 m3 mol1 at 20.3 GPa. The lattice parameters and molar volumes from the simulation are in excellent agreement with the experimental data. Compare with other simulations, our results are better than that of Seetawan et al. [8] and agree to that of Katsumata et al. [7].

1000

Temperature, T K

the Van der Waals attractions are only significant between anionanion pairs. The potential parameters obtained in this way are shown in Table 1.

Present work Experiment 13 Experiment 14 Katsumata et al. 7

3.55 3.50 3.45 3.40 3.35 3.30 3.25 3.20

0

5

10

15

20

Pressure, P GPa Fig. 4. Variation of molar volume as a function of pressure at 298 K.

5.75  1012 Pa1 at 298 K and 1 atm [1], which is to be contrasted with our calculated value of 5.80  1012 Pa1. Our simulation result suggests that it decreases from 5.80  1012 Pa1 at ambient pressure to 3.66  1012 Pa1 at 20.3 GPa. The thermal variation of isothermal compressibility between 298 K to 2000 K is shown in Fig. 9. Our calculation predicts that it will increase from 5.80  1012 Pa1 at room temperature to 8.25  1012 Pa1 at 2000 K. 3.5. Heat capacity The heat capacity at constant pressure, cP is defined as:

cP ¼

  @H @T P

ð2Þ

126

W.F. Goh et al. / Computational Materials Science 60 (2012) 123–129

6.0 1

Pa

Present work Experiment 13

3.5

3.4

3.3

0

5

10

15

Present work

5.5

12

3.6

Isothermal compressibility, T 10

Molar volume, V 10 5 m3 mole

3.7

Experiment 1

5.0

4.5

4.0

3.5

20

0

5

10

Pressure, P GPa

15

20

Pressure, P GPa

Fig. 5. Variation of molar volume as a function of pressure at 381 K.

Fig. 8. Isothermal compressibility as a function of pressure.

1

8.5

3.5 3.4 3.3

0

5

10

15

Present work

8.0

12

Pa

Presen twork Experiment 13

3.6

Experiment 1

Isothermal compressibility, T 10

Molar volume, V 10 5 m3 mole

3.7

7.5 7.0 6.5 6.0

20

500

Pressure, P GPa

1000

1500

2000

Temperature, T K

Fig. 6. Variation of molar volume as a function of pressure at 467 K.

Fig. 9. Isothermal compressibility as a function of temperature.

1.4

140

Heat capacity at constant pressure, Cp J mol K

Linear thermal expansion coefficient, 10 5 K 1

1.6

1.2 1.0

Present work Experiment 12 Experiment 15 Experiment 16 Experiment 17

0.8 0.6 0.4

120 100

Present work Experiment 18

80

Experiment 12 Experiment 19

60

Seetawan et al. 8 500

1000

1500

2000

Temperature, T K

40

800

1000

1200

1400

1600

1800

2000

Temperature, T K

Fig. 7. Linear thermal expansion coefficients as a function of temperature.

Fig. 10. Heat capacity at constant pressure as a function of temperature.

where H = E + PV is the enthalpy of the system, P is the external pressure, V is the volume and E is the system’s internal energy. In order to obtain cP at temperature T, values of H at five temperatures: T, T ± DT and T ± 2DT, where DT  10 K, were evaluated using NPT-constant simulation. The variation of the enthalpy with temperature is linear and hence a straight line can be fitted through these data. The slope of the line yields cP. Fig. 10 shows simulated heat capacities at constant pressure as a function of temperature together with the experimental data [18,12,19]. cP shown here are above Debye temperature (694 K [20]), where simulation is expected to work more efficiently in view of the classical basis of the program. The simulation result of Seetawan et al. [8] shows anomalies at temperature range around 1200 K to 1600 K. They claimed that this is due to the

strontium titanate undergoing structural change to a tetragonal structure and then going back to the cubic structure at higher temperature. However, no evidence can confirm these anomalies. Our simulation result suggests no such transformation but a linear increase that is confirmed by the experiment. 3.6. Thermal conductivity Thermal conductivity calculations were performed using nonequilibrium molecular dynamics method as described by MüllerPlathe [21]. This method is implemented in LAMMPS [10]. The heat flux is imposed through kinetic energy exchange and the temperature gradient is the system’s response. Each of the simulation boxes of size 3  3  60 unit cell up to 3  3  220 unit cell was

127

W.F. Goh et al. / Computational Materials Science 60 (2012) 123–129

W mK

12 10

Thermal conductivity,

8 6 4 2 0

0

50

100

150

200

250

300

Simulation length, L z unit length Fig. 13. Thermal conductivity as a function of simulation length at room temperature. The dashed line is the limit of the fitting function.

W mK

Fig. 11. Instantaneous heat flux as a function of time for 3  3  60 unit cell at 298 K.

14 Present work Müller Plathe Present work Green Kubo Experiment 5

12

Thermal conductivity,

10

Experiment 23 1473K,1h

8

Experiment 23 1673K,2h

6 4 2 0

500

1000

1500

2000

Temperature, T K Fig. 14. Thermal conductivity as a function of temperature.

Fig. 12. Instantaneous thermal conductivity as a function of time for 3  3  60 unit cell at 298 K.

subdivided into 20 bins along z direction. The kinetic energy of the hottest atom in the 1st bin and the coolest atom in the 11th bin (middle bin) were exchanged at fixed interval, so that the 11th bin became the hottest bin and the 1st bin became the coolest bin. This induced a temperature gradient along the z direction. Thus, the thermal conductivity, j can be calculated by:

j¼

Jz dT=dz

ð3Þ

where Jz is the heat flux and dT/dz is the temperature gradient along the z direction. We found out that the kinetic energy exchange after the range of 90–210 time steps yield fairly linear temperature gradient for all time steps, hence a dwell time of 150 time steps has been chosen to compute the thermal conductivity. Instantaneous heat flux and thermal conductivity as a function of time for 3  3  60 unit cell at 298 K are shown in Figs. 11 and 12. Obviously, the instantaneous thermal conductivity converges at 0.3 ns or 300,000 time steps. Therefore, the simulations were equilibrated for 300,000 time steps, and then the kinetic energy transferred and temperature gradient were accumulated for 3,000,000 time steps. Phonon scatterings at boundaries of heat source (11th bin) and heat sink (1st bin) give rise to finite size effect is non-negligible in molecular dynamics, if the simulation length is not significantly

larger than the phonon mean free path. Worth mentioning that the phonon mean free path of strontium titanate is around 15 Å at room temperature [22]. The thermal conductivity calculations were carried out on cells of size 3  3  f unit cell, where f = 60, 80, 100, 120, 180, 200, 220. Fig. 13 shows simulated thermal conductivity as a function of simulation length Lz at room temperature. It shows a trend of thermal conductivities converging at large simulation length. A function of the form.

f ðLz Þ ¼ ð1  eLz =s Þ

ð4Þ

where  and s are fitting parameters, is fitted through these points as shown in Fig. 13. The saturation limit of the fitting function gives thermal conductivity in the bulk state, which is at infinite length. Thus, this method eliminates the finite size effect due to phonon scattering at the boundaries. Fig. 14 shows simulated thermal conductivity as a function of temperature, together with the experimental data [5,23]. Note that the two experimental values of the thermal conductivity are inconsistent. From kinetic theory, thermal conductivity as a function of temperature should follows 1/Tg at high temperature, where g has a value between 1 and 2. The experimental values from Muta et al. [5] fulfilled this trend, but the experimental values from Ito and Matsuda [23] showed almost constant values with increasing temperature. Our simulation result obtained using Müller-Plathe method shows declination from 11.04 W/mK at room temperature to 3.485 W/mK at 2000 K, and follows the trend of 1/Tg with increasing temperature as predicted by the theory. Our result is in good agreement with the experimental values obtained from Muta et al. [5]. It is not clear how Seetawan et al.’s [8] result of thermal conductivity obtained by them agrees so well with the

128

W.F. Goh et al. / Computational Materials Science 60 (2012) 123–129

W mK

80 average

60

xx

Thermal conductivity,

yy

40

zz

20

0

500

20

1000

1500

2000

Time, t ps

Fig. 15. Instantaneous thermal conductivities as a function of time for 8  8  8 unit cell at 298 K. javerage = (jxx + jyy + jzz)/3.

experimental values, whereas the lattice parameters obtained from their potential parameters do not agree with the experimental values. We have carried out the thermal conductivity calculations using alternative method, i.e. Green–Kubo method, which is implemented in LAMMPS [10]. The Green–Kubo method starts by calculating the instantaneous heat flux ~ JðtÞ from the phase space configuration at time t given by the following equation:

" # P~ ! 1 P ~ ~ ~ J ðtÞ ¼ ei v i þ ðf ij  v j Þxij V i i
ð5Þ

where ei is the energy of atom i; ~ f ij is the force on atom i due to neighboring atom j; ~ v i is the velocity of atom i, and ~xij is the distance between atom i and j. The thermal conductivity tensor jij is then obtained by integrating the ensemble average of the time correlation function of the heat flux Ji(t):

jij ¼

V kB T 2

Z

high temperatures. Due to unreliable potential parameters found in the literature, a set of potential parameters for the potential function has been re-derived by fitting the calculated lattice parameters against the experimental values using LAMMPS [10]. With the set of improved potential parameters, we have used the LAMMPS [10] molecular dynamics package to calculate molar volumes and lattice parameters in the temperature range of 298–2000 K, and the pressure range of 0.1 MPa to 20.3 GPa. The variations of lattice parameters and molar volumes as a function of temperature and pressure are in good correspondence with the experimental data. Thermal expansion coefficient and isothermal compressibility as a function of temperature at fixed pressure, and isothermal compressibility as a function of pressure at fixed temperature have been calculated. These thermodynamic properties are compatible with the known experimental data. Heat capacity at constant pressure as a function of temperature has been evaluated. Our result did not show anomalies, which are introduced in the simulation result of Seetawan et al. [8]. Our values correspond to the experimental data. Using non-equilibrium molecular dynamics method, thermal conductivity as a function of temperature has been calculated. The result with corrections made for finite size effect matches the trend predicted from the kinetic theory and is in good agreement with the experimental values. Therefore, it can be concluded that our improved potential parameters are successful in obtaining the physical, thermodynamic and thermal transport properties of strontium titanate. Acknowledgement The authors acknowledge the support from the USM Fellowship. References

1

0

hJ i ð0ÞJ j ðtÞidt

ð6Þ

where i, j = x, y, z, V is the volume, T is the temperature and kB is the Boltzmann’s constant. Details can be found elsewhere [24]. Selection of correlation time is important in Green–Kubo calculation. It should be long enough to obtain the decay of the correlation function, but not too long that the noise added to the correlation signal contaminates the intrinsic value of the stress correlation function [25]. In this calculation, correlation time of 20 ps has been used. Fig. 15 shows instantaneous thermal conductivities as a function of time for 8  8  8 unit cell at 298 K. Obviously, the thermal conductivities converge at 1000 ps or 1,000,000 time steps. One advantage of the Green–Kubo method is the extraction of thermal conductivities along the three directions simultaneously. Note that jxx, jyy and jzz are almost equal at high temperature. This implies the isotropic characteristic of strontium titanate. Simulation result obtained using Green–Kubo method is also shown in Fig. 14. Our result agrees reasonably well with the results from different experiments. 4. Conclusion Molecular dynamics simulation has been carried out on strontium titanate in order to determine the correct parameters for the potential and therefore its thermodynamic and thermal transport properties, i.e. thermal expansion, isothermal compressibility, heat capacity and thermal conductivity. We have pointed out the set of potential parameters published by Seetawan et al. [8] has an unphysical attraction between Sr and O at close range, which causes the simulation to being unable to preserve the number of particles in the simulation and leads to the failure of reproducing the experimental lattice parameters at

[1] R.O. Bell, G. Rupprecht, Phys. Rev. 129 (1963) 90–94, http://dx.doi.org/ 10.1103/PhysRev.129.90. . [2] R.W.G. Wyckoff, Crystal Structures, John Wiley and Sons, New York Sons, New York, 1964. [3] P. Blennow, A. Hagen, K.K. Hansen, L.R. Wallenberg, M. Mogensen, Solid State Ion. 179 (35–36) (2008) 2047–2058, http://dx.doi.org/10.1016/j.ssi.2008.06.023. . [4] H. Muta, K. Kurosaki, S. Yamanaka, J. Alloys Compd. 368 (1–2) (2004) 22–24, http://dx.doi.org/10.1016/j.jallcom.2003.07.016. . [5] H. Muta, K. Kurosaki, S. Yamanaka, J. Alloys Compd. 392 (1–2) (2005) 306–309, http://dx.doi.org/10.1016/j.jallcom.2004.09.005. . [6] D.J. Cumming, J.A. Kilner, S. Skinner, J. Mater. Chem. 21 (13) (2011) 5021–5026, http://dx.doi.org/10.1039/C0JM03680C. . [7] T. Katsumata, Y. Inaguma, M. Itoh, K. Kawamura, Solid State Ion. 108 (1–4) (1998) 175–178, http://dx.doi.org/10.1016/S0167-2738(98)00036-8. . [8] T. Seetawan, G. Wong-ud-dee, C. Thanachayanont, V. Amornkitbumrung, Chinese Phys. Lett. 27 (2) (2010) 026501. . [9] Y. Ida, Phys. Earth Planet. Inter. 13 (2) (1976) 97–104, http://dx.doi.org/ 10.1016/0031-9201(76)90074-1. . [10] S. Plimpton, J. Comput. Phys. 117 (1) (1995) 1–19, http://dx.doi.org/ 10.1006/jcph.1995.1039. . [11] E.N. Bunting, G.R. Shelton, A.S. Creamer, J. Am. Ceram. Soc. 30 (4) (1947) 114– 125, http://dx.doi.org/10.1111/j.1151-2916.1947.tb18879.x. . [12] D. de Ligny, P. Richet, Phys. Rev. B 53 (1996) 3013–3022, http://dx.doi.org/ 10.1103/PhysRevB.53.3013. . [13] M. Guennou, P. Bouvier, J. Kreisel, D. Machon, Phys. Rev. B 81 (2010) 054115, http://dx.doi.org/10.1103/PhysRevB.81.054115. . [14] L. Edwards, R. Lynch, J. Phys. Chem. Solids 31 (3) (1970) 573–574, http:// dx.doi.org/10.1016/0022-3697(70)90098-3. . [15] K. Itoh, K. Ochiai, H. Kawaguchi, C. Moriyoshia, E. Nakamura, Ferroelectrics 159 (1994) 85–90. [16] F.W. Lytle, J. Appl. Phys. 35 (7) (1964) 2212–2215, http://dx.doi.org/10.1063/ 1.1702820. .

W.F. Goh et al. / Computational Materials Science 60 (2012) 123–129 [17] S. Devanarayanan, P.S. Narayanan, Ind. J. Pure Appl. Phys. 6 (1968) 714–716. [18] K. Jacob, G. Rajitha, J. Chem. Thermodynam. 43 (1) (2011) 51–57, http:// dx.doi.org/10.1016/j.jct.2010.08.011. . [19] J.P. Coughlin, R.L. Orr, J. Am. Chem. Soc. 75 (3) (1953) 530–531, http:// dx.doi.org/10.1021/ja01099a006. arXiv:http://pubs.acs.org/doi/pdf/10.1021/ ja01099a006, . [20] H. Ledbetter, M. Lei, S. Kim, Phase Trans. 23 (1) (1990) 61–70, http://dx.doi.org/ 10.1080/01411599008241819. arXiv:http://www.tandfonline.com/doi/pdf/10. 1080/01411599008241819, . [21] F. Müller-Plathe, J. Chem. Phys. 106 (14) (1997) 6082–6085, http://dx.doi.org/ 10.1063/1.473271. .

129

[22] A.M. Hofmeister, J. Appl. Phys. 107 (10) (2010) 103532, http://dx.doi.org/ 10.1063/1.3371815. . [23] M. Ito, T. Matsuda, J. Alloys Compd. 477 (1–2) (2009) 473–477, http://dx.doi.org/ 10.1016/j.jallcom.2008.10.031. . [24] J.W. Lawson, M.S. Daw, J. Charles, W. Bauschlicher, J. Appl. Phys. 110 (8) (2011) 083507, http://dx.doi.org/10.1063/1.3647754. . [25] D. Nevins, F.J. Spera, Mol. Simulat. 33 (15) (2007) 1261–1266, http://dx.doi.org/ 10.1080/08927020701675622. arXiv:http://www.tandfonline.com/doi/pdf/10. 1080/08927020701675622, .