Interatomic Potential Model for Molecular Dynamics Simulation of Lithium Borate Melts/Glasses

Abstract

An improved interatomic potential model was proposed for molecular dynamics simulations of lithium borate melt/glass systems. Charge of ion was reconsidered and a new composition dependent ionic charge model was suggested. A new three-body potential model controlling B-O-B angles was also proposed. The three-body term functioned to avoid square network ring consisted of B-O bonds, without preventing the change of boron coordination number between three and four. The edge-shared tetrahedra of four-coordination boron observed in the previous simulation were cleared by applying this three-body potential model.

1 Introduction

Borosilicate glasses are the major oxide glasses, used for window or drinking glasses. The major components of the borosilicate glasses are boron oxide, silicon oxide, alkali oxides (A₂O) and alkaline earth oxides. The composition of the oxide glasses is an important factor to give its properties, whereas that makes it difficult to select appropriate composition for each application. There are unexplained relations between the glass molecular structure and the several glass properties. Molecular dynamics (MD) simulation is a prospective method to bring some information about the structure or behavior of multicomponent glasses, but it has not yet become an effective enough tool even for binary glasses. To create the realistic random network oxide glass structure by MD simulation, over a thousand atoms should be necessary. The large numbers of MD steps are also needed for relaxation of the glass structure. These requirements suggest the classical MD simulation is more suitable than the ab initio MD simulation for oxide glasses. In classical MD simulations for borate glasses started by Soules in the 80s, the number of atoms in the MD cell was under one thousand [1,2,3]. To this day, the number of atoms treated in MD simulation continues to increase [4,5,6,7,8,9,10]. In this study, we applied the classical MD to study lithium borate glasses with about 5000 atoms.

Alkali borates: xA₂O-(1−x) B₂O₃ (0 < x < 1) (A = Na, Li, K) are vitrifiable over the x range from zero to around 0.5. The density as an example, some properties of alkali borate glasses are not varied monotonously with x, it is known as "boron anomaly" [11]. The one factor of this phenomenon is thought to be the coexistence of two structural units, i.e., three- or four- coordinated boron. In the pure B₂O₃ glass, only three-coordinated boron (B^III) (that is bonded three oxide ions) exists. Possible chemical reactions when A₂O is added to B₂O₃ are shown in chemical equations (1) and (2), respectively.   

(1)

  

(2)

In the equation (1), a B-O bond is broken and the bridging oxygen (BO, Ø) changes to the non-bridging oxygen (NBO, O⁻). The equation (2) is a deformation of B^III to four-coordinated boron (B^IV). Bray et al. reported the variation of B^IV ratio with x of xA₂O-(1−x) B₂O₃, glasses by ¹¹B NMR results [12]. Their results showed the equation (2) entirely occurs at x < 0.4, and the chemical equation (1) or exchange of B^IV to B^III(equation (3)) are started at x ≥ 0.4.   

(3)

In this work, we discuss the speculation of the ionic charge of the element and the suggestion of an interatomic potential model added three-body terms to two-body terms. Lithium borates xLi₂O-(1−x) B₂O₃ (0 < x < 1) was selected as a first target. We are interested in the lithium borate as a candidate material of the separator used in lithium-ion secondary batteries. The MD simulation results of lithium borates are verified by the reproducibility of the ratio of the B^IV for the all boron, or other structural parameters.

2 Interatomic potential model

2.1 Two-body interatomic potential model

In most MD simulation studies for borate glasses, Coulomb potential with short-range repulsive force potential has been applied as the interatomic potential [1,2,3,4,5,6,7,8,9,10]. An example of another interatomic potential type, by Takada et al. was reported − MD simulation that applied the interatomic potential function as function of the coordination number change of boron [8]. In this work, the two-body interatomic potential model (equation (4)) was applied.   

$\begin{array}{c}{U}_{ij}\left(r_{ij}\right)=\frac{z_i{z}_j{e}^2}{4\pi {\epsilon }_0{r}_{ij}}+f_0\left(b_i+b_j\right)\exp \left(\frac{a_i+a_j-r_{ij}}{b_i+b_j}\right)-\frac{c_i{c}_j}{r_{ij}^6}\\ +D_1{}_{ij}\exp \left(-{\beta }_{1ij}{r}_{ij}\right)+D_2{}_{ij}\exp \left(-{\beta }_{2ij}{r}_{ij}\right)\end{array}$(4)

In order from left, the terms in the right side of equation (4) show Coulomb potential, short-range repulsive force potential, van der Waals force potential, and covalent bond force potential that consist of 4^th and 5^th terms. The letters ε_0,e, and r_ij mean electric constant, elementary charge, and distance of i and j atoms, respectively. Other parameters were settled through try-and-error MD simulation for LiB₃O₅ crystal in this work.

2.2 Charge of ions

Xu et al. [6] reported qualitative reproduction of Bray's R_B^IV results [12] by applying the formal charge to all elements of xNa₂O-(1−x) B₂O₃ system. This result indicates the full-ionic potential model is useful to estimate the structure of glasses, but might not be enough in detail. Because B₂O₃ is an acidic oxide on the contrary, Na₂O and Li₂O are basic oxides. The latter substances are ionic crystals, but the former is assumed to be a covalent material with ionic property. Therefore the B-O bond in the borate glasses should be considered to have both of covalency and ionicity. The 4^th and 5^th terms of the equation (4) are added terms as the covalent bond force potential of B-O bond.

After some trial-and-error simulations for LiB₃O₅ crystal [13], 60% of the formula charge is applied to B ion in this work. For the whole range of xLi₂O-(1−x) B₂O₃ system, the charge of oxide ion, z_O should change as a function of x. The z_O is determined from the charge neutrality of the system. The charges of each ion are set as equation (5).   

$\begin{array}{c}{z}_{\text{Li}}=1.0\\ z_{\text{B}}=1.8\\ z_{\text{O}}\left(x\right)=\left(\left(2-2x\right)z_{\text{B}}+2x\cdot z_{\text{Li}}\right)/(3-2x)\end{array}$(5)

The results of MD simulation applied the partial ionic charge defined by equation (5), however, not successfully reproduced the variation of the glass structure with x. In this work, the charges of boron ion, z_B and z_O are arranged as following;   

$\begin{array}{c}{z}_{\text{Li}}=1.0\\ z_{\text{B}}\left(\alpha \right)=1.8-3\alpha \\ z_{\text{O}}\left(\alpha \right)=-1.2-2\alpha \end{array}$(6)

  

$\alpha =\frac{0.4x}{6-5x}$(7)

Here, α is a function of x as equation (7), derived from the charge neutrality of the system. The z_O, z_B and z_Li values used in this work are listed in Table 1 with numbers of atoms in the MD cell.

Table 2.  Potential parameters

Parameters for the two-body potential terms
element	z		a/nm		b/nm		c/J1/2mol−1/2 nm3
O	−1.2 − 2α		0.1861		0.0151		1.7720
B	1.8 − 3α		0.0723		0.0082		0.0
Li	1.0		0.1110		0.0110		0.1940
pair	D1/MJ mol−1		β1/nm−1		D2/MJ mol−1		β2/nm−1
B-O	120.5582		51.0		−7.3926		22.6
Parameters for the three-body potential terms
relation	f/J	pθ/degree−1		θ0/degree		pr/nm−1		r0/nm
B-O-B	5.0 × 10−4	1.70		100.0		40.00		0.100

Table 1.  Number of atoms, N in MD cell and charges, z

x	Ntotal	NO	NB	NLi	zO	zB	zLi
0.00	5000	3000	2000	―	−1.2000	1.8000	―
0.10	5016	2926	1881	209	−1.2145	1.7782	1.0
0.20	5014	2834	1744	436	−1.2320	1.7520	1.0
0.33	5005	2695	1540	770	−1.2615	1.7077	1.0
0.40	5019	2629	1434	956	−1.2800	1.6800	1.0

2.3 Three-body potential term

In the simulated lithium borate glasses applied the interatomic potential model defined by equations (4) and (6), some undesirable structures included that not given by the Zachariasen's rules for glass formation [14]. Figure 1 shows the undesirable structures schematically. There are two suspicious structures. The first is the edge-shared BO₄-tetrahedron pair, and the second is the three boron coordinated oxide ion (O^III). The former can be described in another way, i.e., 2-membered (or square) ring structure consisting of B-O bonds. In this work, a new three-body potential function was introduced to avoid the undesirable 2-membered ring structure.

Figure 1.

 Schematic view of the undesirable local structure observed in MD results applied 2-body potential model with partial ionic charges for Li₂O-B₂O₃ melt/glasses.

Inoue et al. reported the MD simulation study applied the three-body potential for alkali borate glasses [5]. They introduced interaction functions for O-B-O and B-O-B relations. The function for O-B-O was separated into two patterns corresponding to B^III and B^IV. The three-body potential shifted the center of gravity of B-O-B to the outside of this unit. The simulation successfully reproduced the boroxyl ring (6-membered ring).

In this study, a simple three-body potential function was introduced aimed to exclude the formation of the undesirable structures as shown in Figure 1. As a model that does not prevent the exchange between B^III and B^IV, the three-body potential function with parameter θ_BOB was applied. The three-body potential is shown in equation (8).   

$\begin{array}{c}U(r_{{\text{OB}}_1},r_{\text{OB}}{}_2,{\theta }_{\text{BOB}})\\ =-f\left\{({\theta }_{\text{BOB}}-{\theta }_0)-\frac{1}{p_{\theta }}\ln \left[1+\exp \left[p_{\theta }\right({\theta }_{\text{BOB}}-{\theta }_0\left)\right]\right]\right\}\times \\ \sqrt{k_1\left(r_{{\text{OB}}_1}\right)\cdot k_2\left(r_{{\text{OB}}_2}\right)}\end{array}$(8)

The k₁(r_OB1) and k₂(r_OB2) are shown in equation (9) and (10).   

$k_1\left(r_{{\text{OB}}_1}\right)=\frac{1}{\exp \left[p_r\right(r_{{\text{OB}}_1}-r_0\left)\right]+1}$(9)

  

$k_2\left(r_{{\text{OB}}_2}\right)=\frac{1}{\exp \left[p_r\right(r_{{\text{OB}}_2}-r_0\left)\right]+1}$(10)

Here, θ₀ and r₀ are parameters to control variables of θ_BOB and r_OBn, respectively. Parameters p_θ and p_r are coefficients existing inside of exponential functions. Equations (9) and (10) are extracted from the three-body potential model for water and ice [15]. Figure 2 shows the three-body potential, U. The potential becomes positive when r₀ < r_OBn and θ_BOB < θ₀. Under other conditions, this three-body potential does not have any effect. The parameters were determined through the trial-and-error of MD simulation for Li₂O-B₂O₃ melt system by monitoring the existence of 2-membered ring structure. The obtained parameters are listed in Table 2.

Figure 2.

 Three-body potential, U(θ_BOB, r_OB1, r_OB2)

3 Simulation method

MD simulations for xLi₂O-(1−x) B₂O₃ melts/glasses were performed with about 5000 atoms in a MD cell. The initial positions and velocities of atoms were set randomly, then the MD calculation was started at 2000 K or 3000 K. The pressure was kept at 0.1 MPa. First, the calculation for relaxation of the system was performed that changed the atom positions from meaningless distribution to a presumable distribution. About 150000 steps were applied to the relaxation. Continuously, the 500 K of under temperature was specified and the cooling simulation was performed. The calculation for the relaxation at the new temperature was gone on by 100000 steps. This scheme was repeated until the specified temperature led to 500 K. All MD simulations were performed by using MXDORTO [16] that partially changed for the new three-body potential calculation.

4 Results and discussion

All MD results are obtained by using the partial ionic charges defined in equations (6) and (7). Two interatomic potential models were applied to compare the results. The first was the two-body potential model of the equation (4) (IA2), and the second was the three-body potential model of the equation (4) with equation (8) (IA3). Figure 3 shows the pair correlation functions (PCF) of B-O for x = 0.33 at 1500 K. The difference between the result of IA2 and IA3 is particular at over 0.25 nm, although the distance between the nearest neighbors is similar to each potential model result. The first local minimum of PCF means the most far point of the B-O nearest relation is 0.217 nm and 0.216 nm for IA2 and IA3 results, respectively. In this work, the threshold distance for definition of B-O bond was set 0.2 nm in consideration of the universality of this analysis toward other composition glasses. This B-O bond definition affects the calculation of the coordination numbers and analysis of ring structures in the glasses for instance.

Figure 3.

 B-O pair correlation function, P_B-O of Li₂O-B₂O₃ system at 500 K.

Figure 4 shows parts of the MD cells of x = 0.33, simulated at 1500 K by using IA2 potential (a) and by IA3 potential (b), drawing by VESTA [17]. Tetrahedrons and triangles indicate BO₄ and BO₃ units, respectively. In Figure 4 (a), there are structures similar to the undesirable units shown in Figure 1, whereas no such structure shown in Figure 4 (b). The results of the structure analysis are compiled in Table 3.

Figure 4.

 Snapshots of a part of MD cell of 0.33Li₂O-0.67B₂O₃ melt at 1500 K simulated with IA2 model (a) and with IA3 model (b).

Table 3.  Ratio, R_u of structural units (%)

		IA2 potential model					IA3 potential model
Temp./K		x					x
		0	0.1	0.2	0.33	0.4	0	0.1	0.2	0.33	0.4
2000	RBIV	2.65	6.06	31.19	12.65	12.34	0.15	2.50	6.65	7.79	5.02
	RBIII	97.20	93.89	68.8	86.95	87.38	99.85	97.45	93.24	92.09	94.85
	RBII	0.15	0.05	0.40	0.38	0.28	0.00	0.05	0.29	0.13	0.14
	ROIII	1.77	2.97	3.88	5.08	4.6	0.10	0.92	1.59	1.60	1.45
	ROII	98.13	90.77	83.27	69.20	62.19	99.90	92.65	85.00	73.84	65.08
	ROI	0.10	6.25	12.77	24.79	32.03	0.00	6.36	13.02	23.34	31.80
	RO0	0.00	0.00	0.07	0.93	1.18	0.00	0.07	0.39	1.22	1.67
1500	RBIV	4.05	11.85	31.19	19.93	20.99	4.09	2.82	6.65	11.56	11.51
	RBIII	95.95	88.14	68.8	80.07	78.94	95.91	97.18	93.24	88.37	88.50
	RBII	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.12	0.06	0.00
	ROIII	2.70	4.72	10.27	5.90	6.62	0.10	0.85	0.88	1.60	1.94
	ROII	97.30	91.05	83.27	71.43	62.42	99.9	92.99	86.94	75.51	67.82
	ROI	0.00	4.24	6.46	22.26	30.35	0.00	6.12	12.10	22.19	28.45
	RO0	0.00	0.00	0.00	0.41	0.61	0.00	0.03	0.07	0.71	1.79
1000	RBIV	6.20	19.04	32.58	41.82	39.27	0.10	3.45	9.64	16.3	19.25
	RBIII	93.80	80.97	67.43	58.18	60.73	99.9	96.55	90.36	83.69	80.76
	RBII	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	ROIII	4.13	7.38	10.48	11.84	9.43	0.07	1.09	1.80	2.71	1.83
	ROII	95.87	90.33	83.70	71.80	66.76	99.93	92.89	87.12	75.99	71.81
	ROI	0.00	2.29	5.82	16.22	23.24	0.00	6.02	10.9	20.63	25.03
	RO0	0.00	0.00	0.00	0.15	0.57	0.00	0.00	0.18	0.67	1.33
500	RBIV	8.50	22.44	36.81	49.21	51.6	0.35	4.09	11.81	21.30	23.70
	RBIII	91.50	77.57	63.19	50.78	48.40	99.65	95.91	88.19	78.70	76.29
	RBII	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	ROIII	5.67	9.26	12.1	14.25	13.39	0.23	1.40	2.33	2.67	2.74
	ROII	94.33	88.76	83.06	71.21	65.12	99.77	92.69	87.4	78.89	72.08
	ROI	0.00	1.98	4.83	14.40	20.92	0.00	5.91	10.09	17.81	24.19
	RO0	0.00	0.00	0.00	0.15	0.46	0.00	0.00	0.18	0.63	0.99

Distributions of B-O-B angle of IA2 and IA3 potential model results at 500 K are shown in Figure 5 (a) and (b). From Figure 5 (a) a drift of the angle distribution from wide to narrow with increase of x is observed. Another feature is distribution around 90 degrees that is similar to B-O-B angle of 2-membered ring structure.

Figure 5.

 B-O-B angle distribution in xLi₂O-(1−x) B₂O₃ at 500 K simulated by IA2 model (a) and by IA3 model (b).

From Figure 5 (b), there is no B-O-B angle distribution observed around 90 degrees, whereas the distribution around 172 degrees is observed for all glasses. For glasses of x ≥ 0.33, the distribution around 150 degrees increases with x.

Figure 6 shows the x dependence of R_B^IV and R_O^III. Here, the R_x^Y indicates the ratio of the Y-coordinated X element atoms for all X element atoms in the MD cell. The dashed line in Figure 6 shows the theoretical value of R_B^IV_th, calculated from equation (11).   

$R_{\text{B}}^{\text{IV}}{}_{th}=\frac{100x}{1-x},(x\le 0.5).$(11)

Figure 6.

 Ratio of structural units in xLi₂O-(1−x) B₂O₃ at 500 K and ¹¹B NMR results [12].

The value of the equation (11) is equal to the maximum of R_B^IV_th when the borate glass structure satisfies the Zachariasen's rules [14].

The R_B^IV values of IA2 potential model are exceeding the R_B^IV_th. It is the abnormal results that derived from the existence of the 2-membered ring structures as shown in Figure 1. The results of IA3 potential model indicate the glass structure is improved from the simulation by IA2 potential model. The first, the R_B^IV became smaller than the theoretical value that shown as dashed line. The second improved point is the decrease of R_O^III though O^III is not extinct.

The ring structure was investigated to confirm the 2-membered rings and to check the profile of B-O network. The ring is determined when the iteration of searching the bonded atoms that start from a boron reaches the starting boron, again. The ring size is the count of boron atoms in the searching sequence. Figure 7 shows the results of the ring analyses. N_ring is the normalized count of n-membered rings per 100 boron atoms. Figure 7 (a) shows the existence of 2-membered rings at every composition glass simulated by IA2 potential model. Two differences from Figure 7 (a) are observed in the N_ring shown in Figure 7 (b). The first is the very low value of 2-membered ring, and the second is the flat distribution of each ring size for x.

Figure 7.

 Ring size distribution in xLi₂O-(1−x) B₂O₃ at 500 K simulation applied IA2 (a) and IA3 (b) potential model.

The former result indicates the effect of the three-body potential included IA3 potential model, the 2-membered ring was decrease with decreasing of the undesirable structure unit shown in Figure 1. The latter seems to be not necessary. It is necessary to examine the details.

The decrease of B^IV and O^III results in the decrease of density of the system. The densities of glasses at several temperatures are shown in Figure 8. The densities of simulated glasses by IA2 potential model are near to the experimental data [18,19], but it is not acceptable because this density must be brought by the abnormal dense structure units. The densities of IA3 simulations are lower than the experimental, especially for small x compositions. The B₂O₃ rich melts are more viscous, therefore more relaxation steps might be required to get to the stable conditions. For small x compositions, the MD steps in this work are considered to be insufficient and large steps simulation may show a more dense structure.

Figure 8.

 Density of xLi₂O-(1−x) B₂O₃ melt/glasses. Open circle and solid circle are simulated by IA2 and IA3, respectively. Squares are the reference data at 298 K [18] and 1129 K [19].

5 Conclusion

The new three-body potential model was investigated to apply the MD simulation of Li₂O-B₂O₃ melt/glass systems. The analyzed results were compared to that of simulations by using a two-body potential model. The improvement of the melt/glass structures was confirmed, but the density of the system was not consistent with the measured values. The continuous revision of this potential model is necessary for obtaining the realistic melt/glass structures. We already started discussion and improving of the partial ionic charge model applied in this work [20].

References

• 1   T. F. Soules, J. Chem. Phys., 71, 4570 (1979).
• 2   T. F. Soules, J. Chem. Phys., 73, 4032 (1980).
• 3   T. F. Soules, J. Non-Cryst. Solids, 49, 29 (1982).
• 4   K. Hirao, N. Soga, J. Am. Ceram. Soc., 68, 515 (1985).
• 5   H. Inoue, N. Aoki, I. Yasui, J. Am. Ceram. Soc., 70, 622 (1987).
• 6   Q. Xu, K. Kawamura, T. Yokokawa, J. Non-Cryst. Solids, 104, 261 (1988).
• 7   T. F. Soules, J. Non-Cryst. Solids, 123, 48 (1990).
• 8   A. Takada, C. R. A. Catlow, G. D. Price, J. Phys. Condens. Matter, 7, 8693 (1995).
• 9   N. Ohtori, K. Takase, I. Akiyama, Y. Suzuki, K. Handa, I. Sakai, Y. Iwadate, T. Fukunaga, N. Umesaki, J. Non-Cryst. Solids, 293–295, 136 (2001).
• 10   M. Á. González, C. Mondelli, G. D'Angelo, C. Crupi, M. R. Johnson, J. Non-Cryst. Solids, 354, 203 (2008).
• 11   A. K. Varshneya, "Fundamentals of Inorganic Glasses", Academic Press. (1994).
• 12   P. J. Bray, J. G. O' Keefe, Phys. Chem. Glasses, 4, 37 (1963).
• 13   Y. F. Shepelev, R. S. Bubnova, S. K. Filatov, N. A. Sennova, N. A. Pilneva, J. Solid State Chem., 178, 2987 (2005).
• 14   W. H. Zachariasen, J. Am. Chem. Soc., 54, 3841 (1932).
• 15   N. Kumagai, K. Kawamura, T. Yokokawa, Mol. Simul., 12, 177 (1994).
• 16   K. Kawamura, MXDORTO, JCPE, P029 (1996).
• 17   VESTA, K. Momma and F. Izumi, J. Appl. Cryst., 44, 1272 (2011).
• 18   L. Shartsis, W. Capps, S. Spinner, J. Am. Ceram. Soc., 36, 35 (1953).
• 19   X. M. Shi, Q. Wang, C. X. Li, X. J. Niu, F. P. Wang, K. Q. Lu, J. Cryst. Growth, 290, 637 (2006).
• 20   Y. Yamamoto, N. Sawaguchi, M. Sasaki, J. Comput. Chem. Jpn., 13, 173 (2014).