CC BY
14
УДК 538.95
Modelling the Structure of Nanoparticle-embedding Matrices: Molecular Dynamics in Li2B4O7
Alain Marbeuf* Janis Kliava
CPMOH, UMR 5798 Universite Bordeaux I - CNRS 351 Cours de la Liberation, 33405 Talence Cedex,
France
Received 06.01.2010, received in revised form 10.02.2010, accepted 20.02.2010 A new model describing interatomic and angular interactions, taking into account periodic properties in borate-type solid phases, is presented and applied to Li2B4 O7 through simulations at temperatures ranging from 0 K to the melting point and in the pressure range 0 to 10000 MPa. Simulation reproduces quite well cell lengths, atomic positions and distances in boron-oxygen polyhedrons and the polar nature of the crystal structure. An order-disorder type ferro-paraelectric transition of the second kind is found to occur at a Curie point TC ~ 839 K, corresponding to jumping of Li atoms between two lattice sites. By increasing or decreasing the pressure, the total energy and the crystal properties for simulations performed at 300 K show a shoulder at pt ~ 5000 MPa, implying the existence of a reversible second-order phase transition. The cell volume below pt follows a Murnaghan law with the bulk modulus B0 = 15.6 GPa and its first derivative Bo' = 4.31 (at ambient pressure). In contrast to the low-pressure phase where threefold and fourfold boron atoms coexist, in the high-pressure phase all borons are fourfold-coordinated. The present approach can be directly applied to modelling the structure of nanosized systems.
Keywords: lithium tetraborate, molecular dynamics, crystal structure, order-disorder ferroelectric transition, pressure-induced phase transition.
Introduction
As transparent magnets, borate glasses containing nanoparticles represent a new class of nanomaterials, promising for technical applications [1]. The parent crystalline phase Li2B4O7 (LTBc) has received much attention because of piezoelectric and nonlinear optical properties (second harmonic generation, SHG) related to its polar crystal structure [2-11].
In this work, a new model describing interatomic and angular interactions, taking into account periodic properties in borate-type solid phases, is presented and applied to LTBc. Starting from the point of view based on DFT calculations with coordination-dependent boron charges [12-14],a parameter set is chosen from structural considerations and tested through various simulations at temperatures ranging from 0 K to the melting point. The description of LTBc is achieved through simulations in the pressure range 0 to 104 MPa. A special attention is paid to structural and dielectric properties of the different crystal phases.
* a.marbeuf@cpmoh.u-bordeaux1.fr © Siberian Federal University. All rights reserved
1. Molecular Dynamics and Interaction Model in Li2B4O7 System
The potential model is based on two-particle and three-particle potentials. The bond interactions are described by a Born-Mayer repulsion term added to the Coulomb potential:
Vij qi qj e /rij + bij exp( rij /Pij )• (1)
where rij is the interatomic distance between ions i and j, e is the elementary charge, qi and qj are the relative effective charges of i and j ions, bij, pij are the repulsive parameters in this potential model. The infinite sum of all Vj is evaluated by the Ewald sum technique [12] implying real-space calculations inside a sphere of a cutoff radius (ca 12 A) and reciprocal-space calculations outside.
We follow Maslyuk et al. [13-15] who found a good agreement between DFT results and experiments (in crystal structure, elastic properties, optical and phonon spectra) by assuming different effective charges for threefold (B3) and fourfold (B4) boron atoms, qB3 and qB4, but invariant for different ternary oxides of the [(1 — x) B2O3 + x Li2O] system and temperature-independent. The x-invariance of the boron charges implies a variation of both qLi and qo with the concentration ratio r = [B4] / ([B3] + [B4]), satisfying the charge neutrality and adjusted in order to obey to the characteristics of phonon dispersion in ternary crystals [14]. As a consequence, a correlation appears between qLi and qo charges, as shown in Fig. 1, except for Li3BO3 (x = 0.75) where the absence of direct connectivity between B3 results in stronger Li-O interactions.
•0.76'
-0.80 -0.84'
cf
-0.88 -0.92-0.96-1---■----1 1-------1----1 1----1----1-----1----10.58----------------------------------------------------0.60-0.62-0.64-0.66-0.68
Oli
Fig. 1. Correlation between lithium and oxygen effective relative charges for different ternary oxides [(1 — x) B2O3 + xLi2O]; r is the fraction of fourfold boron. The full squares are fitted
values by Maslyuk et al. [14]; the full curve is a guide for the eye
For angular interactions involving an ion i, we choose a truncated harmonic potential:
Vijk = 1 Aijk(-dijk — j)2 exp[—(rij8 + r-ik8)/pijk8] (2)
where Aijk represents the strength of the angular interactions, $ijk is the angle formed by the
ions i, j, k, $jk is the equilibrium angle and pijk is the truncating parameter. In LTBc we distinguish four angular interactions, resulting from:
(i) the presence of (B3) and (B4) involving $jk = 120° and $jk = 109.47° respectively;
(ii) the oxygen angular connectivity which is different between two tetrahedrons (^j ~ 135°) and between a tetrahedron and a triangle ($jk ~ 120°).
The entire parameter set used in the simulations described below is collected in Table 1. All values come from Maslyuk et al. [15] except for the boo coefficient which is adjusted through a MD equilibration at 0 K with respect to the crystallographic data of [16], extrapolated from the 80-300 K range to 0 K. Our boo value (9.1 x 105 kJ/mol) appears less repulsive than that given in [15] (1.37 x 107 kJ/mol).
Table 1. Parameters of the potential model used in MD simulations of Li2B4O7 crystal
Elements Atomic charges in e units (from [14])
Li qLi = 0.666
B3 qB3 = 1.072
B4 qB4 = 1.223
O qo = —0.846
Interactions Parameters and their values (from [ [14]], except boo)
Li - O bLio = 137631 kJ/mol PLio = 0.236 A
B3 - O bB3o = 51846 kJ/mol PB3o = 0.198 A
B4 - O bB4o = 112389 kJ/mol PB4o = 0.198 A
O - O boo = 910000 kJ/mol Poo = 0.183 A
O - B3 - O AoB3o = 0.3580 kJ/mol/rad2 $oB3o = = 120.07° poB3o = 2.4 A
O - B4 - O AoB4o = 0.1587 kJ/mol/rad2 $oB4o = = 109.41° poB4o = 2.4 A
B3 - O - B4 AB3oB4 = 0.00758 kJ/mol/rad2 $B3oB4 = = 120.67° PB3oB4 = 2.4 A
B4 - O - B4 Ab4ob4 = 0.1149 kJ/mol/rad2 $B4oB4 = 135.68° pB4oB4 = 2.4 AA
Simulation boxes of 13000 atoms (5 x 5 x 5 cells with 8 formula units) were prepared using experimental crystallographic data [17,18] based on the tetragonal space group I41cd (Table 2) with c-axis as the polar direction. The MD-simulations were performed with DL_POLY [19] in the Berendsen isothermal and isobaric (NpT or N<rT isostress ensembles) under periodic boundary conditions, with a time step of 0.5 fs. The Berendsen thermostat was used with a time constant of 0.05 ps. The barostat time constant of 5 ps was necessary to ensure stability of the simulations. The structure was relaxed by energy minimization coupled to a heat bath at 0 K in the NvT-ensemble. Then a series of equilibration simulations was performed at atmospheric pressure in NpT conditions, by temperature steps of 20 K below 300 K and of 50 K above 300 K. In this series the experimental room-temperature c/a ratio from [17] was used. Above 400 K,
simulations were also made with the c/a ratio deduced from [16]. Besides, in order to verify that the tetragonal symmetry is conserved, some simulations were made in the NaT ensemble (where b/a and c/a are free to vary). In the NpT simulations, fractional coordinates of oxygen atoms O4 were fixed accordingly to symmetry rules of the space group. The simulations at pressures up to 10 GPa were performed at room temperature.
2. Molecular Dynamics and Interaction Model in Li2B4O7 System
2.1. Structural Results
The NpT and NaT simulations consistently demonstrate the stability of the quadratic cell. Table 2 compares the lattice parameters a, c obtained from the simulations at 300 K with the experimental data. Fig. 2 shows the experimental (a, [17]) and the simulated structure (b) viewed
Table 2. Lattice parameters at 300 K: comparison between experiments [17,18] and simulations
Experimental Simulation
a (A) 9.479 ± 0.003 [17] 9.470 ± 0.006 [18] 9.586(1)
c (A) 10.280 ± 0.004 [17] 10.28 ± 0.01 [18] 10.431(1)
along the c-axis. Firstly, we can see that the simulation reproduces quite well the cell lengths with a small departure (Да/а = 1.12 %, Де/c = 1.47%). Secondly, the calculated distances between B3 and O atoms, d(B3 — O) = 1.279 A, are shorter, while those between B4 and O, d(B4 — O) =1.500 A, are slightly longer in comparison with the mean experimental values [17]: d(B3 — O) =1.373Aand d(B4 — O) = 1.476 A, respectively. Nevertheless, the whole crystal structure is quite well reproduced.
2.2. Thermodynamical Results
Temperature dependences of the total energy (Etot) and the cell volume (u) are shown, respectively, in Figs. 3-4. Both curves can be fitted with second-order polynomials, see the legends of Figs. 3-4.
Thus, no phase transition is found below 700 K, in agreement with X-ray diffractometric [16] and IR spectroscopic studies [2]. Burak and Moroz [20] have claimed having observed an isostructural phase transition at 235 K: no such transitions have been detected in our simulations.
Finally, the melting point is found between 950 and 1050 K. This uncertainty comes from the choice of the starting c/a ratio in the simulations. Nevertheless, these values are not too much below the thermal determinations of Shastry and Hummel [21] (see Fig. 5).
Fig. 2. Experimental (a, in accordance with the crystal structure determined by Krogh-Moe [17]) and simulated (b) room-temperature crystal structure viewed along the c-axis
Fig. 3. Temperature dependence of the total energy in NpT simulations ▲: c/a ratio from [17]; V: c/a from [16]. The full line is the corresponding fit: Etot = —7049.40 + 0.2992 T + 2.20 x 10-5 T2
2.3. Non-linear Properties
For each simulation run, the dielectric polarization vector (P) has been calculated. Its component along the c-axis (Pz) is found to decrease with increasing temperature above 200 K
Table 3. Mean fractional coordinate of Li atoms along the c-axis and its spatial fluctuation, as a function of temperature below and above the phase transition
T (K) ZLi AzLi
300 0.8724 0.0013
400 0.8744 0.0013
800 0.8773 0.0022
850 0.8780 0.0030
910 0.8769 0.0041
Fig. 4. Temperature dependence of the cell volume in NpT simulations. ▲: c/a ratio from [17]; ♦ : experiment [17]; V: c/a from [16]; Y: experiment [16]. The full line is the corresponding fit: u = 913.7 + 0.140 T + 3.6 x 10-5 T2
(Fig. 5), according to the empirical law
Pz 2 = 165 + 0.13 (839 — T). (3)
One can see that Pz does not vanish above T = 839 K (its "remnant" value is Pz « 13cm-2). Meanwhile, Px and Py also have non-zero values which are temperature-independent, so that above T = 839 K all three components of P become comparable. This behaviour is an artifact due to border effects. Indeed, in contrast to the other simulations made (see Section II), P is calculated within the simulation box without applying periodic boundary conditions.
As shown in Fig. 5, Pz2 behaves similarly to the second-harmonic generation coefficient (S33) measured by Furusawa et al. [7]. These authors found a linear temperature dependence of S33 above room temperature and, extrapolating this dependence, suggested the existence of a second-order transition at To « 863 K. Because no soft mode associated with a displacive-type
phase transition was observed by Raman spectroscopy in the temperature range 300-900 K [11], Furusawa et al. concluded that this second-order phase transition could be of order-disorder type. The present simulations show that fluctuations of the fractional coordinate of Li along the c-axis drastically increase in the temperature range 850-950 K (see Table 3). This behaviour can be explained assuming that in this range the Li atoms are jumping between two different sites in a crystallographic supercell. So, our simulations confirm the possibility of a second-order ferroelectric ^ paraelectric phase transition in LTBc, the Curie point TC « 839 K given by eq. (3) corresponding to an order-disorder transition.
Fig. 5. Square of the calculated Pz-component of the dielectric polarization compared with the second-harmonic generation coefficient S33 as a function of temperature. ▲ and V have the same meaning as in the preceding figures ; : S33 values from [7]; : the melting point Tfus = 1190
K [21]. The full line is the fit, see eq. (3); the dashed line shows an extrapolation of high-temperature S33 values
3. Simulations at Room Temperature
3.1. Low-pressure Phase
Pressure dependences of the total energy and the cell volume are shown, respectively, in Figs. 6-7. By increasing the pressure p, all these quantities unambiguously show a leap occurring at pt ~ 5000 MPa, and by decreasing the pressure one passes through exactly the same values. These smooth ‘discontinuities’ and the absence of pressure hysteresis imply the existence of a second-order pressure-induced phase transition, and one can see that this transition is quite reversible.
In the low-pressure phase (called thereafter I - [LTBc]), Etot can be fitted with the following equation:
Etot = -6957.83 - 0.0021p + 1.27 x 10-4p2 (4)
Particularly interesting is the fact that in the I - phase, which is stable at ambient pressure,
Fig. 6. Pressure dependence of the total energy at room temperature. A and △: NpT simulations with increasing and decreasing the pressure, respectively. The full and dash-dotted lines are the corresponding fits, eqs. (4) and (7)
the calculated u(p) dependence can be fitted by the following equation (see Fig. 7):
\
P\
Q 2000 4000 6000 8000 10000
p (MPa)
Fig. 7. Pressure dependence of the cell volume with increasing (a) and decreasing (△) the pressure in NpT simulations at room temperature. ♦: experiment [17]; ▼: experiment [16]. The full and dash-dotted lines are the corresponding fits, eqs. (5) and (8)
u = 958.5 1 +
3610
-0.232
(5)
P
This equation of state is of the Murnaghan type [22]:
— 1
( B0' \ Bo'
u(p) = uo ( 1 + B=0P) (6)
where, in our case, at atmospheric pressure the bulk modulus is Bo = 15.6 GPa and its first derivative is B0' = 4.31.
3.2. High-pressure Phase
In contrast with the low-pressure phase, in the high-pressure phase (called thereafter as II [LTBc]) Etot and u depend linearly on the pressure:
Etot = -6930.2 + 0.0043 p (7)
and
u = 740.3 - 9.37 x 10-3p. (8)
The II-polymorph structure is illustrated in Fig. 8: the mean atomic fractional coordinate values show that the symmetry operations of the space group I41cd are conserved, in agreement with the second-order character of the I ^ II phase transition. This pressure-induced phase transition
• O 0 04 ®B3 • B4 O Li
Fig. 8. Calculated high-pressure crystal structure of Li2B4O7 at room-temperature viewed along the c-axis
obeys to the usual increase of the boron coordination: all B-atoms are in B4 conformation (r = [B4] / ([B3] + [B4]) = 1), as in the II - B2O3 high-pressure form studied by Prewitt and Shannon [23]. The II - [LTBc] phase involves two kinds of boron-oxygen tetrahedrons (see Fig. 9). The former B3 atoms are at the centre of distorted tetrahedrons with an average B-O distance 1.390 A. The former B4 atoms remain at the centre of distorted tetrahedrons, but with larger average B-O distance (1.512 A), nearly equal to its value in the I - [LTBc] polymorph (1.497 A). The oxygen O4 associated with the short B - O distance is coordinated by only two B4 atoms, whereas the other oxygens are coordinated by three borons. The resulting boron-oxygen network has a tetrahedral linkage similar to that in II - B2O3 and produces cages for Li atoms. Note that the B3 to B4 transformation has been obtained without any change in the force field between the low- and high-pressure phases.
Fig. 9. Boron-oxygen tetrahedral linkage in the high-pressure form of Li2B4O7. The atom symbols are the same as in Fig. 2 and Fig. 8
By analogy with ab initio calculations on II - B2O3 [24,25], one can assume that II - [LTBc] is more ionic than I - [LTBc]. Therefore, the Coulomb contribution to the total energy is the driving force for the transformation from I - [LTBc] to II - [LTBc] structures.
4. Conclusions
A model describing interatomic and angular interactions with boron charge depending on the coordination has been used in molecular dynamics of the polar Li2B4O7 crystal phase as a function of temperature under constant pressure and as a function of pressure at constant temperature. The simulations performed under atmospheric pressure provide a good description of the crystal structure and confirm the order-disorder ferroelectric transition of the second kind below the melting point, previously inferred from the second-harmonic generation experiments. The calculations performed at room temperature predict a pressure-induced phase transition preserving the crystal symmetry. The corresponding mechanism is the same as for the B2O3 oxide where all threefold borons become fourfold under pressure.
The consistent approach in molecular dynamics is a powerful tool to predict crystal structure and to assure a continuity of physical properties between crystalline, glassy and nanometric phases.
We acknowledge computational facilities provided by the Universite Bordeaux I with POle M3PEC Mesocentre Regional (http://www.m3pec.u-bordeaux1.fr) supercomputers.
References
[1] J. Kliava, I. Edelman, O. Ivanova, R. Ivantsov, O. Bayukov, E. Petrakovskaja, V. Zaikovskiy, I. Bruckental, Y. Yeshurun, S. Stepanov, J. Appl. Phys., 104(2008), 103917-1-11.
[2] N.D. Zhigadlo, M. Zhang, E.K. Salje, J.Phys.: Condens. Matter, 13(2001), 6551-6561.
[3] A.V. Vdovin, V.N. Moiseenko, V.S. Gorelik, Ya.V. Burak, Phys. Solid State, 43(2001), 16481654.
[4] G.L. Paul, W. Taylor, J. Phys. C: Solid State Phys, 15(1982), 1753-1764.
[5] S-I. Furusawa, S. Tange, Y. Ishibashi, K. Miwa, J. Phys. Soc. Jpn., 59(1990), 1825-1830.
[6] A.E. Elalaoui, A. Maillard, M.D. Fontana, J. Phys.: Condens. Matter, 17(2005), 7441-7454.
[7] S-I. Furusawa, O. Chikagawa, S. Tange, T.Ishidate, H. Orihara, Y. Ishibashi, K.Miwa, J. Phys. Soc. Jpn.,60(1991), 2691-2693.
[8] H.A. Sideck, G.A. Saunders, B.J. James J. Phys. Chem. Sol., 51(1990), 457-465.
[9] T.Y. Kwon, J.J. Ju, H.K. Kim, J.W. Cha, J.N. Kim, M. Cha, S.I. Yun, Mater. Lett, 27(1996), 317-321.
10] H.R.Jung, B.M.Jin, J.W. Cha, J.N. Kim, Mater. Lett, 30(1997), 41-45.
11] I. Martynyuk-Lototska, O. Mys, V. Adamiv, Y.Y. Burak, R.Vlokh, Ukr. J. Phys. Opt., 3(2003), 264-266.
12] P. Ewald, Ann. Phys., 64(1921), 253.
13] V.V. Maslyuk, M.M.Islam, T.Bredow, Phys. Rev, B72(2005), 125101-1-9.
14] V.V. Maslyuk, T.Bredow, H.Pfnur, Eur. Phys. J., B41(2004), 281-287.
15] V.V. Maslyuk, T.Bredow, H.Pfnur, Eur. Phys. J., B42(2004), 461-466.
16] M.D. Mathiews, A.K.Tyagi, P.N.Moorthy, Thermochim. Acta, 320(1998), 89-95.
17] J. Krogh-Moe, Acta Cryst, 15(1962), 190-193; B24(1968), 179-181.
18] S.F.Radaev, L.A. Muradyan, L.F. Malakhova, Y.V. Burak, V.I. Simonov, Kristallografiya, 34(1989,) 1400.
19] W. Smith, T.R. Forester, I.T. Todorov, DL_POLY Code, STFC Daresbury Laboratory, Daresbury, Warrington WA4 4AD, UK, 2002.
20] Y.V. Burak, I.E.Moroz, Phys. Chem. Glasses, 44(2003), 241-243.
21] B.S. R. Sastry, F.A. Hummel, J. Am. Ceram. Soc., 41(1958), 7-17.
22] F.D. Murnaghan, Proc. National Acad. Sci. USA, 30(1944), 244-247.
23] C.T. Prewitt, R.D. Shannon, Acta Cryst., B, 24(1968), 869-874.
24] A. Takada, C.R. A.Catlow, J.S. Lin, G.D. Price, M.H.Lee, V. Milman, M.C. Payne, Phys. Rev., B, 51(1995), 1447-1455.
[25] A. Takada, C.R.A. Catlow, G.D. Price, C.L. Hayward, Phys. Chem. Minerals, 24(1997), 423431.
Моделирование структуры матриц, включающих наночастицы: молекулярная динамика в Li2B4O7
Ален Марбеф Янис Клява
Представлена новая модель описания межатомных и угловых взаимодействий в твердофазных кристаллических боратах с учетом периодичности структуры. Модель применена для моделирования свойств Li2B4O7 в диапазоне температур от 0 K до точки плавления и давлений от 0 до 10000 MPa. Моделирование весьма хорошо воспроизводит размеры ячейки, атомные позиции и расстояния в боро-кислородных полиэдрах, а также полярный характер кристаллической структуры. Предсказан сегнето-параэлектрический переход второго рода типа порядок-беспорядок при температуре Кюри (Tc ~ 839K), , обусловленный прыжками атомов лития между двумя позициями в решетке. С увеличением или уменьшением давления полная энергия и свойства кристалла при 300 K показывают наличие ступени при pt ~ 5000 MPa, предполагающей существование обратимого фазового перехода второго рода. При давлениях ниже pt объем ячейки следует закону Мурнагана (Murnaghan) с объемным модулем В0 = 15.6 GPa и его первой производной Во' = 4.31 (при атмосферном давлении). В отличие от фазы низкого давления, в которой сосуществуют трех- и четырехкоординированные атомы бора, в фазе высокого давления весь бор четырехкратно координирован. Настоящий подход непосредственно применим для моделирования структуры наноразмерных систем.
Ключевые слова: тетраборат лития, молекулярная динамика, кристаллическая структура, се-гнетоэлектрический переход типа порядок-беспорядок, индуцированный давлением фазовый переход.
