- Email: [email protected]
“Force theorem” and elastic constants of solids
Solid State Communications, Vol. 49, No. 7, pp. 701-705, 1984. Printed in Great Britain.
0038-1098/84 $3.00 + .00 Pergamon Press Ltd.
"FORCE THEOREM" AND ELASTIC CONSTANTS OF SOLIDS N.E. Christensen Max-Planck-Institut for Festk6rperforschung, D-7000 Stuttgart-80, West Germany and Physics Laboratory I, The Technical University of Denmark, DK-2800 Lyngby, Denmark
(Received 27 October 1983 by M. Cardona) The so-called "Force Theorem" derived by Andersen describes how the change in the total energy of an electron system can be calculated to first order in a virtual displacement. It is demonstrated here, how this theorem can be applied in calculation of elastic constants of crystals although this in fact requires calculation of total-energy changes accurate to second order in the distortion parameter.
1. INTRODUCTION IN THE LOCAL-DENSITY (LD) approximation, the change in total energy of the system, e.g. a solid or molecule, following a small change, a displacement of a part of the system, can to first-order be calculated as 6 u = ~ ~ . E . + ~Ue~.
(1)
The first term in equation (1) is the change in oneparticle energy calculated for a virtual displacement that is restricted in the following sense. The electronic structure (band structure in a solid) is calculated selfconsistently for the undistorted system using the localdensity scheme. The part of the system (a "cell" in a solid) which has to be moved is virtually cut loose and moved to its new position. The potential inside the "cell" is kept frozen, i.e. equal to the potential which was obtained in the first (self-consistent) calculation. Using these potential functions the electronic structure is again calculated - with the "cell" moved to the new position. This second calculation is not self-consistent. The difference between the new one-particle energy sum and that obtained from the first calculation gives the first term on the r.h.s, of equation (1). The second term, fi Uel, is the change in (classical) electrostatic interaction energy between the "cell" which has been cut out and the rest of the system following this restricted distortion. This theorem, the Force Theorem (FT in the following) was derived by Andersen [1 ] and it was, for instance, used [1] to deduce the pressure relations of Pettifor [2] in a simple way. These relations follow [1 ] in the Atomic Sphere Approximation [3] (ASA) readily from equation (1) for a monatomic solid if the "cells" are taken as (neutral) atomic spheres, i.e. the second term is zero. Skriver [4] recently used equation (1) to 701
calculate structural energy differences for the alkaline earth metals assuming that the second term, ~ Uel, could be neglected. As stressed by Heine [1 ] the derivation of the FT [1 ] helps to explain qualitatively the success in using the one-electron energy sum only to calculate total energy changes. Friedel explained [5] the approximate parabolic behaviour of the cohesive energies through the transition-metal series from such an approximation. Further, this approximation, in connection with tight-binding models, has been successful in explaining anomalies in temperature-dependence of elastic constants and in phonon spectra [6-8]. It is the purpose of the present communication to follow a suggestion by Andersen [9], and explicitly to demonstrate that the FT also can be applied to calculate elastic shear constants, although this in fact requires energy changes to second-order in the deformation parameter. It will further be shown that in this case the cells cannot be considered as (neutral) spheres with a spherically symmetric charge distribution. The second term plays an extremely important role and must be calculated in a model that takes the non-sphericity of the charge and its change to first-order in the deformation into account. Details of the calculations will not be presented here - only the structure of some of the important relations will be discussed. A few numerical results (tetragonal shear constants for Pd and Au) will be given. A detailed discussion together with more extensive numerical results will be presented elsewhere [10].
2. ELASTIC CONSTANTS FROM FORCE THEOREM A shear of a crystal causes a point originally at r to be moved to
r' = I
• r.
(2)
=
If a volume conserving, tetragonal shear is taken as an example - as we do in the actual calculations presented here - the deformation matrix can be expressed
=
I
'
i
I
~
1
0.0
as =
Vol. 49, No. 7
"FORCE THEOREM" AND ELASTIC CONSTANTS OF SOLIDS
702
-0.1
e 7/2 e -~/2 0
001 .
(3)
hE v
~-0.2 uJ
We could in principle calculate the elastic constant c [in this case c ' = ½(Cll --cl2)] as c -
-0.3
1 a2U
(4)
3 a-? '
where the total energy, U, of the crystal is calculated as a function of 3, by means of self-consistent band calculations for a series o f sheared crystals corresponding to different "/-values. We wish to show, however, that only one self-consistent calculation is necessary, and further that we can do this in the simplest case, namely for the unsheared crystal and calculate c from the FT. This procedure avoids complicated - and numerically less accurate - calculations o f electron-electron interactions and their derivatives with respect to 7. Although the fundamental relations are simple in principle, we shall discuss them in some detail. Equation (4) is rewritten as c
1 a2u
1 aF
3 a-/2
3 a-/ '
(5)
-
where the "force" F is F -
aU a-/'
(5a)
FT specifies how to calculate the "force", i.e. a U t o first-order in % Assume that a specific 3, (y = 7o), i.e. a specific lattice shear, has been chosen, then the "force" corresponding to this 7o is
[ aE~ + aU~ ]
F(-/o) = [ 0-/
O-/
(6)
.r=%'
where E ~ ("band structure term") is the one-particle energy sum: EF(~)
Eb6(-/) =
¢O
g ~
N.y(E) EdE.
(6a)
As mentioned previously, the FT assumes that aEbs/a-/]~=,/, is calculated by taking the difference between Ebs(-/o + A-/) and Ebs(-/o) where the potentials in the cells of the crystal sheared to 7o + A-/are the potentials obtained by a self-consistent (LD) band
0./0/~
Tetragonal shear
~
I 2
.Au oPd
,
, 0 •y
I -2
,
1 -4 xlO -2
Fig. 1. Change in one-particle energy (Ebs) as a function o f deformation parameter, -/, for f.c.c. Pd and Au. Tetragonal shear. Circles indicate the calculated values. The curves are parabolas extracted from a polynomial least-squares fit. calculation at 70. If we choose 70 = 0, and calculate fiEbs as a function of A7 we obtain - for small A-/-- a parabola. Actual numerical results are shown for Pd and Au (tetragonal deformation) in Fig. 1. If we choose another value for 7o; i.e. we perform the self-consistent calculation for a sheared lattice and derive a new fiEbs relation, actual numerical calculation have shown that this (70 chosen to 2%) to high degree of precision coincides with the curves obtained by starting from the unsheared crystal. A priori, this can only be shown [9] to be the case if the muffin-tin, the ASA - or a similar approximation for the one-electron potential is used. Thus, as far as the first term in equation (6) is concerned, and its contribution to the elastic constant, it is concluded that the self-consistent calculation can be carried out at -/o = 0. Palladium and gold have on purpose been chosen as examples. The one-particle energy changes shown in Fig. 1 are negative. Therefore, if only the contributions from the one-electron energies are taken into account, i.e. 6Uea in equation (1) is neglected, then both crystals are predicted to be unstable (c' is negative). The calculation of the electrostatic contribution has given the most practical difficulties. In principle it is easy to describe, however, how it can be done. The "cells" which are moved during the restricted virtual displacement are the Wigner-Seitz cells. Although these, for a monoatomic solid, are neutral objects, they
"FORCE THEOREM" AND ELASTIC CONSTANTS OF SOLIDS
Vol. 49; No. 7
interact electrostatically due to their non-spherical shape and due to the fact that they contain a non-spherical charge distribution. The electrostatic interaction energy [9, 10] between a selected cell at R = 0 and the rest of the lattice, assuming that the shear is 3'0, is raral ral ra-ra' (3'o), u (3'o) = Z All' Ql' (3'0) Q~n(3'0) S~+r
(7)
If
rara+
where slm(3'o) = ~] ~ o 4=0
Y~(R%) Rl+l ,
(7a)
7o
is a structure factor, 1~o being the lattice vectors of the (sheared) crystal. The Q's are the multipole-moments [10] of the charge distribution in a cell, and the coefficients A are given in [10]. (Below Q~n, Qrm'An'ram'and Sin+i,m' will be denoted by Qi, Qj, AO and SO). The electrostatic "force" at the shear 3'0 is Fea(7o) = {[Uel(7o + A T ) - Uea(7o)]/A3'}zx~o,
u
AoQ+(7o) Qj(7o)
(8)
[ 33' b: o'
Qi(3'o)Qj(3'o) ~ Qi(0)Qj(0)
(+°
[Q/(3') Qj(3')]
/
,
3,=0
-,.:+,.o - I 3+, j.,.:o + 3'° l and the "force" is, to first-order in 3"0:
1 t
+ Q Qj
(9)
The derivation of equation (9) has been explained in detail in order to avoid a pitfall that would lead to an electrostatic contribution given by the first term in equation (9) alone. Actual calculations show that the second term is extremely important; ~(Qi Qj)/a3"bsiJ/a3" is larger in magnitude than QiQj a2StJ/33"2 and has the opposite sign in the cases which we have considered. The electrostatic term cannot be related to a shear-derivative of a Madelung energy for point charges in a uniform background. An approach reminiscent of the one presented here, and suggested earlier by Andersen [1,9], has been presented by Ashkenazi etal. [12, 13]. As we do not understand their method of estimating the electrostatic interaction term we refrain here from a detailed comparison of our results to those of [ 12, 13 ].
The elastic constant corresponding to a volume conserving shear is calculated by means of the FT as c = c ~ + eel
1 =o
32SiJ I , 33' ~3,=o
where the first term is Fea(0). This is zero, and the electrostatic contribution to the elastic constant is
(10)
where cea is given in equation (9), and 1 [aZE~(7)]
[3S~(7)]
can then be calculated, equations (7)-(8), as a function of 3'o. For small 3'0, Fea(7o) is proportional to 3'o. The "force" is zero for 3'o = 0. For small values of 3'o we can use the following expansions
+ 3'0
= 1~, { --a2SiJ+ a(QiQj) SSiJt 3 ii Aij QiQj a3"2 ~3" -~7 ),),=o"
3. NUMERICAL EXAMPLES
where the moments which are used to calculate Uea(3'0 + ~3') and Uea(3'o) are the same, namely those corresponding to the charge distribution calculated from the (self-consistent) band calculation for the lattice sheared as specified by the mode and by 3'0. The electrostatic "force" Fea(7o) = Z
Cel
703
(10a)
a3': J :o' with E~(3') given by equation (6a). Self-consistent potentials were generated by means of the [3] LMTO (Linear Muffin Tin Orbital) method, and 8Eta(3') was then from these potentials calculated (including the combined correction term [3]) from one LMTO calculation at each shear, 3'. For Pd and Au (tetragonal shear) the 5E~-curves have already been discussed (Fig. 1). The values obtained for ct~ are given in Table 1. The calculation of cea [equation (9)] caused considerable difficulty at the initial stages of the work, because we did not expect the second term [equation (9)] to be significant, and because it turned out that also the deviation from sphericallity of the charge density and its shear-dependence must be included in the calculation. Figure 2 shows for Pd how different approximations affect the results obtained by equation (9). The results are plotted against the maximum angular momentum, lmax, included, i.e. the maximum values of l and l' included in equation (7). The points in Fig. 2 connected with a dotted line were obtained from equation (9) assuming a "muff'm-tin" charge distribution in the cells, i.e. a density which is spherically symmetric inside the inscribed sphere and constant in the interstitial region. The dashed line corresponds to calculations where the muffin-density
704
"FORCE THEOREM" AND ELASTIC CONSTANTS OF SOLIDS
Vol. 40, No. 7
Table 1. Calculated and experimental tetragonal-shear constants for Pd and Au. The four entries for Pd correspond to different approximations used in calculation o f the electrostatic term, c'el [equation (9)]: (MT-MT): Muffin-tin charge density used in Q and ~Q/37. (Ns-MT): Non-spherical charge included in Q, but not in OQ/~7. (Ns-Ns): Non-spherical density used for calculations o f Q as well as OQ/O7. The value o f c'eagiven in the column with the additional label "cubic" was obtained keeping the non-spherical charge fixed (i.e. with cubic symmetry) under shear. The total shear constant Ctot' is the sum of the band-structure term, Cbs, and the electrostatic contribution c'ea.All values are in Mbar Pd(MT-MT)
Pd(Ns-MT)
0.74 0.65
0.55 0.46
1
ebs %1 P Ctot C'exp [14]
-
-
Pd(Ns-Ns) (cubic)
Pd(Ns-Ns)
0.65 0.56
0.36 0.27 0.29
0.09
Au(Ns-Ns)
-
-
O. 16 0.36 0.20 0.16
4. DISCUSSION AND CONCLUSIONS
contribution must be included, and (3) the nonsphericity of the charge distribution and its variation with shear must be taken into account. It was mentioned earlier that the second term in equation (9) is important. If we, in the case of gold, omit this term, the calculated value of c' is -- 1.24 Mbar (!). Admittedly, the present model [11 ] can be further refined - and using it for calculations for all [10] transition metals we expect that larger differences between theory and experiment may be found. However, the calculations are sufficiently reliable to demonstrate the applicability of the FT. Further, the success of the present scheme for calculation of elastic constants suggests that a similar method may be applied to cases of more complicated deformations (e.g. phonons) and of more physical interest. Some of the advantages of calculating total-energy differences by means of the FT instead of performing a series of self-consistent calculations have been mentioned (Section 2). These are in particular obvious in the cases where the ASA approximation is sufficient, and the choice of cutting along neutral spheres implies that the electrostatic interaction term vanishes. It might be argued that, since this term cannot be neglected in our case, although we still displace neutral cells, the advantages of the FT are less clear, or we might have chosen to move "cells" defined in a different way. If we, for instance, had decided to consider the "cells" as containing only the nuclei our calculation would correspond to an application of the Hellmann-Feynman theorem. This would, however, require very accurate calculations of the core states which, on the contrary, is not necessary when the cutting surface is chosen to lie in the regions of space where the charge density is low - the Wigner-Seitz cell.
The numerical results given in Table 1 show that: (1) The Force Theorem [1] can be used in first-principles calculations of elastic constants, (2) the electrostatic
Acknowledgements - The author is grateful to O.K. Andersen, who initiated this work, for supplying many of the basic ideas used, and to D. G16tzel for a good collaboration in the initial stages. W. Temmerman is
1.0
'
I
~
'
I
'
I
Pd
E O
-~
0.5
~.ii.." ."
iv"
•
",
.. ..'~
U
0
, 0
I ~
I
I
8
Imclx
,
I 12
Fig. 2. Electrostatic contribution to the tetragonal shear constant in Pd [equation (9)], plotted against the maximum angular momentum lmax of the multipole moments [equation (7)]. Dotted line: Muffin-tin charge density. Dashed line: Muffin-tin density in calculation of ~Q/~7, non-spherical density in the calculation of Q. Full line: Non-spherical charge in the Q as well as in the aQ/a7 calculations. was used in the moment derivatives [the last term of equation (9)] whereas the Q's themselves were calculated with the non-spherical charge distribution. The results illustrated in Fig. 2 with the solid lines were obtained with the moments and their shear-derivatives all calculated for non-spherical charge distributions. The actual methods for the calculation of these quantities will be described elsewhere [10, 11 ]. The three values obtained for lrnax = 12 are given, for comparison, in Table 1.
Vol. 49, No. 7
"FORCE THEOREM" AND ELASTIC CONSTANTS OF SOLIDS
thanked for discussions concerning the calculation of the non-spherical charge.
9. 10. 11.
REFERENCES 1.
2. 3. 4. 5. 6. 7. 8.
A.R. Mackintosh & O.K. Andersen, Electrons at the Fermi Surface (Edited by M. Springford), Section 5.3.11. Cambridge University Press, Cambridge (1979). For a further discussion, see V. Heine, Solid State Phys. 35, 1 (1980) Section 16.b. D.G. Pettifor, Commun. Phys. 1,141 (1976). O.K. Andersen, Phys. Rev. BI2, 3060 (1975). H.L. Skriver, Phys. Rev. Lett. 49, 1768 (1982). J. Friedel, The Physics of Metals, (Edited by J.M. Ziman), Cambridge University Press (1969). J. Ashkenazi, M.- Dacorogna, M. Peter, Y. Talmor, E. Walter & S. Steinemann, Phys. Rev. BI8, 4120 (1978). P. Bujard, R. Sanjines, E. Walker, J. Ashkenazi & M. Peter, J. Phys. FI1,775 (1981). C.M. Varma & W. Weber, Phys. Rev. BI9, 6142 (1979).
12. 13. 14.
705
O.K. Andersen (private communication); see also V. Heine [ 1]. N.E. Christensen & D. G16tzel (to be published). N.E. Christensen (to be published) and O.K. Andersen (private communication). The nonspherical charge distribution was obtained from the eigenvectors of the LMTO calculation. Outside the inscribed sphere, the wavefunctions are made up of pseudo-MTO's (see [3]), and it is assumed that, for s- and p-states, the non-spherical part of the charge density also inside the spheres could be represented sufficiently accurately by the contributions from these functions. The d-states inside the sphere are treated separately. J. Ashkenazi, Phys. Rev. B2, 1512 (1982). M. Dacorogna, J. Ashkenazi & M. Peter, Phys. Rev. B26, 1527 (1982). G. Simmons & H. Wang, Single Crystal Elastic Constants and Calculated Aggregate Properties: A. Handbook, MIT Press, Cambridge, Massachussetts (1971).
0038-1098/84 $3.00 + .00 Pergamon Press Ltd.
"FORCE THEOREM" AND ELASTIC CONSTANTS OF SOLIDS N.E. Christensen Max-Planck-Institut for Festk6rperforschung, D-7000 Stuttgart-80, West Germany and Physics Laboratory I, The Technical University of Denmark, DK-2800 Lyngby, Denmark
(Received 27 October 1983 by M. Cardona) The so-called "Force Theorem" derived by Andersen describes how the change in the total energy of an electron system can be calculated to first order in a virtual displacement. It is demonstrated here, how this theorem can be applied in calculation of elastic constants of crystals although this in fact requires calculation of total-energy changes accurate to second order in the distortion parameter.
1. INTRODUCTION IN THE LOCAL-DENSITY (LD) approximation, the change in total energy of the system, e.g. a solid or molecule, following a small change, a displacement of a part of the system, can to first-order be calculated as 6 u = ~ ~ . E . + ~Ue~.
(1)
The first term in equation (1) is the change in oneparticle energy calculated for a virtual displacement that is restricted in the following sense. The electronic structure (band structure in a solid) is calculated selfconsistently for the undistorted system using the localdensity scheme. The part of the system (a "cell" in a solid) which has to be moved is virtually cut loose and moved to its new position. The potential inside the "cell" is kept frozen, i.e. equal to the potential which was obtained in the first (self-consistent) calculation. Using these potential functions the electronic structure is again calculated - with the "cell" moved to the new position. This second calculation is not self-consistent. The difference between the new one-particle energy sum and that obtained from the first calculation gives the first term on the r.h.s, of equation (1). The second term, fi Uel, is the change in (classical) electrostatic interaction energy between the "cell" which has been cut out and the rest of the system following this restricted distortion. This theorem, the Force Theorem (FT in the following) was derived by Andersen [1 ] and it was, for instance, used [1] to deduce the pressure relations of Pettifor [2] in a simple way. These relations follow [1 ] in the Atomic Sphere Approximation [3] (ASA) readily from equation (1) for a monatomic solid if the "cells" are taken as (neutral) atomic spheres, i.e. the second term is zero. Skriver [4] recently used equation (1) to 701
calculate structural energy differences for the alkaline earth metals assuming that the second term, ~ Uel, could be neglected. As stressed by Heine [1 ] the derivation of the FT [1 ] helps to explain qualitatively the success in using the one-electron energy sum only to calculate total energy changes. Friedel explained [5] the approximate parabolic behaviour of the cohesive energies through the transition-metal series from such an approximation. Further, this approximation, in connection with tight-binding models, has been successful in explaining anomalies in temperature-dependence of elastic constants and in phonon spectra [6-8]. It is the purpose of the present communication to follow a suggestion by Andersen [9], and explicitly to demonstrate that the FT also can be applied to calculate elastic shear constants, although this in fact requires energy changes to second-order in the deformation parameter. It will further be shown that in this case the cells cannot be considered as (neutral) spheres with a spherically symmetric charge distribution. The second term plays an extremely important role and must be calculated in a model that takes the non-sphericity of the charge and its change to first-order in the deformation into account. Details of the calculations will not be presented here - only the structure of some of the important relations will be discussed. A few numerical results (tetragonal shear constants for Pd and Au) will be given. A detailed discussion together with more extensive numerical results will be presented elsewhere [10].
2. ELASTIC CONSTANTS FROM FORCE THEOREM A shear of a crystal causes a point originally at r to be moved to
r' = I
• r.
(2)
=
If a volume conserving, tetragonal shear is taken as an example - as we do in the actual calculations presented here - the deformation matrix can be expressed
=
I
'
i
I
~
1
0.0
as =
Vol. 49, No. 7
"FORCE THEOREM" AND ELASTIC CONSTANTS OF SOLIDS
702
-0.1
e 7/2 e -~/2 0
001 .
(3)
hE v
~-0.2 uJ
We could in principle calculate the elastic constant c [in this case c ' = ½(Cll --cl2)] as c -
-0.3
1 a2U
(4)
3 a-? '
where the total energy, U, of the crystal is calculated as a function of 3, by means of self-consistent band calculations for a series o f sheared crystals corresponding to different "/-values. We wish to show, however, that only one self-consistent calculation is necessary, and further that we can do this in the simplest case, namely for the unsheared crystal and calculate c from the FT. This procedure avoids complicated - and numerically less accurate - calculations o f electron-electron interactions and their derivatives with respect to 7. Although the fundamental relations are simple in principle, we shall discuss them in some detail. Equation (4) is rewritten as c
1 a2u
1 aF
3 a-/2
3 a-/ '
(5)
-
where the "force" F is F -
aU a-/'
(5a)
FT specifies how to calculate the "force", i.e. a U t o first-order in % Assume that a specific 3, (y = 7o), i.e. a specific lattice shear, has been chosen, then the "force" corresponding to this 7o is
[ aE~ + aU~ ]
F(-/o) = [ 0-/
O-/
(6)
.r=%'
where E ~ ("band structure term") is the one-particle energy sum: EF(~)
Eb6(-/) =
¢O
g ~
N.y(E) EdE.
(6a)
As mentioned previously, the FT assumes that aEbs/a-/]~=,/, is calculated by taking the difference between Ebs(-/o + A-/) and Ebs(-/o) where the potentials in the cells of the crystal sheared to 7o + A-/are the potentials obtained by a self-consistent (LD) band
0./0/~
Tetragonal shear
~
I 2
.Au oPd
,
, 0 •y
I -2
,
1 -4 xlO -2
Fig. 1. Change in one-particle energy (Ebs) as a function o f deformation parameter, -/, for f.c.c. Pd and Au. Tetragonal shear. Circles indicate the calculated values. The curves are parabolas extracted from a polynomial least-squares fit. calculation at 70. If we choose 70 = 0, and calculate fiEbs as a function of A7 we obtain - for small A-/-- a parabola. Actual numerical results are shown for Pd and Au (tetragonal deformation) in Fig. 1. If we choose another value for 7o; i.e. we perform the self-consistent calculation for a sheared lattice and derive a new fiEbs relation, actual numerical calculation have shown that this (70 chosen to 2%) to high degree of precision coincides with the curves obtained by starting from the unsheared crystal. A priori, this can only be shown [9] to be the case if the muffin-tin, the ASA - or a similar approximation for the one-electron potential is used. Thus, as far as the first term in equation (6) is concerned, and its contribution to the elastic constant, it is concluded that the self-consistent calculation can be carried out at -/o = 0. Palladium and gold have on purpose been chosen as examples. The one-particle energy changes shown in Fig. 1 are negative. Therefore, if only the contributions from the one-electron energies are taken into account, i.e. 6Uea in equation (1) is neglected, then both crystals are predicted to be unstable (c' is negative). The calculation of the electrostatic contribution has given the most practical difficulties. In principle it is easy to describe, however, how it can be done. The "cells" which are moved during the restricted virtual displacement are the Wigner-Seitz cells. Although these, for a monoatomic solid, are neutral objects, they
"FORCE THEOREM" AND ELASTIC CONSTANTS OF SOLIDS
Vol. 49; No. 7
interact electrostatically due to their non-spherical shape and due to the fact that they contain a non-spherical charge distribution. The electrostatic interaction energy [9, 10] between a selected cell at R = 0 and the rest of the lattice, assuming that the shear is 3'0, is raral ral ra-ra' (3'o), u (3'o) = Z All' Ql' (3'0) Q~n(3'0) S~+r
(7)
If
rara+
where slm(3'o) = ~] ~ o 4=0
Y~(R%) Rl+l ,
(7a)
7o
is a structure factor, 1~o being the lattice vectors of the (sheared) crystal. The Q's are the multipole-moments [10] of the charge distribution in a cell, and the coefficients A are given in [10]. (Below Q~n, Qrm'An'ram'and Sin+i,m' will be denoted by Qi, Qj, AO and SO). The electrostatic "force" at the shear 3'0 is Fea(7o) = {[Uel(7o + A T ) - Uea(7o)]/A3'}zx~o,
u
AoQ+(7o) Qj(7o)
(8)
[ 33' b: o'
Qi(3'o)Qj(3'o) ~ Qi(0)Qj(0)
(+°
[Q/(3') Qj(3')]
/
,
3,=0
-,.:+,.o - I 3+, j.,.:o + 3'° l and the "force" is, to first-order in 3"0:
1 t
+ Q Qj
(9)
The derivation of equation (9) has been explained in detail in order to avoid a pitfall that would lead to an electrostatic contribution given by the first term in equation (9) alone. Actual calculations show that the second term is extremely important; ~(Qi Qj)/a3"bsiJ/a3" is larger in magnitude than QiQj a2StJ/33"2 and has the opposite sign in the cases which we have considered. The electrostatic term cannot be related to a shear-derivative of a Madelung energy for point charges in a uniform background. An approach reminiscent of the one presented here, and suggested earlier by Andersen [1,9], has been presented by Ashkenazi etal. [12, 13]. As we do not understand their method of estimating the electrostatic interaction term we refrain here from a detailed comparison of our results to those of [ 12, 13 ].
The elastic constant corresponding to a volume conserving shear is calculated by means of the FT as c = c ~ + eel
1 =o
32SiJ I , 33' ~3,=o
where the first term is Fea(0). This is zero, and the electrostatic contribution to the elastic constant is
(10)
where cea is given in equation (9), and 1 [aZE~(7)]
[3S~(7)]
can then be calculated, equations (7)-(8), as a function of 3'o. For small 3'0, Fea(7o) is proportional to 3'o. The "force" is zero for 3'o = 0. For small values of 3'o we can use the following expansions
+ 3'0
= 1~, { --a2SiJ+ a(QiQj) SSiJt 3 ii Aij QiQj a3"2 ~3" -~7 ),),=o"
3. NUMERICAL EXAMPLES
where the moments which are used to calculate Uea(3'0 + ~3') and Uea(3'o) are the same, namely those corresponding to the charge distribution calculated from the (self-consistent) band calculation for the lattice sheared as specified by the mode and by 3'0. The electrostatic "force" Fea(7o) = Z
Cel
703
(10a)
a3': J :o' with E~(3') given by equation (6a). Self-consistent potentials were generated by means of the [3] LMTO (Linear Muffin Tin Orbital) method, and 8Eta(3') was then from these potentials calculated (including the combined correction term [3]) from one LMTO calculation at each shear, 3'. For Pd and Au (tetragonal shear) the 5E~-curves have already been discussed (Fig. 1). The values obtained for ct~ are given in Table 1. The calculation of cea [equation (9)] caused considerable difficulty at the initial stages of the work, because we did not expect the second term [equation (9)] to be significant, and because it turned out that also the deviation from sphericallity of the charge density and its shear-dependence must be included in the calculation. Figure 2 shows for Pd how different approximations affect the results obtained by equation (9). The results are plotted against the maximum angular momentum, lmax, included, i.e. the maximum values of l and l' included in equation (7). The points in Fig. 2 connected with a dotted line were obtained from equation (9) assuming a "muff'm-tin" charge distribution in the cells, i.e. a density which is spherically symmetric inside the inscribed sphere and constant in the interstitial region. The dashed line corresponds to calculations where the muffin-density
704
"FORCE THEOREM" AND ELASTIC CONSTANTS OF SOLIDS
Vol. 40, No. 7
Table 1. Calculated and experimental tetragonal-shear constants for Pd and Au. The four entries for Pd correspond to different approximations used in calculation o f the electrostatic term, c'el [equation (9)]: (MT-MT): Muffin-tin charge density used in Q and ~Q/37. (Ns-MT): Non-spherical charge included in Q, but not in OQ/~7. (Ns-Ns): Non-spherical density used for calculations o f Q as well as OQ/O7. The value o f c'eagiven in the column with the additional label "cubic" was obtained keeping the non-spherical charge fixed (i.e. with cubic symmetry) under shear. The total shear constant Ctot' is the sum of the band-structure term, Cbs, and the electrostatic contribution c'ea.All values are in Mbar Pd(MT-MT)
Pd(Ns-MT)
0.74 0.65
0.55 0.46
1
ebs %1 P Ctot C'exp [14]
-
-
Pd(Ns-Ns) (cubic)
Pd(Ns-Ns)
0.65 0.56
0.36 0.27 0.29
0.09
Au(Ns-Ns)
-
-
O. 16 0.36 0.20 0.16
4. DISCUSSION AND CONCLUSIONS
contribution must be included, and (3) the nonsphericity of the charge distribution and its variation with shear must be taken into account. It was mentioned earlier that the second term in equation (9) is important. If we, in the case of gold, omit this term, the calculated value of c' is -- 1.24 Mbar (!). Admittedly, the present model [11 ] can be further refined - and using it for calculations for all [10] transition metals we expect that larger differences between theory and experiment may be found. However, the calculations are sufficiently reliable to demonstrate the applicability of the FT. Further, the success of the present scheme for calculation of elastic constants suggests that a similar method may be applied to cases of more complicated deformations (e.g. phonons) and of more physical interest. Some of the advantages of calculating total-energy differences by means of the FT instead of performing a series of self-consistent calculations have been mentioned (Section 2). These are in particular obvious in the cases where the ASA approximation is sufficient, and the choice of cutting along neutral spheres implies that the electrostatic interaction term vanishes. It might be argued that, since this term cannot be neglected in our case, although we still displace neutral cells, the advantages of the FT are less clear, or we might have chosen to move "cells" defined in a different way. If we, for instance, had decided to consider the "cells" as containing only the nuclei our calculation would correspond to an application of the Hellmann-Feynman theorem. This would, however, require very accurate calculations of the core states which, on the contrary, is not necessary when the cutting surface is chosen to lie in the regions of space where the charge density is low - the Wigner-Seitz cell.
The numerical results given in Table 1 show that: (1) The Force Theorem [1] can be used in first-principles calculations of elastic constants, (2) the electrostatic
Acknowledgements - The author is grateful to O.K. Andersen, who initiated this work, for supplying many of the basic ideas used, and to D. G16tzel for a good collaboration in the initial stages. W. Temmerman is
1.0
'
I
~
'
I
'
I
Pd
E O
-~
0.5
~.ii.." ."
iv"
•
",
.. ..'~
U
0
, 0
I ~
I
I
8
Imclx
,
I 12
Fig. 2. Electrostatic contribution to the tetragonal shear constant in Pd [equation (9)], plotted against the maximum angular momentum lmax of the multipole moments [equation (7)]. Dotted line: Muffin-tin charge density. Dashed line: Muffin-tin density in calculation of ~Q/~7, non-spherical density in the calculation of Q. Full line: Non-spherical charge in the Q as well as in the aQ/a7 calculations. was used in the moment derivatives [the last term of equation (9)] whereas the Q's themselves were calculated with the non-spherical charge distribution. The results illustrated in Fig. 2 with the solid lines were obtained with the moments and their shear-derivatives all calculated for non-spherical charge distributions. The actual methods for the calculation of these quantities will be described elsewhere [10, 11 ]. The three values obtained for lrnax = 12 are given, for comparison, in Table 1.
Vol. 49, No. 7
"FORCE THEOREM" AND ELASTIC CONSTANTS OF SOLIDS
thanked for discussions concerning the calculation of the non-spherical charge.
9. 10. 11.
REFERENCES 1.
2. 3. 4. 5. 6. 7. 8.
A.R. Mackintosh & O.K. Andersen, Electrons at the Fermi Surface (Edited by M. Springford), Section 5.3.11. Cambridge University Press, Cambridge (1979). For a further discussion, see V. Heine, Solid State Phys. 35, 1 (1980) Section 16.b. D.G. Pettifor, Commun. Phys. 1,141 (1976). O.K. Andersen, Phys. Rev. BI2, 3060 (1975). H.L. Skriver, Phys. Rev. Lett. 49, 1768 (1982). J. Friedel, The Physics of Metals, (Edited by J.M. Ziman), Cambridge University Press (1969). J. Ashkenazi, M.- Dacorogna, M. Peter, Y. Talmor, E. Walter & S. Steinemann, Phys. Rev. BI8, 4120 (1978). P. Bujard, R. Sanjines, E. Walker, J. Ashkenazi & M. Peter, J. Phys. FI1,775 (1981). C.M. Varma & W. Weber, Phys. Rev. BI9, 6142 (1979).
12. 13. 14.
705
O.K. Andersen (private communication); see also V. Heine [ 1]. N.E. Christensen & D. G16tzel (to be published). N.E. Christensen (to be published) and O.K. Andersen (private communication). The nonspherical charge distribution was obtained from the eigenvectors of the LMTO calculation. Outside the inscribed sphere, the wavefunctions are made up of pseudo-MTO's (see [3]), and it is assumed that, for s- and p-states, the non-spherical part of the charge density also inside the spheres could be represented sufficiently accurately by the contributions from these functions. The d-states inside the sphere are treated separately. J. Ashkenazi, Phys. Rev. B2, 1512 (1982). M. Dacorogna, J. Ashkenazi & M. Peter, Phys. Rev. B26, 1527 (1982). G. Simmons & H. Wang, Single Crystal Elastic Constants and Calculated Aggregate Properties: A. Handbook, MIT Press, Cambridge, Massachussetts (1971).