- Email: [email protected]
GAMMEL: A program for the calculation of the pressure dependence of the elastic constants of ionic crystals
Pergamon
Compurers Chem. Vol. 20, No. 4, pp. 485-487, 1996 Copyright 0 1996 Elscvicr Science Ltd Printed in Great Britain. All rights reserved 0097-8485/96 $15.00 + 0.00
0097-8485(95)ooo9o-9
SOFTWARE NOTE GAMMEL: A PROGRAM FOR THE CALCULATION OF THE PRESSURE DEPENDENCE OF THE ELASTIC CONSTANTS OF IONIC CRYSTALS* M. R. SORIANO, J. A. 0. BRUNO and A. BATANAt Grupo de Quimica Tebrica, Departamento de Quimica Inorginica, Analitica y Quimica Fisica, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires, Buenos Aires, Argentina (Received I1 September 1995; accepted in final form 27 November 1995)
Ahatrati-Program GAMMEL has been developed for the calculation of different elastic properties of ionic solids, in particular their pressure dependence. Sodium chloride, caesium chloride and the fluorite structures were studied; an analysis of central and non-central forces was allowed for, considering up to second-neighbour interactions using different potential forms and parametrizations. The program has a flexible structure and allows for a check with experimental elastic data and a cross-checking with an anharmonic property (thermal expansion in the limit of low temperatures). Copyright 0 1996 Elsevier
Science Ltd
the pressure dependence of the second-order elastic constants with different existing models (see below) with interionic pair potentials of seven different forms
1. INTRODUCI’ION In the study of the anharmonic properties of ionic crystals it is important to have an adequate interaction model. This can be looked for through a good fit of the elastic data and its pressure dependence in the limit of low temperatures where the solid can be considered as an elastic continuum. Unfortunately, experimental data at high pressures are subject to considerable uncertainty (Barron, in press), which makes this comparison less reliable, and besides there are very few elastic data at low temperatures. On the other hand, experimental data of low temperature thermal expansion have been obtained with increasing precision (Barron et al., 1980; Barron, in press) so that it can be used as an indirect check of the model. At low temperatures the quasiharmonic approximation can be applied; in this temperature range the results are sensitive to the interaction models chosen (Barron, in press). Thermal expansion is usually studied through the Griineisen function y lb(T), which in the limit of T -+0 K should tend to the yg value calculated from elastic data (either theoretical or experimental). Program GAMMEL was developed in FORTRAN77 language allowing the evaluation of
(Dutt et al., 1985) and giving the chance of trying different parametrizations; central and non-central forces are considered, including up to secondneighbour interactions. Comparison with experimental data is allowed for through the Griineisen parameter in the limit of T + 0 K (ri’ and yt). Different intermediate quantities (SUMCIJ, see below) are defined as a further check of the model and parametrizations used. Ionic crystals of sodium chloride, caesium chloride and fluorite structures can be studied. The program runs on personal computers and was successfully compiled without special options with several compilers. The code is easily portable to other systems. 2. METHOD OF CALCULATION Thermal expansion is usually studied through the Griineisen function
(1)
*A copy of the program is available on request from: Dra. Alicia Batana, Departamento de Quimica InorgzWca, Analitica y Qulmica Fisica, Fact&ad de Ciencias Exactas y Naturales, Pahellbn II, Ciudad Universitaria, 1428 Buenos Aires, Argentina. t Author to whom all correspondence should be addressed.
where V is the molar volume, x1 is the adiabatic compressibility and C,, is the molar heat capacity at constant pressure. The tineisen function y(T)in the low temperature limit (hereafter noted as yO)can be calculated independently using experimental data of the elastic constants cYand their pressure derivatives and these results can be compared with those obtained from cg and &,/ap calculated with theoretical models.
485
486
Software Note parametrization used because its value can be calculated from experimental data of cii and (ac@p):
For cubic crystals:
” = -
a In 0, a In v
SUMCIJ = -2
ac, 6 2XTC44 ap
=- I+-_ I
=
YI +
Y2 +
(2)
Y3
where 0, is the Debye temperature, f(.s, t) is the anisotropy factor and is a function of the independent variables s and t, where s = (c,i - c~)/(+ + c~) and t = (c,~ - c,.,)/Q; f(s, t) is calculated using the de Launay’s tables (de Launay, 1956, 1959) by interpolation (Akima, 1974). The interaction models used are those compatible with previous shell model calculations (Woods et al., 1962; Cowley, 1962; Cowley et al., 1963), particularly for those considering up to second-neighbour interactions, including radial and tangential forces between ions of type l-2, l-l and 2-2 (1 and 2 designate cations and anions, respectively). The use of an effective ionic charge is included, instead of a full formal one. The calculation of (i$/ap) is performed using the following equations: ac,, - -11f(SUMCll ap3
-c,,+Kll)
(3)
where: SUMCll
=~;‘;,+~(~T;,+0;‘;,+9~~,+~~,) 2
Kll = 3.834 $ 0
ac12 - $(SUMClZ -=
- 2c,, + K12)
ap
(4)
where: SUMC12 = -2d;;,+
a[(+;;,
+ &,)
- 5(4;;, + &?,)I K12 = -0.3389
ah -=
ap
0
2
F
-% (SUMC44 - 2c, + K44)
(5)
where: SUMC44 = 24 ;;, + ,&J
;:;, + &,) + 3($;:;, +
K44= -3.8341
$
3 XT
G,)l
’
0 based on the formalism of cii and (&&?p) previously given (Soriano & Batana, in press). The intermediate quantity SUMCIJ can also be used as a check of the interaction model and the
+ 2~ + KIJ.
(6)
aP
Seven different forms can be used for the interionic potential 4(r) as was mentioned before, namely (i) the Born-Land& potential (Ar -“), (ii) the BomMayer potential (B exp( -r/p)), (iii) the Hellmann potential ((B/r)exp( - r/p,)), (iv) the Wasastjema potential (Cr’ exp( -fir)), (v) the Varshni-Shukla potential (1 exp( -kr2)), (vi) the Modified Varshni-Shukla potential (1, exp( -kr3/2)) and (vii) a potential built as follows: 4r2 = Arm”, &, and +22= A(d2r)-6. Here, d’, 4” and 4”’ are the derivatives of 4(r) with respect to r.
3. DESCRIP’IION OF THE PROGRAM Program GAMMEL consists of a main module linked with 18 subroutines. Data and results are read and written to external files. A brief description of the typical input data is as follows: l LX: number
of input grid points in the xcoordinate (must be two or greater); s in de Launay’s tables. l LY: number of input grid points in the ycoordinate (must be two or greater); t in de Launay’s tables. l LTX: an integer number that identifies the lattice type crystal, namely 1 = fluorite structure, 2 = CsCl structure and 3 = NaCl structure. l TEMP: absolute temperature. l LATYPE: lattice structure. 0 Cl 1, C12, C44: experimental elastic constants (cgs units). l DC1 1, DC12, DC44: experimental pressure derivatives of the elastic constants. l RO: cation-cation distance. lZ1, 22: ionic charges. l DTXINP: inverse volumetric isothermal compressibility derived with respect to pressure. lY2, ZK2, DF122, D2FI22: parameters for the model. lAl2, Al 1, Bll, A22, B22, B: radial and tangential force constants between 1-2, l-l or 2-2 next-nearest neighbours (1 = cation, 2 = anion). If the corresponding model does not need these parameters, zeros should be entered anyway. During the preparation of the input data file, several switches need to be included, allowing (for example) to choose the number of calculations to be done for each crystal, identification of the model to be employed, inclusion or not of non-central interactions, different combinations of experimental data according to their availability, type of interaction potential, and so on.
Software Note
487
Table 1. Some example elastic properties results obtained with GAMMEL. The elastic constants are given in units of 10-‘2dyncn-‘. ya in experimental columns are the values calculated with experimental data LiF
NaCl
Nal
Exp.*
Calc’d
Exp.t
Calc’d
Exp.f
Calc’d
CII CIZ &dp
1.1370 0.4760 9.92 0.6370
1.1276 0.4760 9.781 0.6469
0.5733 0.1228 11.47 0.1323
0.5991 0.1326 11.79 0.1380
0.3590 0.0750 11.83 0.0768
0.3524 0.0794 10.866 0.0766
wap Yo
2.72 1.7192
2.820 1.901
1.91 1.1675
2.01 I 1.3155
2.30 1.3327
2.67 1.3573
* From Namjoshi ef al. (1971). t From Lewis er al. (1967). $ From Woods ef al. (1962).
4. EXAMPLE
In Table 1 the results for a calculation of the elastic constants and their pressure dependence are given and compared with experimental data. In this case, dC,/dP was introduced as input data (this is optional). y,, was calculated both with experimental elastic data (columns 2,4 and 6) and with calculated data (columns 3, 5 and 7). The models used were: LiF and NaI a non-central-force model for negative next-nearest neighbours with effective charge and potential VII (see Section 2); for NaCl similar to the previous one but including also non-central forces for next-nearest positive-ion interactions. Considering the experimental uncertainty of the magnitudes under study, the results obtained with GAMMEL show a good agreement with the experimental values. As expected from the cation-anion relative-size analysis the inclusion of next-nearest positive-ion interactions was necessary for NaCl in order to improve the quantitative agreement with experiment. 5. CONCLUSIONS GAMMEL performs the calculation of elastic constants and their pressure derivatives for sodium chloride, caesium chloride and fluorite-structured crystals. An aspect that we found particularly useful is that it allows, with the use of relatively simple models, a search for: (i) most adequate interaction potential model; (ii) most adequate interionic potential form; and (iii) most adequate parametrization, with partial results to find a failure of a given parameters set (as with SUMCIJ). The expressions used within the program allow for a future use of independent parametrizations (e.g. ub initio or semiempirical). Good results were obtained using simple models with few parameters (see Table 1). Another interesting capability of the program is that it includes the evaluation of y0 (Griineisen function at T -+ 0 K) allowing: (i) the cross-checking of the model used, owing to the more reliable experimental data of yt;
(ii) the suggestion of a starting point with respect to the interaction models to be used during a lattice dynamical calculation and in the study of anharmonic properties, taking the sensitivity, e.g. of y(T), to the potential used within the low temperature region; and (iii) the separation of different contributions to y,,, particularly y2 and y3 [see equation (3)], whose relative values appear to be characteristic of the structure considered and may therefore allow for a check in the prediction of y,, values for compounds with no experimental data available. It is worth taking into account that this method for the evaluation of y,, avoids the integration over the reciprocal space that is very often a source of errors in these calculations (Barron & Batana, 1968; Barron et al., 1980). Acknowledgements-The authors gratefully thank Dr T. H. K. Barron for his very helpful comments and dis-
cussions. The authors are grateful to the referees for their valuable comments, which have been found very useful for revising the manuscript. We also acknowledge financial support from Universidad de Buenos Aires and CONICET.
REFERENCES Akima H. (1974) Comm. ACM 17(l), 26. Barron T. H. K. (in press) Generalised Theory of Thermal Expansion of Solids, Chapter 1 of Vol. I-4. CINDAS Series on Thermophysical Properties. Barron T. H. K. & Batana A. (1968) Whys. Reu. 167, 814. Barron T. H. K., Collins J. G. & White G. K. (1980) Adv. Phys. 29, 609. Cowley R. A. (1962) Proc. Roy. Sot. A268, 121. Cowley R. A., .&&ran W., Brockhouse B. N. 8r Woods A. D. B. (1963) Phvs. Rev. 131. 1030. Dutt N., Agrawal’G. G. & Shanker J. (1985) Phys. Stat. Sol. (b) 132, 99. de Launay J. (1956) J. Chem. Phys. 24, 1071. de Launay J. (1959) J. Chem. Phys. 30, 91. Lewis J. T., Lehozky A. & Briswe C. V. (1967) Phys. Rev. 161, 877. Namjoshi K. V., Mitra S. S. KLVetelino J. F. (1971) Phys. Rev. B 12, 4398. Soriano M. R. & Batana A. (in press) An. Asoc. Quim. Arg. Woods A. D. B., Cochran W. & Brockhouse B. N. (1962) Phys. Rev. 119, 980.
Compurers Chem. Vol. 20, No. 4, pp. 485-487, 1996 Copyright 0 1996 Elscvicr Science Ltd Printed in Great Britain. All rights reserved 0097-8485/96 $15.00 + 0.00
0097-8485(95)ooo9o-9
SOFTWARE NOTE GAMMEL: A PROGRAM FOR THE CALCULATION OF THE PRESSURE DEPENDENCE OF THE ELASTIC CONSTANTS OF IONIC CRYSTALS* M. R. SORIANO, J. A. 0. BRUNO and A. BATANAt Grupo de Quimica Tebrica, Departamento de Quimica Inorginica, Analitica y Quimica Fisica, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires, Buenos Aires, Argentina (Received I1 September 1995; accepted in final form 27 November 1995)
Ahatrati-Program GAMMEL has been developed for the calculation of different elastic properties of ionic solids, in particular their pressure dependence. Sodium chloride, caesium chloride and the fluorite structures were studied; an analysis of central and non-central forces was allowed for, considering up to second-neighbour interactions using different potential forms and parametrizations. The program has a flexible structure and allows for a check with experimental elastic data and a cross-checking with an anharmonic property (thermal expansion in the limit of low temperatures). Copyright 0 1996 Elsevier
Science Ltd
the pressure dependence of the second-order elastic constants with different existing models (see below) with interionic pair potentials of seven different forms
1. INTRODUCI’ION In the study of the anharmonic properties of ionic crystals it is important to have an adequate interaction model. This can be looked for through a good fit of the elastic data and its pressure dependence in the limit of low temperatures where the solid can be considered as an elastic continuum. Unfortunately, experimental data at high pressures are subject to considerable uncertainty (Barron, in press), which makes this comparison less reliable, and besides there are very few elastic data at low temperatures. On the other hand, experimental data of low temperature thermal expansion have been obtained with increasing precision (Barron et al., 1980; Barron, in press) so that it can be used as an indirect check of the model. At low temperatures the quasiharmonic approximation can be applied; in this temperature range the results are sensitive to the interaction models chosen (Barron, in press). Thermal expansion is usually studied through the Griineisen function y lb(T), which in the limit of T -+0 K should tend to the yg value calculated from elastic data (either theoretical or experimental). Program GAMMEL was developed in FORTRAN77 language allowing the evaluation of
(Dutt et al., 1985) and giving the chance of trying different parametrizations; central and non-central forces are considered, including up to secondneighbour interactions. Comparison with experimental data is allowed for through the Griineisen parameter in the limit of T + 0 K (ri’ and yt). Different intermediate quantities (SUMCIJ, see below) are defined as a further check of the model and parametrizations used. Ionic crystals of sodium chloride, caesium chloride and fluorite structures can be studied. The program runs on personal computers and was successfully compiled without special options with several compilers. The code is easily portable to other systems. 2. METHOD OF CALCULATION Thermal expansion is usually studied through the Griineisen function
(1)
*A copy of the program is available on request from: Dra. Alicia Batana, Departamento de Quimica InorgzWca, Analitica y Qulmica Fisica, Fact&ad de Ciencias Exactas y Naturales, Pahellbn II, Ciudad Universitaria, 1428 Buenos Aires, Argentina. t Author to whom all correspondence should be addressed.
where V is the molar volume, x1 is the adiabatic compressibility and C,, is the molar heat capacity at constant pressure. The tineisen function y(T)in the low temperature limit (hereafter noted as yO)can be calculated independently using experimental data of the elastic constants cYand their pressure derivatives and these results can be compared with those obtained from cg and &,/ap calculated with theoretical models.
485
486
Software Note parametrization used because its value can be calculated from experimental data of cii and (ac@p):
For cubic crystals:
” = -
a In 0, a In v
SUMCIJ = -2
ac, 6 2XTC44 ap
=- I+-_ I
=
YI +
Y2 +
(2)
Y3
where 0, is the Debye temperature, f(.s, t) is the anisotropy factor and is a function of the independent variables s and t, where s = (c,i - c~)/(+ + c~) and t = (c,~ - c,.,)/Q; f(s, t) is calculated using the de Launay’s tables (de Launay, 1956, 1959) by interpolation (Akima, 1974). The interaction models used are those compatible with previous shell model calculations (Woods et al., 1962; Cowley, 1962; Cowley et al., 1963), particularly for those considering up to second-neighbour interactions, including radial and tangential forces between ions of type l-2, l-l and 2-2 (1 and 2 designate cations and anions, respectively). The use of an effective ionic charge is included, instead of a full formal one. The calculation of (i$/ap) is performed using the following equations: ac,, - -11f(SUMCll ap3
-c,,+Kll)
(3)
where: SUMCll
=~;‘;,+~(~T;,+0;‘;,+9~~,+~~,) 2
Kll = 3.834 $ 0
ac12 - $(SUMClZ -=
- 2c,, + K12)
ap
(4)
where: SUMC12 = -2d;;,+
a[(+;;,
+ &,)
- 5(4;;, + &?,)I K12 = -0.3389
ah -=
ap
0
2
F
-% (SUMC44 - 2c, + K44)
(5)
where: SUMC44 = 24 ;;, + ,&J
;:;, + &,) + 3($;:;, +
K44= -3.8341
$
3 XT
G,)l
’
0 based on the formalism of cii and (&&?p) previously given (Soriano & Batana, in press). The intermediate quantity SUMCIJ can also be used as a check of the interaction model and the
+ 2~ + KIJ.
(6)
aP
Seven different forms can be used for the interionic potential 4(r) as was mentioned before, namely (i) the Born-Land& potential (Ar -“), (ii) the BomMayer potential (B exp( -r/p)), (iii) the Hellmann potential ((B/r)exp( - r/p,)), (iv) the Wasastjema potential (Cr’ exp( -fir)), (v) the Varshni-Shukla potential (1 exp( -kr2)), (vi) the Modified Varshni-Shukla potential (1, exp( -kr3/2)) and (vii) a potential built as follows: 4r2 = Arm”, &, and +22= A(d2r)-6. Here, d’, 4” and 4”’ are the derivatives of 4(r) with respect to r.
3. DESCRIP’IION OF THE PROGRAM Program GAMMEL consists of a main module linked with 18 subroutines. Data and results are read and written to external files. A brief description of the typical input data is as follows: l LX: number
of input grid points in the xcoordinate (must be two or greater); s in de Launay’s tables. l LY: number of input grid points in the ycoordinate (must be two or greater); t in de Launay’s tables. l LTX: an integer number that identifies the lattice type crystal, namely 1 = fluorite structure, 2 = CsCl structure and 3 = NaCl structure. l TEMP: absolute temperature. l LATYPE: lattice structure. 0 Cl 1, C12, C44: experimental elastic constants (cgs units). l DC1 1, DC12, DC44: experimental pressure derivatives of the elastic constants. l RO: cation-cation distance. lZ1, 22: ionic charges. l DTXINP: inverse volumetric isothermal compressibility derived with respect to pressure. lY2, ZK2, DF122, D2FI22: parameters for the model. lAl2, Al 1, Bll, A22, B22, B: radial and tangential force constants between 1-2, l-l or 2-2 next-nearest neighbours (1 = cation, 2 = anion). If the corresponding model does not need these parameters, zeros should be entered anyway. During the preparation of the input data file, several switches need to be included, allowing (for example) to choose the number of calculations to be done for each crystal, identification of the model to be employed, inclusion or not of non-central interactions, different combinations of experimental data according to their availability, type of interaction potential, and so on.
Software Note
487
Table 1. Some example elastic properties results obtained with GAMMEL. The elastic constants are given in units of 10-‘2dyncn-‘. ya in experimental columns are the values calculated with experimental data LiF
NaCl
Nal
Exp.*
Calc’d
Exp.t
Calc’d
Exp.f
Calc’d
CII CIZ &dp
1.1370 0.4760 9.92 0.6370
1.1276 0.4760 9.781 0.6469
0.5733 0.1228 11.47 0.1323
0.5991 0.1326 11.79 0.1380
0.3590 0.0750 11.83 0.0768
0.3524 0.0794 10.866 0.0766
wap Yo
2.72 1.7192
2.820 1.901
1.91 1.1675
2.01 I 1.3155
2.30 1.3327
2.67 1.3573
* From Namjoshi ef al. (1971). t From Lewis er al. (1967). $ From Woods ef al. (1962).
4. EXAMPLE
In Table 1 the results for a calculation of the elastic constants and their pressure dependence are given and compared with experimental data. In this case, dC,/dP was introduced as input data (this is optional). y,, was calculated both with experimental elastic data (columns 2,4 and 6) and with calculated data (columns 3, 5 and 7). The models used were: LiF and NaI a non-central-force model for negative next-nearest neighbours with effective charge and potential VII (see Section 2); for NaCl similar to the previous one but including also non-central forces for next-nearest positive-ion interactions. Considering the experimental uncertainty of the magnitudes under study, the results obtained with GAMMEL show a good agreement with the experimental values. As expected from the cation-anion relative-size analysis the inclusion of next-nearest positive-ion interactions was necessary for NaCl in order to improve the quantitative agreement with experiment. 5. CONCLUSIONS GAMMEL performs the calculation of elastic constants and their pressure derivatives for sodium chloride, caesium chloride and fluorite-structured crystals. An aspect that we found particularly useful is that it allows, with the use of relatively simple models, a search for: (i) most adequate interaction potential model; (ii) most adequate interionic potential form; and (iii) most adequate parametrization, with partial results to find a failure of a given parameters set (as with SUMCIJ). The expressions used within the program allow for a future use of independent parametrizations (e.g. ub initio or semiempirical). Good results were obtained using simple models with few parameters (see Table 1). Another interesting capability of the program is that it includes the evaluation of y0 (Griineisen function at T -+ 0 K) allowing: (i) the cross-checking of the model used, owing to the more reliable experimental data of yt;
(ii) the suggestion of a starting point with respect to the interaction models to be used during a lattice dynamical calculation and in the study of anharmonic properties, taking the sensitivity, e.g. of y(T), to the potential used within the low temperature region; and (iii) the separation of different contributions to y,,, particularly y2 and y3 [see equation (3)], whose relative values appear to be characteristic of the structure considered and may therefore allow for a check in the prediction of y,, values for compounds with no experimental data available. It is worth taking into account that this method for the evaluation of y,, avoids the integration over the reciprocal space that is very often a source of errors in these calculations (Barron & Batana, 1968; Barron et al., 1980). Acknowledgements-The authors gratefully thank Dr T. H. K. Barron for his very helpful comments and dis-
cussions. The authors are grateful to the referees for their valuable comments, which have been found very useful for revising the manuscript. We also acknowledge financial support from Universidad de Buenos Aires and CONICET.
REFERENCES Akima H. (1974) Comm. ACM 17(l), 26. Barron T. H. K. (in press) Generalised Theory of Thermal Expansion of Solids, Chapter 1 of Vol. I-4. CINDAS Series on Thermophysical Properties. Barron T. H. K. & Batana A. (1968) Whys. Reu. 167, 814. Barron T. H. K., Collins J. G. & White G. K. (1980) Adv. Phys. 29, 609. Cowley R. A. (1962) Proc. Roy. Sot. A268, 121. Cowley R. A., .&&ran W., Brockhouse B. N. 8r Woods A. D. B. (1963) Phvs. Rev. 131. 1030. Dutt N., Agrawal’G. G. & Shanker J. (1985) Phys. Stat. Sol. (b) 132, 99. de Launay J. (1956) J. Chem. Phys. 24, 1071. de Launay J. (1959) J. Chem. Phys. 30, 91. Lewis J. T., Lehozky A. & Briswe C. V. (1967) Phys. Rev. 161, 877. Namjoshi K. V., Mitra S. S. KLVetelino J. F. (1971) Phys. Rev. B 12, 4398. Soriano M. R. & Batana A. (in press) An. Asoc. Quim. Arg. Woods A. D. B., Cochran W. & Brockhouse B. N. (1962) Phys. Rev. 119, 980.