- Email: [email protected]
Molecular dynamics simulation of heavy metal fluoride glasses: comparison of Buckingham and BHM potentials
JOURNA
ELSEVIER
L OF
Journal of Non-CrystallineSolids 184 (1995) 356-362
Molecular dynamics simulation of heavy metal fluoride glasses: comparison of Buckingham and BHM potentials Steven Gruenhut *, Douglas R. MacFarlane Department of Chemistry, Monash University, WellingtonRoad, Clayton, Victoria 3168, Australia
Abstract
The suitability for fluoride glasses of a Buckingham-potential-based molecular dynamics program, FUNGUS, has been investigated. In comparing the Buckingham potential against the traditional Born-Higgins-Mayer potential, it has been shown that the Buckingham potential can reproduce the radial distribution data of Ciccotti et al. and Lucas et al. with a discrepancy of less than 6%. The pressure of Lucas et al. has been reproduced close to the desired zero atmospheric pressure. The diffusion coefficients of Ciccotti et al. have been reproduced with a discrepancy of less than 33%. These discrepancies are consistent with slight differences in the application of the Ewald summation, different analysis techniques (Einstein versus integration of the velocity auto-correlation functions) as well as inherent statistical variations.
1. Introduction
Molecular dynamics (MD) has been growing in popularity as a result of the fact that it provides a direct route from the microscopic details of a system (e.g., atomic mass, atomic interactions, molecular geometry) to the macroscopic properties of experimental interest (e.g., equation of state, transport coefficients, structural order coefficients). Molecular dynamics has, for example, been used to analyse ion transport in solid-state ionic materials [2,3] and to study crystal growth [4] and structural relaxations near the glass transition [5,6]. MD has also been used to study the ice/water interface [7] as well as the action of catalysts [8]. A variety of substances in solid, liquid and gaseous form have been simulated by molecular dynamics. Some examples are C60 [9],
* Corresponding author. Tel: + 61-3 905 4544. Telefax: + 61-3 905 4597. E-mail: [email protected].
zeolites [8], fluoride glasses [12,29,30,2,32], semiconductors [10] and metals [4]. There are currently two main forms of potential used in simulating ionic solids, the B o r n - H u g g i n s Mayer (BHM) potential and the Buckingham potential, the vast majority of fluoride glass simulations using the former. Both these potentials are based on the Born-Mayer (BM) potential derived in 1932 [13]. This potential was later named the Buckingham potential and has the form, excluding the Coulombic part [14],
U/BUCK A/BUCKe x p [ =
riJp]
- (cBUCK/r~),
(1) where A/BucK is the depth of the potential well and includes the characteristic ionic radii, Pauling factor and p (see Eq. (2)). Additionally, p controls the steepness of the exponential repulsion and -cBUCK ij is the dipole-dipole van der Waals or dispersion term [16-18]. Although Buckingham [19] stated that the
0022-3093/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0022-3093(94)00634-2
S. Gruenhut, D.R. Macfarlane /Journal of Non-Crystalline Solids 184 (1995) 356-362 only important van der Waals term was the dipoledipole term, proportional to 1 / r 6, Mayer found that the BM potential gave polarizabilities which were too small and concluded that the dipole-quadrapole term should be also be used [20]. Later, in recalculating lattice energies and equilibrium distances between ions for alkali halide crystals, Huggins used a modified form of the BM potential which included the dipole-quadrapole van der Waals term [21]. This form of the potential function (the BHM potential) is (ignoring the Coulombic part) Ui~ m~ = A D H M b exp BHM
Cij
r6
o'i + orj-- r ij ] P
/)BHM --t) r 8.
,
(2)
where A ~ M are the Pauling factors depending on the ionic charge, o"i and o) are the characteristic radii, r/j is the inter-ionic distance and p and b are the Born repulsive parameters. The parameter b influences the depth of the potential minimum and p controls the steepness of the exponential repulsion and is also known as the hardness parameter [21,15,16]. The usage and derivation of the BM-based potentials has required a number of assumptions and approximations. The following were made in the derivation of the BM potential. (a) A BUCK was assumed to vary very slowly with r and was thus assumed to be a constant. This approximation was only considered adequate, remaining true over a small range of r, its true form being a polynomial in r [14]. (b) In arriving at the expression for the van der Waals forces it is assumed that the distance between the two atoms is sufficiently large that there is no electron interchange (no wave function overlap) [22,23]. (c) The fact that differences in structure and sizes of the ions effect the additivity of the size parameters is ignored [21]. The use of the BM-based potentials in MD relies on the interaction in ionic salts being represented by pairwise additive potentials despite the presence of many-body effects [24]. Thus, these potentials are only 'effective' potentials as they include only aver-
357
age effects of higher-body forces [25]. As a result, despite being able to calculate initial values for the potential parameters via a least-squares fit (requiring lattice parameter, expansion coefficient and compressibility and its temperature and pressure derivatives [26]), the values thus obtained should only be considered a general guide, with adaptation required to fit the desired solid as accurately as the simplified potential allows. In this contribution, a comparison is made between previous results obtained by the standard BHM potential, used with and without van der Waals terms, and the Buckingham potential via a molecular dynamics simulation using 'FUNGUS'. FUNGUS is a molecular dynamics program for ionic solids included in the suite of molecular simulation and modelling software offered by Biosym Technologies. FUNGUS is limited to a variety of forms of the Buckingham potential.
2. Computational details Two reports were chosen from the literature against which a comparison could be made. The first paper, that of Ciccotti et al., reported simulations of lithium fluoride, using a BHM potential that included the van der Waals parameters [11]. In the second paper, Lucas et al. [12] reported simulations of a zirconium-barium (ZB)-based glass (ZrTBa4F36) , ignoring the van der Waals terms in the BHM potential. In order to perform a comparison of these two potentials, the input parameters for the BHM potential, taken from the two chosen papers, were converted into their equivalent Buckingham form. In the case of the ZB glass, a BHM potential was used which had no van der Waals parameters and, being purely exponential, the parameter conversion was performed by direct calculation, i.e., A~? cK = A B H M b exp[( cri + t r j ) / p ] . (3) In the case of LiF, with the addition of van der Waals parameters in the BHM potential, obtaining a suitable set of Buckingham potential parameters is more difficult as a result of there being no dipolequadrapole van der Waals term [17]. Initially, the Buckingham potential parameters were obtained graphically for the exponential portion of the BHM
358
S. Gruenhut, D.R. Macfarlane /Journal of Non-Crystalline Solids 184 (1995) 356-362
14-
?~"
12-
10-
......+_
BHM potentials der Waal terms. Hand Tuned Buckingham potentials including Van der Waal term.
inc. Van
Li-Li
,,
>
ki-F
t-
*5
Q.
-F
0 '
I
'
'
0.5
0.0
Radius
I
1.0 / Angstroms
I
I
1.5
2.0
Fig. 1. BHM potential and hand-tuned Buckingham potential for LiF.
potential. The complete BHM potential, i.e., including the van der Waals terms, was then plotted and, using an iterative procedure, the parameters for the Buckingham potential were selected by hand such that the two functions were coincident (see Fig. 1). The short-range cut-off, used in calculating the potential, was set to half a box length (no values
were given by Ciccotti et al. [11] or Lucas et al. [12]). This is the largest spherical cut-off one can use and has the effect of increasing the simulation time as well as the accuracy of the potential calculation (a rectangular cut-off will be larger but may introduce errors in the potential calculation). The MD simulations were performed in the nor-
Table 1 FUNGUS simulation parameters
Table 2 Buckingham potential parameters for ZB and LiF
Run type
Temperature (K)
Temperature decre-
Time (ps)
ment (K) ZB simulation ~ Equilibration Quench Low T
1500 1500-300 300
LiF simulation b Equilibration Simulation
1287 1287
100 -
Box length (A)
10 13 10
14.76 14.76 14.76
2 10
13.908 13.908
a Time step 1 fs; 141 ions - 21 Zr, 12 Ba, 108 F. b Time step 1 fs;.216 ions - 108 Li, 108 F.
Pair function
A~y cK (eV)
p (~.)
ZB F-F F-Zr F-Ba Zr-Zr Zr-Ba Ba-Ba
10l 7.3 3227.98 7288.5 7952.7 19188.1 45387.1
0.266 0.266 0.266 0.266 0.266 0.266
LiF F-F F-Li Li-Li
495.0 309.0 120.2
0.303950 0.284090 0.285710
Ci~ucK (eV ,~6) 0.0 0.0 0.0 0.0 0.0 0.0
22.3 1.79 0.16
S. Gruenhut, D.R. Macfarlane / Journal of Non-Crystalline Solids 184 (1995) 356-362 Table 3 FUNGUS simulation ionic bond length data for LiF; hand calculated bond lengths for a FCC structure are in parentheses
359
Table 4 Comparison of diffusion coefficients in LiF; the simulation temperature was 1287 K and period 10000 fs
Radial distribution function
Averageionic bond length (A)
Error (%)
From Ref. [11]
FUNGUS/Buckingham potential
Li-Li Li-F F-F
3.04 (2.9) 1.80 (2.06) 3.04 (2.9)
+ 5.0 + 5.0 _+5.0
D (10 -5 cm2/s) 13.6 11.3
Error (%)
D (10 -2
Error (%)
_+3 _+3
10.9 7.7
+ 10 _+10
Li F
mal manner using FUNGUS. Simulation run parameters and potential parameters are given in Tables 1 and 2. Four complete simulations based on different starting configurations were performed for LiF and ZB so that the values obtained could be averaged and some figure of accuracy and repeatability obtained [27].
3. Results
cm2/s)
All values quoted are the average of four runs. The ionic bond lengths that appear in Tables 1 and 4 were obtained by determining the position of the first maximum of the relevant pair distribution function. Diffusion data were obtained by first fitting a straight line to the limiting slope of the mean squared diffusion curve. The gradient of this line gives the diffusion coefficient, D, using the Einstein relation [1] 6 D t = ( ] r i ( t ) - ri(O) ] 2),
The results obtained for the LiF MD simulations are listed in Tables 3 and 4 as well as Figs. 2 and 3.
(4)
where r i (t) is the position of the ith ion at time, t, and r/(0) is the original position of the ith ion.
3.5
3.0 Li-F
2.5
F-F
/ ~2.0 v ID3 I.I_
a
n- 1.5 LI-LI 1.0
0.5
/ /
//
,,'/ 0.0 ¸
3
4 Radius / Angstrom
5
Fig. 2. Radial distributionfunctionsfor LiF at 1287 K.
6
7
360
S. Gruenhut, D.R. Macfarlane/Journal of Non-Crystalline Solids 184 (1995) 356-362
Table 5 Average FUNGUSsimulationpressures (atm) for ZB
Table 6 FUNGUS simulationionic bond length data for ZB glass
Temperature (K)
Average value
Error
- 15 k - 7.7 k
+ 3750 + 7315
Radial distribution function
Averagedionic bond length (,~)
Error (%)
300 1500
F-F F-Zr F-Ba Zr-Zr Zr-Ba Ba-Ba
2.45 2.01 2.60 4.07 a 4.41 4.61
± 5.0 ± 5.0 ± 5.0 ± 5.0 ± 5.0 ± 5.0
The results obtained from the low-temperature (300 K) ZB runs are listed in Tables 5 and 6 as well as Fig. 4. All values quoted are the average of four runs. Errors reported are % relative standard deviations for the four replicates having different starting configurations.
a Small peak at 3.67 ,~
urations can produce very different diffusion rates (differences of greater than 20% were observed in the results obtained). Additionally, differences typically of up to 20% are observed between diffusion coefficients calculated via the Einstein method and via auto-correlation functions as per Ciccotti et al. Finally, it should be noted that any errors propagate exponentially with time. These propagating errors will affect the dynamic data in a much greater way than the static data [1]. Thus, if there was a substan-
4. Discussion The LiF structure shown in the radial distribution functions given by the Buckingham potential is in complete agreement with the BHM potential parameters given by Ciccotti et al. [11]. The diffusion data obtained are also very close to those obtained by Ciccotti et al. The data sets can be considered not significantly different since different starting config70
F Curve Fit: Y=-O.89(±8%)+O.OO4646(±O.4%)*X, R=0.996
60
Li C u r v e Fit: Y=-I .22(±9%)+0.007089(±0.3%)*X, R=0.998
50
Li
euCt)
E °(/)4 0 SD F
t-
,< E3 3 0 CO
20 ¸
10 ¸
0
0
2000
4000
6000
8000
10000
Time / fs
Fig. 3. Mean squared diffusionand standard deviationcurves for LiF at 1500 K. The straight line fit is also indicated.
-Zr
S. Gruenhut, D.R. Macfarlane /Journal of Non-Crystalline Solids 184 (1995) 356-362
18 16 14
Zr-Zr
/
12 U_ O
w
361
8 6
F
F
.
i~ F-Ba t /~,\ /
4
~a-Ba
•
/ /~ l ~
Zr-Ba
/
2 0
2
3
4 Radius /
5
6
7
Angstrom
Fig. 4. Radial distribution data for ZB glass at 300 K.
tial difference in the potential parameters, there could be very large differences in the diffusion data (assuming that the Coulombic component is not dominating the potential). The bond length data produced by the Buckingham potential for ZB is close to those obtained by other MD simulations and X-ray diffraction data on glasses of similar composition [28-32]. Lucas et al. [12] obtained two °adjacent Z r - Z r RDF peaks, oa major peak at 4.3 A and a smaller peak at 3.85 A. The Buckingham potential also produces two peaks, one at 4.07 A and a smaller second peak at 3.65 A. Although both peaks produced by the Buckingham potential are within experimental error of the peaks produced by Lucas et al., the smaller peak is much smaller in the FUNGUS simulations. Additionally, the smaller peak was seen in only three out of the four FUNGUS simulations. The pressure obtained in the ZB equilibration stage at 1500 K is close to zero, which is the desired pressure. Both the RDF and pressure differences observed are caused by the potential being too attractive, causing the bond lengths to be shorter and thus producing negative pressures
at 1500 K (the 300 K pressures will be smaller than the 1500 K ones and thus should always be negative). These minor differences are consistent with slight differences in the application of the Ewald summation.
5. C o n c l u s i o n
Simulations using the Buckingham potential have been compared against two, previously reported, MD simulations using the Born-Huggins-Mayer potential. It has been shown for the ZB system that the Buckingham potential can reproduce the RDF data with a discrepancy of less than 6% and to within 7000 atmospheres of the pressure obtained by Lucas et al. [12]. In the case of Ciccotti et al., the RDF data produced by the Buckingham potential are in complete agreement with the BHM potential, while the diffusion coefficients have been reproduced with a discrepancy of less than 33% [11]. These differences are consistent with slight differences in the application of the Ewald summation, different analysis tech-
362
S. Gruenhut, D.R. Macfarlane /Journal of Non-Crystalline Solids 184 (1995) 356-362
niques (Einstein versus integration of the velocity autocorrelation functions) as well as inherent statistical variations. Thus the Buckingham potential has successfully reproduced two fluoride-based simulations using the Born-Huggins-Mayer potential and shown itself to be suitable potentially for use in fluoride glass simulations.
References [1] M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids (Clarendon, Oxford, 1990). [2] C.R.A. Catlow, Solid State lonics 53-56 (1992) 955. [3] Y. Kaneko and A. Ueda, J. Phys. Soc. Jpn. 57 (1988) 3064. [4] K. Tsumuraya, M.S. Watanabe and S.K. Ikeda, Trans. ISIJ 28 (1988) 869. [5] G.F. Signorini, J. Barrat and M.L. Klein, J. Chem. Phys. 92 (1990) 1294. [6] C.A. Angell, J.H.R. Clarke and L.V. Woodcock, Adv. Chem. Phys. 48 (1981) 397. [7] O.A. Karim and A.D.J. Haymet, J. Chem. Phys. 89 (1988) 6889. [8] C.R.A. Catlow and J.M. Thomas, in: New Methods for Modelling Processes within Solids and at their Surfaces (Oxford University, Oxford, 1993) p. 61. [9] A. Cheng, M. Klein, M. Parrinello and M. Sprik, in: New Methods for Modelling Processes within Solids and at their Surfaccs (Oxford University, Oxford, 1993) p. 133. [10] U. Landman, W.D. Luedtke, M.W. Ribarsky, R.N. Barnett and C.L. Cleveland, Phys. Rev. B37 (1988) 4637.
[11] C. Ciccotti, G. Jacucci and I.R. McDonald, Phys. Rev. A13 (1976) 426. [12] J. Lucas, C.A. Angell and S. Tamaddon, Mater. Res. Bull. 19 (1984) 945. [13] M. Born and J.E. Mayer, Z. Phys. 75 (1932) 1. [14] R.A. Buckingham, Proc. R. Soc. London A168 (1938) 264. [15] D.J. Adams and I.R. McDonald, J. Phys. C7 (1974) 2761. [16] V.P.S. Nain and M.P. Sakensa, Chem. Phys. Lett. 1 (1967) 125. [17] C.R.A. Catlow, K.M. Diller and M.J. Norgett, J. Phys. C10 (1977) 1395. [18] J.O. Hirschfelder, C.F. Curtiss and R.B. Bird, Molecular Theory of Gases and Liquids (Wiley, New York, 1954). [19] R.A. Buckingham, Proc. R. Soc. London A160 (1938) 113. [20] J.E. Mayer, J. Chem. Phys. I (1993) 270. [21] M.L. Huggins, J. Chem. Phys. 5 (1937) 143. [22] J.C. Slater and J.G. Kirkwood, Phys. Rev. 37 (1931) 682. [23] H. Margenau, Rev. Mod. Phys. 11 (1939) 1. [24] L.V. Woodcock and K. Singer, Trans. Faraday Soc. 67 (1970) 12. [25] D. Fincham and D.M. Heyes, Adv. Chem. Phys. 63 (1985) 493. [26] M.J.L. Sangester and M. Dixon, Adv. Phys. 25 (1976) 247. [27] Catalysis User Guide, version 2.3.0 (Biosym Technologies, San Diego, 1993). [28] Y. Kowamoto and T. Horisaka, J. Non-Cryst. Solids 56 (1983) 39. [29] I. Yasui and H. Inoue, J. Non-Cryst. Solids 71 (1985) 39. [30] J. Lucas, D. Louer and C.A. Angell, Mater. Sci. Forum 56 (1985) 449. [31] L.T. Hamill and J.M. Parker, Phys. Chem. Glasses 26 (1985) 52. [32] C.C. Phifer, C.A. Angell, J.P. Lavel and J. Lucas, J. NonCryst. Solids 94 (1987) 315.
ELSEVIER
L OF
Journal of Non-CrystallineSolids 184 (1995) 356-362
Molecular dynamics simulation of heavy metal fluoride glasses: comparison of Buckingham and BHM potentials Steven Gruenhut *, Douglas R. MacFarlane Department of Chemistry, Monash University, WellingtonRoad, Clayton, Victoria 3168, Australia
Abstract
The suitability for fluoride glasses of a Buckingham-potential-based molecular dynamics program, FUNGUS, has been investigated. In comparing the Buckingham potential against the traditional Born-Higgins-Mayer potential, it has been shown that the Buckingham potential can reproduce the radial distribution data of Ciccotti et al. and Lucas et al. with a discrepancy of less than 6%. The pressure of Lucas et al. has been reproduced close to the desired zero atmospheric pressure. The diffusion coefficients of Ciccotti et al. have been reproduced with a discrepancy of less than 33%. These discrepancies are consistent with slight differences in the application of the Ewald summation, different analysis techniques (Einstein versus integration of the velocity auto-correlation functions) as well as inherent statistical variations.
1. Introduction
Molecular dynamics (MD) has been growing in popularity as a result of the fact that it provides a direct route from the microscopic details of a system (e.g., atomic mass, atomic interactions, molecular geometry) to the macroscopic properties of experimental interest (e.g., equation of state, transport coefficients, structural order coefficients). Molecular dynamics has, for example, been used to analyse ion transport in solid-state ionic materials [2,3] and to study crystal growth [4] and structural relaxations near the glass transition [5,6]. MD has also been used to study the ice/water interface [7] as well as the action of catalysts [8]. A variety of substances in solid, liquid and gaseous form have been simulated by molecular dynamics. Some examples are C60 [9],
* Corresponding author. Tel: + 61-3 905 4544. Telefax: + 61-3 905 4597. E-mail: [email protected].
zeolites [8], fluoride glasses [12,29,30,2,32], semiconductors [10] and metals [4]. There are currently two main forms of potential used in simulating ionic solids, the B o r n - H u g g i n s Mayer (BHM) potential and the Buckingham potential, the vast majority of fluoride glass simulations using the former. Both these potentials are based on the Born-Mayer (BM) potential derived in 1932 [13]. This potential was later named the Buckingham potential and has the form, excluding the Coulombic part [14],
U/BUCK A/BUCKe x p [ =
riJp]
- (cBUCK/r~),
(1) where A/BucK is the depth of the potential well and includes the characteristic ionic radii, Pauling factor and p (see Eq. (2)). Additionally, p controls the steepness of the exponential repulsion and -cBUCK ij is the dipole-dipole van der Waals or dispersion term [16-18]. Although Buckingham [19] stated that the
0022-3093/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0022-3093(94)00634-2
S. Gruenhut, D.R. Macfarlane /Journal of Non-Crystalline Solids 184 (1995) 356-362 only important van der Waals term was the dipoledipole term, proportional to 1 / r 6, Mayer found that the BM potential gave polarizabilities which were too small and concluded that the dipole-quadrapole term should be also be used [20]. Later, in recalculating lattice energies and equilibrium distances between ions for alkali halide crystals, Huggins used a modified form of the BM potential which included the dipole-quadrapole van der Waals term [21]. This form of the potential function (the BHM potential) is (ignoring the Coulombic part) Ui~ m~ = A D H M b exp BHM
Cij
r6
o'i + orj-- r ij ] P
/)BHM --t) r 8.
,
(2)
where A ~ M are the Pauling factors depending on the ionic charge, o"i and o) are the characteristic radii, r/j is the inter-ionic distance and p and b are the Born repulsive parameters. The parameter b influences the depth of the potential minimum and p controls the steepness of the exponential repulsion and is also known as the hardness parameter [21,15,16]. The usage and derivation of the BM-based potentials has required a number of assumptions and approximations. The following were made in the derivation of the BM potential. (a) A BUCK was assumed to vary very slowly with r and was thus assumed to be a constant. This approximation was only considered adequate, remaining true over a small range of r, its true form being a polynomial in r [14]. (b) In arriving at the expression for the van der Waals forces it is assumed that the distance between the two atoms is sufficiently large that there is no electron interchange (no wave function overlap) [22,23]. (c) The fact that differences in structure and sizes of the ions effect the additivity of the size parameters is ignored [21]. The use of the BM-based potentials in MD relies on the interaction in ionic salts being represented by pairwise additive potentials despite the presence of many-body effects [24]. Thus, these potentials are only 'effective' potentials as they include only aver-
357
age effects of higher-body forces [25]. As a result, despite being able to calculate initial values for the potential parameters via a least-squares fit (requiring lattice parameter, expansion coefficient and compressibility and its temperature and pressure derivatives [26]), the values thus obtained should only be considered a general guide, with adaptation required to fit the desired solid as accurately as the simplified potential allows. In this contribution, a comparison is made between previous results obtained by the standard BHM potential, used with and without van der Waals terms, and the Buckingham potential via a molecular dynamics simulation using 'FUNGUS'. FUNGUS is a molecular dynamics program for ionic solids included in the suite of molecular simulation and modelling software offered by Biosym Technologies. FUNGUS is limited to a variety of forms of the Buckingham potential.
2. Computational details Two reports were chosen from the literature against which a comparison could be made. The first paper, that of Ciccotti et al., reported simulations of lithium fluoride, using a BHM potential that included the van der Waals parameters [11]. In the second paper, Lucas et al. [12] reported simulations of a zirconium-barium (ZB)-based glass (ZrTBa4F36) , ignoring the van der Waals terms in the BHM potential. In order to perform a comparison of these two potentials, the input parameters for the BHM potential, taken from the two chosen papers, were converted into their equivalent Buckingham form. In the case of the ZB glass, a BHM potential was used which had no van der Waals parameters and, being purely exponential, the parameter conversion was performed by direct calculation, i.e., A~? cK = A B H M b exp[( cri + t r j ) / p ] . (3) In the case of LiF, with the addition of van der Waals parameters in the BHM potential, obtaining a suitable set of Buckingham potential parameters is more difficult as a result of there being no dipolequadrapole van der Waals term [17]. Initially, the Buckingham potential parameters were obtained graphically for the exponential portion of the BHM
358
S. Gruenhut, D.R. Macfarlane /Journal of Non-Crystalline Solids 184 (1995) 356-362
14-
?~"
12-
10-
......+_
BHM potentials der Waal terms. Hand Tuned Buckingham potentials including Van der Waal term.
inc. Van
Li-Li
,,
>
ki-F
t-
*5
Q.
-F
0 '
I
'
'
0.5
0.0
Radius
I
1.0 / Angstroms
I
I
1.5
2.0
Fig. 1. BHM potential and hand-tuned Buckingham potential for LiF.
potential. The complete BHM potential, i.e., including the van der Waals terms, was then plotted and, using an iterative procedure, the parameters for the Buckingham potential were selected by hand such that the two functions were coincident (see Fig. 1). The short-range cut-off, used in calculating the potential, was set to half a box length (no values
were given by Ciccotti et al. [11] or Lucas et al. [12]). This is the largest spherical cut-off one can use and has the effect of increasing the simulation time as well as the accuracy of the potential calculation (a rectangular cut-off will be larger but may introduce errors in the potential calculation). The MD simulations were performed in the nor-
Table 1 FUNGUS simulation parameters
Table 2 Buckingham potential parameters for ZB and LiF
Run type
Temperature (K)
Temperature decre-
Time (ps)
ment (K) ZB simulation ~ Equilibration Quench Low T
1500 1500-300 300
LiF simulation b Equilibration Simulation
1287 1287
100 -
Box length (A)
10 13 10
14.76 14.76 14.76
2 10
13.908 13.908
a Time step 1 fs; 141 ions - 21 Zr, 12 Ba, 108 F. b Time step 1 fs;.216 ions - 108 Li, 108 F.
Pair function
A~y cK (eV)
p (~.)
ZB F-F F-Zr F-Ba Zr-Zr Zr-Ba Ba-Ba
10l 7.3 3227.98 7288.5 7952.7 19188.1 45387.1
0.266 0.266 0.266 0.266 0.266 0.266
LiF F-F F-Li Li-Li
495.0 309.0 120.2
0.303950 0.284090 0.285710
Ci~ucK (eV ,~6) 0.0 0.0 0.0 0.0 0.0 0.0
22.3 1.79 0.16
S. Gruenhut, D.R. Macfarlane / Journal of Non-Crystalline Solids 184 (1995) 356-362 Table 3 FUNGUS simulation ionic bond length data for LiF; hand calculated bond lengths for a FCC structure are in parentheses
359
Table 4 Comparison of diffusion coefficients in LiF; the simulation temperature was 1287 K and period 10000 fs
Radial distribution function
Averageionic bond length (A)
Error (%)
From Ref. [11]
FUNGUS/Buckingham potential
Li-Li Li-F F-F
3.04 (2.9) 1.80 (2.06) 3.04 (2.9)
+ 5.0 + 5.0 _+5.0
D (10 -5 cm2/s) 13.6 11.3
Error (%)
D (10 -2
Error (%)
_+3 _+3
10.9 7.7
+ 10 _+10
Li F
mal manner using FUNGUS. Simulation run parameters and potential parameters are given in Tables 1 and 2. Four complete simulations based on different starting configurations were performed for LiF and ZB so that the values obtained could be averaged and some figure of accuracy and repeatability obtained [27].
3. Results
cm2/s)
All values quoted are the average of four runs. The ionic bond lengths that appear in Tables 1 and 4 were obtained by determining the position of the first maximum of the relevant pair distribution function. Diffusion data were obtained by first fitting a straight line to the limiting slope of the mean squared diffusion curve. The gradient of this line gives the diffusion coefficient, D, using the Einstein relation [1] 6 D t = ( ] r i ( t ) - ri(O) ] 2),
The results obtained for the LiF MD simulations are listed in Tables 3 and 4 as well as Figs. 2 and 3.
(4)
where r i (t) is the position of the ith ion at time, t, and r/(0) is the original position of the ith ion.
3.5
3.0 Li-F
2.5
F-F
/ ~2.0 v ID3 I.I_
a
n- 1.5 LI-LI 1.0
0.5
/ /
//
,,'/ 0.0 ¸
3
4 Radius / Angstrom
5
Fig. 2. Radial distributionfunctionsfor LiF at 1287 K.
6
7
360
S. Gruenhut, D.R. Macfarlane/Journal of Non-Crystalline Solids 184 (1995) 356-362
Table 5 Average FUNGUSsimulationpressures (atm) for ZB
Table 6 FUNGUS simulationionic bond length data for ZB glass
Temperature (K)
Average value
Error
- 15 k - 7.7 k
+ 3750 + 7315
Radial distribution function
Averagedionic bond length (,~)
Error (%)
300 1500
F-F F-Zr F-Ba Zr-Zr Zr-Ba Ba-Ba
2.45 2.01 2.60 4.07 a 4.41 4.61
± 5.0 ± 5.0 ± 5.0 ± 5.0 ± 5.0 ± 5.0
The results obtained from the low-temperature (300 K) ZB runs are listed in Tables 5 and 6 as well as Fig. 4. All values quoted are the average of four runs. Errors reported are % relative standard deviations for the four replicates having different starting configurations.
a Small peak at 3.67 ,~
urations can produce very different diffusion rates (differences of greater than 20% were observed in the results obtained). Additionally, differences typically of up to 20% are observed between diffusion coefficients calculated via the Einstein method and via auto-correlation functions as per Ciccotti et al. Finally, it should be noted that any errors propagate exponentially with time. These propagating errors will affect the dynamic data in a much greater way than the static data [1]. Thus, if there was a substan-
4. Discussion The LiF structure shown in the radial distribution functions given by the Buckingham potential is in complete agreement with the BHM potential parameters given by Ciccotti et al. [11]. The diffusion data obtained are also very close to those obtained by Ciccotti et al. The data sets can be considered not significantly different since different starting config70
F Curve Fit: Y=-O.89(±8%)+O.OO4646(±O.4%)*X, R=0.996
60
Li C u r v e Fit: Y=-I .22(±9%)+0.007089(±0.3%)*X, R=0.998
50
Li
euCt)
E °(/)4 0 SD F
t-
,< E3 3 0 CO
20 ¸
10 ¸
0
0
2000
4000
6000
8000
10000
Time / fs
Fig. 3. Mean squared diffusionand standard deviationcurves for LiF at 1500 K. The straight line fit is also indicated.
-Zr
S. Gruenhut, D.R. Macfarlane /Journal of Non-Crystalline Solids 184 (1995) 356-362
18 16 14
Zr-Zr
/
12 U_ O
w
361
8 6
F
F
.
i~ F-Ba t /~,\ /
4
~a-Ba
•
/ /~ l ~
Zr-Ba
/
2 0
2
3
4 Radius /
5
6
7
Angstrom
Fig. 4. Radial distribution data for ZB glass at 300 K.
tial difference in the potential parameters, there could be very large differences in the diffusion data (assuming that the Coulombic component is not dominating the potential). The bond length data produced by the Buckingham potential for ZB is close to those obtained by other MD simulations and X-ray diffraction data on glasses of similar composition [28-32]. Lucas et al. [12] obtained two °adjacent Z r - Z r RDF peaks, oa major peak at 4.3 A and a smaller peak at 3.85 A. The Buckingham potential also produces two peaks, one at 4.07 A and a smaller second peak at 3.65 A. Although both peaks produced by the Buckingham potential are within experimental error of the peaks produced by Lucas et al., the smaller peak is much smaller in the FUNGUS simulations. Additionally, the smaller peak was seen in only three out of the four FUNGUS simulations. The pressure obtained in the ZB equilibration stage at 1500 K is close to zero, which is the desired pressure. Both the RDF and pressure differences observed are caused by the potential being too attractive, causing the bond lengths to be shorter and thus producing negative pressures
at 1500 K (the 300 K pressures will be smaller than the 1500 K ones and thus should always be negative). These minor differences are consistent with slight differences in the application of the Ewald summation.
5. C o n c l u s i o n
Simulations using the Buckingham potential have been compared against two, previously reported, MD simulations using the Born-Huggins-Mayer potential. It has been shown for the ZB system that the Buckingham potential can reproduce the RDF data with a discrepancy of less than 6% and to within 7000 atmospheres of the pressure obtained by Lucas et al. [12]. In the case of Ciccotti et al., the RDF data produced by the Buckingham potential are in complete agreement with the BHM potential, while the diffusion coefficients have been reproduced with a discrepancy of less than 33% [11]. These differences are consistent with slight differences in the application of the Ewald summation, different analysis tech-
362
S. Gruenhut, D.R. Macfarlane /Journal of Non-Crystalline Solids 184 (1995) 356-362
niques (Einstein versus integration of the velocity autocorrelation functions) as well as inherent statistical variations. Thus the Buckingham potential has successfully reproduced two fluoride-based simulations using the Born-Huggins-Mayer potential and shown itself to be suitable potentially for use in fluoride glass simulations.
References [1] M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids (Clarendon, Oxford, 1990). [2] C.R.A. Catlow, Solid State lonics 53-56 (1992) 955. [3] Y. Kaneko and A. Ueda, J. Phys. Soc. Jpn. 57 (1988) 3064. [4] K. Tsumuraya, M.S. Watanabe and S.K. Ikeda, Trans. ISIJ 28 (1988) 869. [5] G.F. Signorini, J. Barrat and M.L. Klein, J. Chem. Phys. 92 (1990) 1294. [6] C.A. Angell, J.H.R. Clarke and L.V. Woodcock, Adv. Chem. Phys. 48 (1981) 397. [7] O.A. Karim and A.D.J. Haymet, J. Chem. Phys. 89 (1988) 6889. [8] C.R.A. Catlow and J.M. Thomas, in: New Methods for Modelling Processes within Solids and at their Surfaces (Oxford University, Oxford, 1993) p. 61. [9] A. Cheng, M. Klein, M. Parrinello and M. Sprik, in: New Methods for Modelling Processes within Solids and at their Surfaccs (Oxford University, Oxford, 1993) p. 133. [10] U. Landman, W.D. Luedtke, M.W. Ribarsky, R.N. Barnett and C.L. Cleveland, Phys. Rev. B37 (1988) 4637.
[11] C. Ciccotti, G. Jacucci and I.R. McDonald, Phys. Rev. A13 (1976) 426. [12] J. Lucas, C.A. Angell and S. Tamaddon, Mater. Res. Bull. 19 (1984) 945. [13] M. Born and J.E. Mayer, Z. Phys. 75 (1932) 1. [14] R.A. Buckingham, Proc. R. Soc. London A168 (1938) 264. [15] D.J. Adams and I.R. McDonald, J. Phys. C7 (1974) 2761. [16] V.P.S. Nain and M.P. Sakensa, Chem. Phys. Lett. 1 (1967) 125. [17] C.R.A. Catlow, K.M. Diller and M.J. Norgett, J. Phys. C10 (1977) 1395. [18] J.O. Hirschfelder, C.F. Curtiss and R.B. Bird, Molecular Theory of Gases and Liquids (Wiley, New York, 1954). [19] R.A. Buckingham, Proc. R. Soc. London A160 (1938) 113. [20] J.E. Mayer, J. Chem. Phys. I (1993) 270. [21] M.L. Huggins, J. Chem. Phys. 5 (1937) 143. [22] J.C. Slater and J.G. Kirkwood, Phys. Rev. 37 (1931) 682. [23] H. Margenau, Rev. Mod. Phys. 11 (1939) 1. [24] L.V. Woodcock and K. Singer, Trans. Faraday Soc. 67 (1970) 12. [25] D. Fincham and D.M. Heyes, Adv. Chem. Phys. 63 (1985) 493. [26] M.J.L. Sangester and M. Dixon, Adv. Phys. 25 (1976) 247. [27] Catalysis User Guide, version 2.3.0 (Biosym Technologies, San Diego, 1993). [28] Y. Kowamoto and T. Horisaka, J. Non-Cryst. Solids 56 (1983) 39. [29] I. Yasui and H. Inoue, J. Non-Cryst. Solids 71 (1985) 39. [30] J. Lucas, D. Louer and C.A. Angell, Mater. Sci. Forum 56 (1985) 449. [31] L.T. Hamill and J.M. Parker, Phys. Chem. Glasses 26 (1985) 52. [32] C.C. Phifer, C.A. Angell, J.P. Lavel and J. Lucas, J. NonCryst. Solids 94 (1987) 315.