PHYSICAL REVIEW B 67, 094301 共2003兲

Thermodynamic quantities of metals investigated by an analytic statistical moment method

K. Masuda-Jindo,1 Vu Van Hung,2 and Pham Dinh Tam2
1
Department of Materials Science and Engineering, Tokyo Institute of Technology, Nagatsuta 4259, Midori-ku,
Yokohama 226-8503, Japan
2
Department of Physics, Hanoi National Pedagogic University, km8 Hanoi-Sontay Highway, Hanoi, Vietnam
共Received 18 June 2002; published 5 March 2003兲
The thermodynamic properties of metals are studied by including explicitly the anharmonic effects of the
lattice vibrations going beyond the quasiharmonic approximations. The free energy, thermal lattice expansion
coefficients, mean-square atomic displacements, and specific heats at the constant volume and those at the
constant pressure, C v and C p , are derived in closed analytic forms in terms of the power moments of the
atomic displacements. The analytical formulas give highly accurate values of the thermodynamic quantities,
which are comparable to those of the molecular dynamics or Monte Carlo simulations for a wide temperature
range. The present formalism is well suited to calculate the thermodynamic quantities of metals and alloys by
including the many body electronic effects and by combining it with the first-principles approaches.

DOI: 10.1103/PhysRevB.67.094301 PACS number共s兲: 64.10.⫹h, 65.40.⫺b

I. INTRODUCTION for genetic studies of trends among a given class of metallic

materials. Therefore, they do not account for mostly impor-
The first-principles determination of the thermodynamic tant many-body electronic effects in metallic systems, and
quantities of metals and alloys are now of great importance they cannot be relied on for properties of real materials.
for the understanding of structural phase transformations as A number of theoretical approaches have been proposed
well as for the phase diagrams computations.1–3 So far, the to overcome the limitations of the QH theories. The first
first-principles density functional theories4 – 8 have been used calculation of the lowest-order anharmonic contributions to
extensively for the calculations of the ground state properties the atomic mean-square displacement 具 u 2 典 or the Debye-
of various metal systems at the absolute zero temperature. In Waller factor was done by Maradudin and Flinn18 in the
the phase transformations occurring in metals and alloys at leading-term approximation for a nearest-neighbor central-
finite temperatures 共under pressure P兲, the thermal lattice vi- force model. Since then, many anharmonic calculations in-
brations 共anharmonicity effects兲 play an essentially important cluding the lowest-order anharmonic contributions have been
role.9,10 However, most of the first-principles calculations for performed for metal systems.19,20 The method requires ac-
the structural phase transformations and alloy phase diagram knowledge of a number of Brillouin-zone sums14 and the
computations have been done with the use of the lattice vi- calculations are performed for the central-force model crys-
bration theory in the quasiharmonic 共QH兲 approx- tals. Recently, some attempts have been made to take into
imation.11–15 For the alloy phase diagram calculations, there account the bond length dependence of bond stiffness tensors
have been difficulties in accounting for the anharmonicity of in the calculations of the free energy of the substitutional
thermal lattice vibrations, especially for the higher tempera- alloys.21,22 The anharmonic effects of lattice vibrations on
ture region than the Debye temperature because the thermal the thermodynamic properties of the materials have also
lattice expansion plays an important role and cannot be ne- been studied by employing the first-order quantum-statistical
glected. The martensitic phase transformation in substitu- perturbation theory23–25 as well as by the first-order self-
tional alloys such as the Nix Al1⫺x system has also been stud- consistent 共SC兲 phonon theories.26 –31 The theories have been
ied with the QH approximation, and the temperature region used to analyze, e.g., the temperature dependence of ex-
treated by the QH theory is usually lower than the Debye tended x-ray absorption fine-structure 共EXAFS兲 spectra and
temperature.16 the phonon frequencies. However, the previous anharmonic-
The systems considered at high temperatures and high ity theories are still incomplete and have some inherent
pressures require the allowance for anharmonic effects which drawbacks and limitations.
are very essential in these regions. The simplest way is to use In the present study, we use the finite-temperature mo-
the QH Debye-Grüneisen theory.10 However, the results ob- ment expansion technique to derive the Helmholtz free ener-
tained in such a way are not always satisfactory. It is noted gies of metal systems, going beyond the QH approximations.
that the Debye form of the harmonic approximation is rather The thermodynamic quantities, mean-square atomic dis-
crude theory. The applicability of the QH method to the placements, specific heats, and elastic moduli are determined
study of particular metals is often restricted by the isotropic from the explicit expressions of the Helmholtz free energies.
Debye mode and the assumption of the mean sound velocity The Helmholtz free energy of the system at a given tempera-
v . 17 The temperature dependence of the lattice constant and ture T will be determined self-consistently with the equilib-
the linear thermal expansion coefficient are calculated by rium thermal lattice expansions of the crystal.
minimizing the free energy with respect to the volume of the We will use the electronic many-body potentials, i.e.,
system. Due to their simplicity, pair potentials are often used second-moment tight-binding 共TB兲 potentials,32– 40 for the

0163-1829/2003/67共9兲/094301共14兲/$20.00 67 094301-1 ©2003 The American Physical Society

K. MASUDA-JINDO, VU VAN HUNG, AND PHAM DINH TAM PHYSICAL REVIEW B 67, 094301 共2003兲

evaluation of the internal energy of the system. In metals the normal-mode frequencies 共phonon spectrum兲.48 This scheme
long-range Coulomb interaction and the partially filled va- is called as the QH approximation.
lence bands lead to interatomic forces that are inherently In the present study the thermodynamic quantities are cal-
many-body in nature. For more than a decade, the culated with the use of the electronic many-body potentials
embedded-atom method 共EAM兲41– 45 and the second-moment or the potentials derived by EAM. We note that the present
approximation 共SMA兲 of the TB scheme have been the two analytic formulation is quite useful when we combine it with
most common approaches, able to overcome the major limi- the ab initio theoretical scheme by numerically evaluating
tations of two-body pair potentials.18,46,47 The physical basis the harmonic k and anharmonic ␥ 1 and ␥ 2 parameters which
of EAM models makes them valid, especially for normal or will be defined in the subsequent derivations. The SMA TB
noble metals, whereas SMA is a priori well suited for tran- scheme is well suited to describe the cohesion of transition
sition elements 共with narrow d-band bonding兲. metals since they are elements with a partially filled narrow
In Sec. II, we will make a general derivation of the ther- d band superimposed on a broad free-electron-like s-p band.
mal lattice expansion and Helmholtz free energy of the The narrowness of the d band, especially in the 3d series, is
monoatomic cubic metals based on the fundamental prin- a consequence of the relative constriction of the d orbitals
ciples of quantum-statistical mechanics. The thermodynamic compared with the outer s and p orbitals. As one moves
quantities of the metals are then derived in terms of the across the periodic table, the d band is gradually being filled.
power moments of the atomic displacements from the Helm- Most of the properties of the transition metals are character-
holtz free energy of the system. Section III includes our main ized by the filling of the d band and ignoring the sp electrons.
calculation results of the thermodynamic quantities of some This constitutes Friedel’s d-band model which further as-
cubic metals. Finally, Sec. IV summarizes the present study. sumes a rectangular approximation for the density of states
␳ i (E) such that the bonding energy of the solid is primarily
due to the filling of the d band and proportional to its width.
II. THEORY
In the SMA, the bonding energy is then proportional to the
We derive the thermodynamic quantities of metals, taking root of the second moments 冑␮ (2) i . In metals, an important
into account the higher- 共fourth-兲 order anharmonic contribu- contribution to the structure comes from the repulsive term
tions in the thermal lattice vibrations going beyond the QH represented as a sum of pair potentials accounting for the
approximation. The basic equations for obtaining thermody- short-range behavior of the interaction between ions. There-
namic quantities of the given crystals are derived in a fol- fore, the cohesive energy of a transition metal consists of
lowing manner: The equilibrium thermal lattice expansions
are calculated by the force balance criterion and then the E coh⫽E rep⫹E bond . 共2兲
thermodynamic quantities are determined for the equilibrium
lattice spacings. The anharmonic contributions of the ther- The SMA has been used to suggest various functional
modynamic quantities are given explicitly in terms of the form for interatomic potentials in transition metals such as
power moments of the thermal atomic displacements. the Finnis-Sinclair potential,34 the closely related embedded
Let us first define the lattice displacements. We denote uil atom potential, and the TB SMA, also referred in the litera-
the vector defining the displacement of the ith atom, in the ture as to Gupta potential.33 The functional form we adopted
lth unit cell, from its equilibrium position. The potential en- here for elemental metals is that of the many-body SMA
ergy of the whole crystal U(uil ) is expressed in terms of the potential
positions of all the atoms from the sites of the equilibrium

兺冠 兺 冋 冉 冊册
lattice. We may assume that this function has a minimum N N
1 rij
when all the uil are zero, for the perfect lattice is presumably E ci ⫽ A exp ⫺ p ⫺1
a configuration of stable equilibrium. We use the theory of N i⫽1 j⫽i r0

再 兺 冋 冉 冊 册冎 冡
small atomic vibrations, and expand the potential energy U N 1/2
as a power series in the Cartesian components, u ilj , of the rij
⫺ ␰ 2i j exp ⫺2q ⫺1 , 共3兲
displacement vector uil around this point j⫽i r0

冋 册 冋 册
⳵U ⳵ 2U which has five parameters: ␧ 0 , ␰ i j 共for pure metals, ␰ i j
U⫽U 0 ⫹ 兺 ⳵ u ilj
uilj ⫹ 兺 j⬘
uilj ui ⬘⬘l ⬘
j
⫽ ␰ 0 ), p, q, and r 0 . The total cohesive energy E c of the
⬘ ⬘ ⳵ uil ⳵ ui ⬘ l ⬘
j
i,l, j
eq ⬘
ii ,l,ll , j j
eq
system is then written as the sum of the E ci . The parameters
A, ␰ 0 , p, and q are fitted to reproduce some experimental
⫹¯ , 共1兲 quantities at zero temperature 共cohesive energy E c , lattice
parameter a, bulk modulus and elastic constants兲. In the sum-
where U 0 denotes the internal 共cohesive兲 energy of the sys- mations over the index j in Eq. 共3兲 are either limited to the Z 1
tem. If we truncate the above expansion of Eq. 共1兲 up to the first neighbors only, and in that case we use the parameters
second-order terms, then the full interatomic potential is re- A, ␨ 0 , p, and q determined by Rosato, Guillope, and
placed by its quadratic expansion about the equilibrium Legrand,35 or extended up to the fifth neighbors, and in that
atomic positions. The system is then equivalent to a collec- case we use the parameters of Cleri and Rosato.36 Cleri and
tion of harmonic oscillators, and diagonalization of the cor- Rosato36 fitted these parameters to experimental data for 16
responding dynamical matrix yields the squares of the fcc and hexagonal-close-packed 共hcp兲 transition metals.

094301-2
THERMODYNAMIC QUANTITIES OF METALS . . . PHYSICAL REVIEW B 67, 094301 共2003兲

TABLE I. Parameters of the second moment TB potentials for cubic metals.

A 共eV兲 ␰ 共eV兲 p q E c 共eV/atom兲 a 共Å兲

Al共1兲a 0.1221 1.316 8.612 2.516 ⫺3.339 4.050

Al共2兲b 0.0334 0.7981 14.6147 1.112 ⫺3.339 4.050
Ni 0.1368 1.7558 10.00 2.70 ⫺4.435 3.523
Cu 0.0993 1.3543 10.08 2.56 ⫺3.544 3.615
Rh 0.0629 1.660 18.450 1.867 ⫺5.752 3.803
Pd 0.1746 1.718 10.867 3.742 ⫺3.936 3.887
Ag共1兲a 0.1028 1.1780 10.928 3.139 ⫺2.960 4.085
Ag共2兲b 0.1231 1.2811 10.12 3.37 ⫺2.960 4.085
Au 0.2061 1.790 10.229 4.036 ⫺3.779 4.079
Pt 0.2975 2.695 10.612 4.004 ⫺5.853 3.924
Li 0.0333 0.3249 7.75 0.737 ⫺1.63 3.49
Na 0.0159 0.2910 10.13 1.30 ⫺1.13 4.29
K 0.0205 0.2625 10.58 1.34 ⫺0.93 5.24
Rb 0.0194 0.2464 10.48 1.40 ⫺0.85 5.61
Cs 0.0205 0.2421 9.62 1.45 ⫺0.80 6.04
a
indicates parameters taken from Ref. 36.
b
indicates parameters taken from other sources: Al共2兲 from Ref. 60 and Ag共2兲 from Ref. 35.

The SMA TB potentials have been further extended and The supplementary forces ␣ i act in the direction of the gen-
revised not only for bulk metal systems but also for nanos- eralized coordinates q i . The thermodynamic quantities of the
cale materials. For Rh clusters, Chein, Blaston-Barojas, and harmonic crystal 共harmonic Hamiltonian兲 will be treated in
Pederson38 proposed the size-dependent parameters of the the Einstein approximation. In this respect, the present for-
SMA TB potentials, on the basis of their generalized gradient mulation is similar conceptually to the treatment of quantum
approximation 共GGA兲 calculations. A different parametriza- Monte Carlo method by Frenkel.52,53
tion strategy was introduced by Sigalas and After the action of the supplementary forces ␣ i the system
Papaconstantopoulos39 in which the parameters were fitted to passes into a new equilibrium state. For obtaining the statis-
local density approximation 共LDA兲 calculations of the total tical average of an thermodynamic quantity 具 q k 典 a for the
energy as a function of lattice constant. Li, Barojas, and new equilibrium state, we use the general formula for the
Papaconstantopoulos40 fitted the SMA TB potential param- correlation. Specifically, we use a recurrence formula54 based
eters to a LDA database that consists of the total energy as a on the density matrix in the quantum statistical mechanics
function of the lattice constant for both bcc and fcc lattices, 共for more details see Appendix A兲
rather than the fitting procedure to experimental quantities.
To simulate the long-range nature of the metallic bonding by
sp electrons in alkali metals, the interactions up to 12th- ⳵ 具 K̂ n 典 a
具 K̂ n⫹1 典 a ⫽ 具 K̂ n 典 a 具 q̂ n⫹1 典 a ⫹ ␪
neighbor shells 共228 atoms in bcc crystal兲 are taken into ⳵ ␣ n⫹1

冉 冊 冓 冔
account.40 Their potentials fitted to the first-principles LDA ⬁
B 2m iប 2m
⳵ K̂ 共n2m 兲
results are available for various metals, and more refined
nonorthogonal basis TB schemes39 are also proposed for the
⫺␪ 兺
m⫽0 共 2m 兲 ! ␪ ⳵ ␣ n⫹1
, 共5兲
a
quantitative calculations. The present thermodynamic formu-
lation is well suited to couple with any kind of TB schemes
where ␪ ⫽k B T, m is the atomic mass, and K̂ n is the correla-
mentioned above. The SMA TB potential parameters used in
tion operator of the nth order:
the present calculations are given in Table I.
We now consider a quantum system, which is influenced
by supplemental forces ␣ i in the space of the generalized 1
K̂ n ⫽ 关 ... 关 q̂ 1 ,q̂ 2 兴 ⫹ q̂ 3 兴 ⫹ ...] ⫹ q̂ n ] ⫹ . 共6兲
coordinates q i . 49–51 For simplicity, we only discuss mon- 2 n⫺1
atomic metallic systems, and hereafter omit the indices l on
the sublattices. Then, the Hamiltonian of the crystalline sys- In Eq. 共5兲 above, the symbol 具¯典 expresses the thermal av-
tem is given by eraging over the equilibrium ensemble with the Hamiltonian
Ĥ and B 2n denotes the Bernoulli numbers. 关 q i ,q j 兴 ⫹ repre-
Ĥ⫽Ĥ 0 ⫺ 兺i ␣ i q̂ i , 共4兲 sents the anticommutation relation. The general decoupling
formula of Eq. 共5兲 enables us to get all moments of the lattice
system and to investigate the nonlinear thermodynamic prop-
where Ĥ 0 denotes the crystalline Hamiltonian without the erties of the materials, taking into account the anharmonicity
supplementary forces ␣ i and the carets represent operators. of the thermal lattice vibrations. The Helmholtz free energy

094301-3
K. MASUDA-JINDO, VU VAN HUNG, AND PHAM DINH TAM PHYSICAL REVIEW B 67, 094301 共2003兲

of the system can then be obtained by taking into account the

higher-order moments 共up to fourth order兲.
The atomic force acting on a given ith atom in the lattice
␥ i⫽
1
6 冋冉 冊 冉
⳵ 4 E ci
⳵ u i4␣ eq
⫹6
⳵ 4 E ci
⳵ u i2␤ ⳵ u i2␥ 冊册
eq
⬅ 61 共 ␥ 1i ⫹6 ␥ 2i 兲 ,

can be evaluated by taking derivatives of the internal energy 共12兲

of the ith atomic site and evaluating the power moments of
the atomic displacements. If the ith atom in the lattice is respectively. In the SMA TB scheme, the parameters k i , ␥ 1i ,
affected by a supplementary force ␣ ␤ , then the total force and ␥ 2i are composed of two contributions 共band structure
acting on it must be zero, and one gets the force balance and repulsive energies兲 and k i is given by the following
relation as form:

兺␣ 冉 冊
⳵ 2 E ci
⳵ u i␣⳵ u i␤
具 u i ␣ 典 ⫹ 21 兺
eq ␣,␥
冉 ⳵ 3 E ci
⳵ u i␣⳵ u i␤⳵ u i␥ 冊 具 u i␣u i␥典
eq
k i⫽
q
冋 冉 冊 册
␩ 共 2 兲⫺ 2
r0 i
q 共 3 兲 ⫺1/2
␩ ␮ 2i ⫺A 共 p/r 0 兲
r0 i 兺j 冋 1⫺ᐉ 2i j
rij

兺 冉 冊 冉 冊册
1 ⳵ E ci4
⫹ 具 u i ␣ u i ␥ u i ␩ 典 ⫺ ␣ ␤ ⫽0. p
3! ␣ , ␥ , ␩ ⳵ u i ␣ ⳵ u i ␤ ⳵ u i ␥ ⳵ u i ␩ eq
⫺ᐉ 2i j exp关 ⫺ p 共 r i j /r 0 ⫺1 兲兴 , 共13兲
r0
共7兲
Here, the subscript eq indicates evaluation at equilibrium. where ␩ (2)
i and ␩ (3)
i are defined, respectively, as
The thermal averages of the atomic displacements 具 u i ␣ u i ␥ 典
and 具 u i ␣ u i ␥ u i ␩ 典 共called second- and third-order moments兲 at
given site Ri can be expressed in terms of the first moment
具 u i ␣ 典 with the aid of Eq. 共5兲 as
␩ 共i 2 兲 ⫽ 兺
j
冋 册 1⫺l 2i j
rij
␰ 2i j exp关 ⫺2q 共 r i j /r 0 ⫺1 兲兴 , 共14兲

⳵ 具 u i␣典 a
具 u i␣u i␥典 a⫽ 具 u i␣典 a具 u i␥典 a⫹ ␪ ␩ 共i 3 兲 ⫽ 兺 l 2i j ␰ 2i j expb ⫺2q 共 r i j /r 0 ⫺1 兲 c , 共15兲
⳵␣␥
j

⫹
ប ␦ a␥
2m ␻
coth
ប␻
2␪
⫺
␪ ␦ a␥
m␻2 冉 冊
, 共8兲
with

具 u i ␣ u i ␥ u i ␩ 典 a ⫽ 具 u i ␣ 典 a 具 u i ␥ 典 a 具 u i ␩ 典 a ⫹ ␪ P ␣␥ ␩ 具 u i ␣ 典 a
⳵ 具 u i␥典 a
⳵␣␩ li j⫽ 冉 冊
⳵rij
⳵x
⫽ 共 x j ⫺x i 兲 /r i j .

⫹␪2
⳵ 2 具 u i ␣ 典 a ប 具 u i ␩ 典 a ␦ ␣␥
⳵ ␣ ␥⳵ ␣ ␩
⫹
2m ␻
coth
ប␻
2␪ 冉 冊 After a bit of algebra, ␥ 1i defined by Eq. 共12兲 is given by

具 u i ␩ 典 a ␦ ␣␥
⫺␪
m␻2
.

Here, P ␣␥ ␩ is 1 ( ␣ ⫽ ␥ ⫽ ␩ ) or 0 共otherwise兲 depending on ␣,

共9兲
␥ 1i ⫽ 冉 冊冋
q
r0
⳵ 2 ␩ 共i 2 兲
⳵x2
⫺2
⳵ 2 ␩ 共i 3 兲 q
⳵x2 r0 冉 冊册 ⫺1/2
␮ 2i

␥, and ␩ 共Cartesian component兲 and ␻ is the atomic vibration

frequency similar to that defined in the Einstein model,
⫺ 冉 冊冋 q
r0
2
␩ 共i 2 兲 ⫺2 ␩ 共i 3 兲 冉 冊册
q
r0
2
⫺3/2
␮ 2i

冉 冊兺 冋
which will be given by Eq. 共11兲. Then Eq. 共7兲 is transformed
into the new differential equation p 3 共 1⫺6l 2i j ⫹5l 4i j 兲
⫹A
r0 j r 3i j

␥ i␪ 2
d2y
d␣2
⫹3 ␥ i ␪ y
dy
d␣
⫹ ␥ i y 3 ⫹k i y ⫹
3 共 1⫺6l 2i j ⫹5l 4i j 兲 p
r 2i j r0
⫺ 冉冊
6l 2i j p
rij r0 冉冊 2

冉 冊册
␪ 4
⫹ ␥ i 共 X coth X⫺1 兲 y⫺ ␣ ␤ ⫽0, 共10兲 p
k ⫹l 4i j exp兵 ⫺ p 共 r i j /r 0 兲 ⫺1 其 . 共16兲
r0

where X⬅ប ␻ /2␪ and y⬅ 具 u i 典 . Here, k i and ␥ i are second-

and fourth-order derivatives of E ci and defined by the fol- The second derivatives of ␩ (2)i and ␩ (3)
i appearing in the first
lowing formulas: term of the right-hand side of Eq. 共16兲 are also given explic-
itly in terms of the TB potential parameters and the direction

k i⫽ 冉 冊
⳵ 2 E ci
⳵ u i2␣ eq
⬅m␻ 2 , 共11兲
cosines l i j and m i j between the central atom i and its neigh-
boring atoms j 共see Appendix B兲. ␥ 2i is expressed explicitly
as

094301-4
THERMODYNAMIC QUANTITIES OF METALS . . . PHYSICAL REVIEW B 67, 094301 共2003兲

TABLE II. Lattice sums appearing in the harmonic k 1 and anharmonic ␥ 1 and ␥ 2 parameters in cubic
metals. 兺 1 ⬅ 兺 j⫽i 1⫺6l i2j ⫹5l i4j , 兺 2 ⬅ 兺 j⫽i 1⫺3l i2j ⫺3m i2j ⫹15l i2j m i2j , 兺 3 ⬅ 兺 j⫽i l i2j ⫹m i2j ⫺6l i2j m i2j .

Crystal structure Neighbors 1 2 3 4 5

fcc Zi 12 6 24 12 24
Distance 1 & ) 2 冑5

兺l
j⫽i
2
ij 4 2 8 4 8

兺l
j⫽i
4
ji 2 2 4 2 164/25

兺l m
j⫽i
2
ij
2
ij 1 0 2 1 18/25

兺 1
⫺2 4 ⫺4 ⫺2 44/5

兺 2
3 ⫺6 6 3 ⫺66/5

兺 3
2 4 4 2 292/25

bcc Zi 8 6 12 24 8
Distance 1 2/) 2 冑6/3 冑11/3 2

兺l
j⫽i
2
ij 8/3 2 4 8 8/3

兺l
j⫽i
4
ji 8/9 2 2 664/121 8/9

兺l m
j⫽i
2
ij
2
ij 8/9 0 1 152/121 8/9

兺 1
⫺32/9 4 ⫺2 416/121 ⫺32/9

兺 2
16/3 ⫺6 3 ⫺624/121 16/3

兺 3
0 4 2 1024/121 0

冉 冊冋 ⳵ 2 ␩ 共i 2 兲 ⳵ 2 ␩ 共i 3 兲 q
冉 冊册 冋 册冉 冊
⳵ ␩ 共i 1 兲
冉 冊冋 冉 冊册
2 2 2 2
q ⫺1/2
q ⫺3/2
q q ⫺5/2
␥ 2i ⫽ ⫺2 ␮ 2i ⫺2 ␮ 2i ⫹ ␩ 共i 2 兲 ⫺2 ␩ 共i 3 兲 ␮ 2i
r0 ⳵y2 ⳵y2 r0 ⳵y r0 r0 r0

⫹A冉 冊兺 冋 p
r0 j
1⫺3l 2i j ⫺3m 2i j ⫹15l 2i j m 2i j
r 3i j
⫹
1⫺3l 2i j ⫺3m 2i j ⫹15l 2i j m 2i j p
r 2i j r0 冉冊
⫺
l 2i j ⫹m 2i j ⫺6l 2i j m 2i j
rij 冉冊
p
r0
2
⫹l 2i j m 2i j 冉 冊册p
r0
3
exp兵 ⫺ p 共 r i j /r 0 ⫺1 兲 其 , 共17兲

where ␩ (1) is defined by In determining the atomic displacement 具 u i 典 , the symme-

i
try property appropriate for cubic crystals is used

␩ 共i 1 兲 ⫽ 兺 l i j ␰ 2i j exp关 ⫺2q 共 r i j /r 0 ⫺1 兲兴 . 共18兲

具 u i␣典 ⫽ 具 u i␥典 ⫽ 具 u i␩典 ⬅ 具 u i典 . 共19兲
j
Then, the solutions of the nonlinear differential equation of
Eq. 共10兲 can be expanded in the power series of the supple-
Here, we note that ␥ 1i and ␥ 2i depend sensitively on the mental force ␣ as
structure of crystals through factors including direction co-
sines as can be seen in Eqs. 共16兲 and 共17兲. The factors in- y⫽⌬r⫹A 1 ␣ ⫹A 2 ␣ 2 . 共20兲
cluding direction cosines for cubic crystals are presented in
Table II. The derivatives of ␩ (1)
i , ␩ i , and ␩ i
(2) (3)
with respect Here, ⌬r is the atomic displacement in the limit of zero of
to the y variable are given in Appendix B. supplemental force ␣. Substituting the above solution of Eq.

094301-5
K. MASUDA-JINDO, VU VAN HUNG, AND PHAM DINH TAM PHYSICAL REVIEW B 67, 094301 共2003兲

共20兲 into the original differential equation Eq. 共10兲, one can where the second term denotes the harmonic contribution to
get the coupled equations on the coefficients A 1 and A 2 , the free energy.
from which the solution of ⌬r is given as With the aid of the free energy formula ⌿⫽E⫺TS, one
can find the thermodynamic quantities of metal systems. The
共 ⌬r 兲 2 ⬇ 关 ⫺C 2 ⫹ 冑C 22 ⫺4C 1 C 3 兴 /2C 1 , 共21兲 specific heats and elastic moduli at temperature T are directly
where derived from the free energy ⌿ of the system. For instance,
the isothermal compressibility ␹ T is given by
C 1 ⫽3 ␥ i ,

C 2 ⫽3k i 1⫹ 冋 ␥ i␪
k 2i
共 X coth X⫹1 兲 , 册 共22兲 ␹ T ⫽3 共 a/a 0 兲 3 冒冋 2 P⫹
1 & ⳵ 2⌿
3N a ⳵ r 2冉 冊册 T
, 共28兲

C 3 ⫽⫺
2 ␥ i␪ 2
k 2i
冉 1⫹
X coth X
2
. 冊 where

Using Eqs. 共8兲 and 共21兲, it can be shown that mean square
atomic displacement 共second moment兲 in cubic crystals is
given by
⳵ 2⌿
⳵r2
⫽3N
1 ⳵ 2U 0
6 ⳵r2 再⫹ ␪
2k ⳵r2冋
X coth X ⳵ 2 k

␪
具 u 2 典 ⫽ X coth X⫹
k
2 2 ␥␪ 2
3 k3
共 1⫹X coth X/2兲 ⫺
1 ⳵k
冉 冊冉
4k 2 ⳵ r
2
X coth X⫹
X2
sinh2 X 冊 册冎 . 共29兲

2 ␥ 2␪ 3
⫹ 共 1⫹X coth X 兲共 1⫹X coth X/2兲 . 共23兲
k5 On the other hand, the specific heats at constant volume C v
is
Once the thermal expansion ⌬r in the lattice is found, one
can get the Helmholtz free energy of the system in the fol-
lowing form:

⌿⫽U 0 ⫹⌿ 0 ⫹⌿ 1 , 共24兲
C v ⫽3Nk B 再 X2
2
2␪
⫹ 2
sinh X k 冋冉 2 ␥ 2⫹ 冊
␥ 1 X 3 coth X
3 sinh2 X

where ⌿ 0 denotes the free energy in the harmonic approxi-

mation and ⌿ 1 the anharmonicity contribution to the free
⫹
␥1
3 冉
1⫹
X2
sinh2 X
⫺ ␥ 2 冊 冉
X4
sinh4 X
⫹
2X 4 coth2 X
sinh2 X 冊 册冎 .
energy.38 – 40 We calculate the anharmonicity contribution to
the free energy ⌿ 1 by applying the general formula 共30兲

⌿⫽U 0 ⫹⌿ 0 ⫹ 冕具 ␭

0
V̂ 典 ␭ d␭, 共25兲 The specific heat at constant pressure C p is determined from
the thermodynamic relations
where ␭V̂ represents the Hamiltonian corresponding to the
anharmonicity contribution. It is straightforward to evaluate
9TV ␣ T2
the following integrals analytically C p ⫽C v ⫹ , 共31兲
␹T

I 1⫽ 冕0
␥1
具 u 4i 典 d ␥ 1 , I 2⫽ 冕 0
␥2
具 u 2i 典 ␥2 1 ⫽0 d ␥ 2 . 共26兲
where ␣ T denotes the linear thermal expansion coefficient
Then the free energy of the system is given by and ␹ T the isothermal compressibility. In Eqs. 共27兲, 共29兲, and
共30兲 above, the suffices i for the parameters k, ␥ 1 and ␥ 2 are
⌿⫽U 0 ⫹3N ␪ 关 X⫹ln共 1⫺e ⫺2X 兲兴 omitted because each atomic site is equivalent in a mono-

再 冋 冉 冊册
atomic cubic crystal with primitive structure. The relation-
␪2 2 X coth X
⫹3N ␥ X 2 coth2 X⫺ ␥ 1 1⫹ , ship between the isothermal and adiabatic compressibilities,
k2 2 3 2 ␹ T and ␹ s , is simply given by

⫹ 4
k 3 冋
2␪3 4 2
␥ 2 X coth X 1⫹
X coth X
2 冉 冊 Cv

冉 冊 册冎
␹ s⫽ ␹ . 共32兲
X coth X Cp T
⫺2 ␥ 1 共 ␥ 1 ⫹2 ␥ 2 兲 1⫹ 共 1⫹X coth X 兲 ,
2
共27兲 One can also find ‘‘thermodynamic’’ Grüneisen constant as

094301-6
THERMODYNAMIC QUANTITIES OF METALS . . . PHYSICAL REVIEW B 67, 094301 共2003兲

FIG. 1. Comparison of linear thermal expansion coefficients ␣ T of 共a兲 Cu, 共b兲 Pd, 共c兲 Ag, and 共d兲 Mo, calculated by using the Morse
potentials. Solid and dot-dashed lines show the results of self-consistent 共SC兲 and non-self-consistent 共NSC兲 treatments of the statistical
moment method, respectively, while the dashed ones are the results of the QH theory by Moruzzi, Janak, and Schwarz 共MJS兲.

␥ G⫽ 冋 册
V ⳵S
C␷ ⳵V T
⫽
␣ TB SV
CP
, 共33兲
energy F(T,V) on the system volume V can be explored by
homogeneous scaling of the atomic potentials 兵 R 0i 其 . Then,
for each temperature T the equilibrium volume V is obtained
where B S ⬅ ␹ ⫺1 by minimizing Helmholtz energy F with respect to V. In Fig.
S denotes the adiabatic bulk modulus.
1, we present the linear thermal expansion coefficients ␣ T of
Cu, Pd, Ag, and Mo metals, calculated by the present theory,
III. RESULTS AND DISCUSSIONS together with those of the QH theory by Moruzzi, Janak, and
A. Comparison with the quasiharmonic theory
Schwarz 共MJS model兲.10 The linear thermal expansion coef-
ficients ␣ T by the present statistical moment theory and those
Firstly, we compare the thermodynamic quantities of met- of the QH theory by Moruzzi et al. are referred to as SMM
als calculated by the present statistical moment method and MJS, respectively. In order to allow the direct compari-
共SMM兲 with those by the QH theory.10 The basic idea of the son between the two different schemes, the linear thermal
QH approximation is that the explicit dependence of the free expansion coefficients ␣ T of the cubic metals are calculated

094301-7
K. MASUDA-JINDO, VU VAN HUNG, AND PHAM DINH TAM PHYSICAL REVIEW B 67, 094301 共2003兲

with the use of the same Morse type of potentials, exactly

identical forms as used in the QH calculations by MJS.10 The
four metals Cu, Pd, Ag, and Mo are chosen simply because
the linear thermal expansion coefficients ␣ T are well repro-
duced by the two-body Morse potentials as demonstrated by
them.10
The solid lines in Fig. 1 show the linear thermal expan-
sion coefficients ␣ T calculated by the self-consistent 共SC兲
treatments of the present SMM scheme, while the dot-dashed
ones are obtained by the non-self-consistent 共NSC兲 treat-
ments. In the SC treatments, the characteristic parameters k,
␥ 1 , and ␥ 2 are determined self-consistently with the lattice
constants a T at given temperature T. However, in the NSC
treatments, the harmonic k, and anharmonic ␥ 1 , and ␥ 2 pa-
rameters are fixed to those values evaluated at the appropri-
ate reference temperature T 0 共e.g., absolute zero temperature
or some reference temperature; here T 0 is chosen to be 0 K
and taken to be constant for the whole temperature region兲.
The calculated linear thermal expansion coefficients ␣ T by
the present theory are in good agreement with those by QH
theory for the lower temperature region below the Debye
temperature and the agreement is better for the SC calcula-
tions. This indicates that the thermal lattice expansion gives
rise to the significant reduction in the parameters k, ␥ 1 , and
␥ 2 , and thereby changes the thermodynamic quantities ap-
preciably even for the lower temperatures.

B. Thermodynamic quantities of metals by second moment

TB potentials
With the use of the analytic expressions presented in Sec.
II, it is straightforward to calculate the thermodynamic quan-
tities of metals and alloys at the thermal equilibrium. Firstly,
the equilibrium lattice spacings are determined, using Eqs.
共20兲 and 共21兲, in the SC treatment including temperature
共bond length兲 -dependent k, ␥ 1 , and ␥ 2 values. The thermal
lattice expansion can also be calculated by standard proce-
dure of minimizing the Helmholtz energy of the system: We
have checked that both calculations give almost identical re-
sults on the thermal lattice expansions. We calculate the ther- FIG. 2. 共a兲 The linear thermal expansion coefficient ␣ T 共a兲 and
mal lattice expansion and mean-square atomic displacements 共b兲 mean-square atomic displacements 具 u 2 典 of Cu crystal calculated
of some fcc 共transition兲 metals and bcc alkali metals, for by the present method. The corresponding experimental values are
which the reliable many-body potentials are available, and presented by symbols 䊊.
compare them with those by the molecular dynamics 共MD兲
and Monte Carlo 共MC兲 simulations. So far, a number of the
SMA base TB potentials have been proposed for fcc metals. differences in the calculated quantities when we use the
Specifically, we use the SMA TB potentials by Rosato Lennard-Jones 共LJ兲 type of pair potentials. The bold line in
et al.35 and by Cleri and Rosato36 for fcc metals, which are Fig. 2共a兲 represents the calculated ␣ T by the present SMM,
known to give good descriptions of cohesive properties of while the dashed line ␣ T values by the Lennard-Jones type of
fcc elements. For alkali metals Li, Na, K, Rb, and Cs, we use potential; ␸ (r)⫽D 0 兵 (r 0 /r) n ⫺(n/m)(r 0 /r) m 其 , with n
the potential parameters proposed recently by Li et al.40 ⫽9.0, m⫽5.5, r 0 ⫽2.5487 Å, and D 0 ⫽4125.7 K 共0.35553
In the TB scheme by Rosato et al.,35 the interaction range eV兲, respectively. The overall agreement between the calcu-
is limited to the first nearest neighbors, while in the TB lated and experimental ␣ T values is better for the calcula-
scheme by Cleri and Rosato,36 it is extended to the fifth tions by the SMA TB potential, although LJ potential param-
neighbors. In Fig. 2, we present the linear thermal expansion eters are not best fitted to reproduce the experimental ␣ T
coefficients ␣ T and mean-square atomic displacements 具 u 2 典 values. We note that the classical MD simulation,59 shown by
of Cu crystal, together with the experimental values 共by sym- the dot-dashed curve in Fig. 2共a兲, do not reproduce the cor-
bols 䊊兲.55–58 For this calculation, the electronic many-body rect curvature of the linear thermal expansion coefficient ␣ T ,
potentials are used for Cu crystal, but there are no large and is qualitatively incorrect due to the absence of the quan-

094301-8
THERMODYNAMIC QUANTITIES OF METALS . . . PHYSICAL REVIEW B 67, 094301 共2003兲

differ significantly from those results by MD simulations,

especially for the lower temperature region, i.e., below the
Debye temperature. This is due to the fact that in the classi-
cal MD simulations the quantum mechanical vibration ef-
fects are not taken into account. One sees that the quantum
mechanical zero point vibrations give main contributions at
lower temperature region T⭐100 K. The agreement between
the present calculation and the experimental results is fairly
good for the whole temperature region, from zero to ⬃800
K, much higher than the Debye temperature. In Fig. 3共b兲, we
show the mean-square atomic displacements 具 u 2 典 of Ag
crystal calculated by the present SMM using the SMA TB
potentials of Refs. 35 and 36, together with the experimental
results.62 One sees in Fig. 3共b兲 that TB parameters by Rosato,
Guillope, and Legrand35 共first-neighbor TB potential兲 leads
to larger mean-square atomic displacement 具 u 2 典 in Ag crys-
tal compared to those results by using the TB parameters by
Cleri and Rosato36 共5th neighbor TB potential兲. The similar
tendency is also found for the thermal expansion coefficients
␣ T of Ag crystal, larger ␣ T values by TB potential by Rosato,
Guillope, and Legrand.35 In the present formalism, the ther-
mal lattice expansion and mean-square atomic displacements
are characterized by the harmonic k and anharmonic ␥ pa-
rameters. In particular, the thermal lattice expansion 共mate-
rial dependence兲 is predicted by a ratio of ␥ /k 2 and the
mean-square displacement 具 u 2 典 by ␥ /k 2 共and also by ␥ 2 /k 5 )
parameter as well. The ratios ␥ /k 2 of Cu crystal calculated
by using the TB potential by Rosato, Guillope, and
Legrand35 are in fact larger than those results by Cleri and
Rosato36 for whole temperature region. The mean square
atomic displacement 具 u 2 典 in Ag crystal by the fifth-neighbor
TB potential36 are in fairly good agreement with the experi-
mental results for the whole temperature region, and they are
in good agreement with the MD simulation results for high
temperature region.
The calculated mean-square atomic displacements 具 u 2 典 of
Ag crystal by the present method is also compared with
those by the cluster variation method 共CVM兲. As is well
known, CVM63– 65 is an analytical statistical method that di-
rectly gives us the free energy of a system. The CVM was
FIG. 3. Mean-square atomic displacements 具 u 2 典 of 共a兲 Al and originally designed for the statistical mechanics of the Ising
共b兲 Ag crystals as a function of temperature. In 共a兲, the dashed line model on a fixed lattice, and extended recently to treat sys-
shows the results of MD simulations by Papanicolaou et al. 共Ref. tems with continuous degrees of freedom, such as the lattice
60兲, while the solid circles are the experimental values. site distortion, due to thermal vibrations, thermal dilatation,
and mixture of atoms of different sizes. In general, in CVM
treatments the correlations in the atomic displacements are
tum mechanical vibration effects. One also sees in Fig. 2共b兲 taken into account within the small atomic clusters 共e.g.,
that the agreements between the calculated and experimental small clusters such as pair, tetrahedron, or octahedron clus-
results of the mean square atomic displacements 具 u 2 典 in Cu ters兲. Finel and Tétot gave the first application of the Gauss-
crystal are quite excellent for the SMA TB calculations, com- ian CVM65 for the thermodynamic quantities of some transi-
pared to those by two-body potentials. This implies that the tion metals. It has been demonstrated that Gaussian CVM
present SMM scheme with SMA TB potentials provides us gives the excellent results of the thermodynamic quantities
fully quantitative estimates for the thermodynamic quantities of metals 共the CPU time is several orders of magnitude
of elemental metals. smaller than the one needed for numerical MD or MC simu-
We show in Fig. 3共a兲 the mean-square atomic displace- lations兲. The thin dot-dashed and thin dashed curves in Fig.
ments 具 u 2 典 of Al crystal as a function of temperature T, 3共b兲 represent the mean-square atomic displacement 具 u 2 典 of
together with those values by the MD simulation60 and ex- Ag crystal obtained by the Gaussian CVM65 using the SMA
perimental data.61 The present calculations by using SMM TB potentials of Refs. 35 and 36, respectively. Both CVM

094301-9
K. MASUDA-JINDO, VU VAN HUNG, AND PHAM DINH TAM PHYSICAL REVIEW B 67, 094301 共2003兲

TABLE III. Bulk modulus, linear thermal expansion, and Grüneisen constant calculated with the use of
the SMA TB potentials. Experimental values of Na*共RT兲 are those values for 250 K.

B T (GPa) ␣ (10⫺6 K⫺1 ) ␥G

Calc.

Element T⫽0 RT Expt. Calc. Expt. Calc. Expt.

Al 87 75 72 24.5 23.6 2.09 2.19

Cu 153 137 137 15.9 16.7 2.21 2.00
Ni 190 182 184 14.7 12.7 2.01 1.88
Ag 114 96 101 23.5 19.7 2.78 2.36
Rh 306 280 271 10.9 8.2 2.19 2.43
Pd 204 171 181 14.3 11.6 2.22 2.18
Au 185 164 173 17.2 14.2 3.21 3.04
Pt 301 259 278 11.2 8.9 3.06 2.56
Li 16.8 12.4 11.6 65.4 56.0 1.18 1.18
Na* 6.5 4.3 6.8 83.9 71.0 1.53 1.31
K 5.3 3.6 3.2 98.7 83.0 1.54 1.37
Rb 4.0 2.8 3.1 104.6 90.0 1.65 1.67
Cs 2.9 2.1 2.0 108.8 97.0 1.55 1.44

calculations of ␣ T are generally in agreement with the ex- lated Grüneisen constants ␥ G for low temperatures are well
perimental results. We note that for 具 u 2 典 calculations of Ag compared with the experimental values which are deduced
crystal, however, the present analytic SMM gives much effi- from the low 共room-兲 temperature specific heats.
cient analytic calculations and much better results compared The lattice specific heats C v and C p at constant volume
to those by CVM calculations. and at constant pressure are calculated using Eqs. 共30兲 and
The calculated thermodynamic quantities of cubic metals, 共31兲, respectively. However, the evaluations by Eqs. 共30兲 and
fcc 共in addition to Cu, Ag, and Al presented above兲 and alkali 共31兲 are the lattice contributions, and their values may not be
共bcc兲 metals, by the present method are summarized in Table directly compared with the corresponding experimental val-
III. In the present calculations, we use the TB potential pa- ues. We do not include the contributions of lattice vacancies
rameters by Li, Barojas, and Papaconstantopoulos40 for al- and electronic parts of the specific heats C v , which are
kali metals Li, Na, K, Rb, and Cs. This TB model takes into
known to give significant contributions in metals for higher
account the interatomic interactions up to 12th neighbors,
temperature region near the melting temperature. In particu-
i.e., 228 atoms in bcc lattice. The relative magnitudes of
lar, it has been demonstrated that lattice vacancies make a
linear thermal expansion coefficients of fcc 共transition兲 met-
als are in good agreement with the experimental results. large contribution to the specific heats for the high-
However, the thermal lattice expansion coefficients ␣ of al- temperature region.66 The electronic contribution to the spe-
kali metals are systematically larger 共⬃10%兲 than those of cific heat at constant volume C ele
v is proportional to the tem-
experimental results, although their relative magnitudes are perature T and given by C elev ⫽ ␥ e T, ␥ e being the electronic
in good agreement with the experimental results. The calcu- specific heat constant.56,66 The electronic specific heats C ele
v
lated Grüneisen constants and elastic moduli are also pre- values are estimated to be 0.8 –13.4% of C latv for metals con-
sented in Table III. The anharmonicity of the lattice vibra- sidered here by the free-electron model.56 Therefore, the
tions is well described by the Grüneisen constant ␥ G . The present formulas of the lattice contribution to the specific
material of larger value of ␥ G may be regarded as the mate- heats, both C v and C p , for the cubic metals tend to under-
rial with higher lattice anharmonicity. So, the evaluation of estimate the specific heats for higher temperature region,
the Grüneisen constant is of great significance for the assess- when compared with the experimental results. The lattice
ment of anharmonic thermodynamic properties of metals and contribution of specific heats C p calculated for Cu crystal is
alloys. The experimental Grüneisen constants ␥ G of fcc met- shown in Fig. 4, together with the experimental results58 and
als are larger than 2 except for Ni, while those of alkali those of MD simulation results. As expected from above
metals are less than 2 and take values around ⬃1.5. The mentioned reasonings, the calculated C p values of solid Cu
calculated Grüneisen constants ␥ G of fcc metals are also are smaller than the experimental values at high tempera-
larger than 2, while those values of alkali metals are less than tures. However, the temperature dependence 共curvature兲 of
2, in agreement with the experimental results. The calculated C p of Cu crystal by the present method is in good agreement
␥ G values by the present method have the weak temperature with the experimental results, in contrast to the MD simula-
dependence, i.e., show the slight increase with increasing tion results. In the MD simulations, the heat capacities per
temperature as in the calculations by QH theory.10 The tabu- atom at constant pressure C p can be obtained for metals by

094301-10
THERMODYNAMIC QUANTITIES OF METALS . . . PHYSICAL REVIEW B 67, 094301 共2003兲

the evaluations of k, ␥ 1 , and ␥ 2 , and thus for the calcula-

tions of thermodynamic quantities of the present study.

IV. CONCLUSIONS

We have presented an analytic formulation for obtaining

the thermodynamic quantities of metals and alloys based on
the finite temperature moment expansion technique in the
statistical physics. The thermal lattice expansion of mon-
atomic crystals 共fourth-order anharmonic contribution兲 is de-
rived explicitly in terms of the three characteristic param-
eters, k 1 , ␥ 1 , and ␥ 2 . The present formalism takes into
account the quantum-mechanical zero-point vibrations as
well as the higher-order anharmonic terms in the atomic dis-
placements and it enables us to derive the various thermody-
namic quantities of metals and alloys for a wide temperature
range. We are able to calculate the thermodynamic quantities
quite efficiently and accurately by using the analytic formu-
las and taking into account the many-body electronic effects
in metallic systems. The calculated thermodynamic quanti-
ties of metals are in good agreement with the experimental
FIG. 4. Specific heats per atom at constant pressure C p plotted results as well as with those by MD and MC simulations 共in
against temperature T for Cu crystal, in unit of Boltzmann’s con- some cases, better results by the present method兲.
stant k B . The experimental data 共Ref. 58兲 are shown as the open Although in this paper we only used many-body elec-
circles. tronic potentials, the extension to coupling the present SMM
scheme with the ab initio density functional theories is
straightforward. This can be done by evaluating three char-
taking the numerical derivative of the internal energy with acteristic parameters k, ␥ 1 , and ␥ 2 for cubic systems. It can
respect to temperature.59 The MD simulations by Mei, Dav- also be applied directly for the composition-temperature
phase diagrams calculations of alloys for the full temperature
enport, and Fernando59 give reasonable values of C p for
range from absolute zero to the melting temperatures T m .
higher temperature region when compared with the experi-
mental data. However, it should be noted that the MD simu-
lations are only adequate above the room temperature, and ACKNOWLEDGMENTS
the calculated C p value deviates from the experimental data
at low temperatures because quantum effects are not taken The authors would like to thank Professor S. Tsuneyuki of
Institute for Solid State Physics of the University of Tokyo
into account in the classical MD simulations.
for valuable discussions. The support of the supercomputing
The bulk moduli B T of cubic metals are evaluated at ab-
facility of Institute for Solid State Physics, the University of
solute zero temperature and at room temperature 共RT兲 and
Tokyo is also acknowledged.
presented also in Table III. The ratios of bulk moduli,
B T /B 0 , with respect to those of the absolute zero tempera-
ture are calculated to be 0.85–0.90 共fcc metals兲 and ⬃0.7 APPENDIX A: MOMENT DEVELOPMENT BY DENSITY
共alkali metals兲 which are favorably compared with the ex- MATRIX FORMALISM
perimental results. In general, the calculated bulk moduli
To derive the mean-square atomic displacement 共second
B T (RT) and B 0 are in good agreement with the experimental
moment兲 and higher-order power moments of the thermal
results as well as QH calculations.
lattice vibrations, we use the formalism based on the density
As a final remark of this section, we note that the present
matrix ␳ˆ , which is defined by
statistical moment method can be incorporated in a straight-

冉 冊
forward manner with the first-principles density functional
theory, by simply evaluating three kinds of derivatives 共one ⌿⫺Ĥ
for harmonic and two for anharmonic contributions兲 of the ␳ˆ ⫽exp , 共A1兲
␪
atomic total energy with respect to the Cartesian coordinates.
The density functional TB67 and TBTE68 共tight-binding total
energy兲 methods with Slater-Koster parameters derived from where ⌿ and Ĥ denote the Helmholz free energy and Hamil-
the first-principles theories can be readily applied to evaluate tonian of the system, respectively. In the presence of the
k, ␥ 1 , and ␥ 2 values, on the basis of Hellman-Feynman theo- constant supplemental forces ␣ 1 , ␣ 2 ,..., ␣ N in the system,
rem. The full density functional theories such as the linear- the Hamiltonian Ĥ is given by Ĥ⫽Ĥ 0 ⫺ 兺 i ␣ i q̂ i . The density
response approach by Giannozzi et al.69 and real-space matrix ␳ˆ is normalized so as to satisfy the condition Tr␳ˆ
finite-element density matrix method70 can also be used for ⫽1 and given by the solution of the Liouville equation

094301-11
K. MASUDA-JINDO, VU VAN HUNG, AND PHAM DINH TAM PHYSICAL REVIEW B 67, 094301 共2003兲

⳵␳ˆ 具 共 q̂⫺ 具 q̂ 典 兲共 F̂⫺ 具 F̂ 典 兲 典

iប ⫽ 关 Ĥ, ␳ˆ 兴 ⫺ . 共A2兲

冉 冓 冔冊
⳵t
⳵ 具 F̂ 典 ⳵ F̂
⫽⫺ ␪ ⫺
We use the following identities on the derivatives of an op- ⳵␣ ⳵␣

冉 冊冓 冔
erator function Ê ␭ composed of two different operators  ⬁
Bm iប m
⳵ F̂ 共 m 兲
and B̂: ⫹␪ 兺
m⫽1
共 ⫺1 兲 m
m!
⫺
␪ ⳵␣
. 共A9兲

Ê ␭ 共 ␶ ,Â,B̂ 兲 ⬅exp关 ␶ 共 Â⫹␭B̂ 兲兴 共A3兲 Then, one gets the decoupling formula

⳵ Ê ␭ 共 ␶ ,Â,B̂ 兲
⳵␭ 再
⫽ ␶ B̂⫹
⬁
␶ n⫹1 共 iប 兲 n
兺
n⫽1 共 n⫹1 兲 !
关 Â⫹␭B̂ 关 Â⫹␭B̂... 关 Â
1
2
具 关 F̂,q̂ k 兴 ⫹ 典 a ⫺ 具 F̂ 典 a 具 q̂ k 典 a

冎
⫹␭B̂... 兴兴兴 Ê ␭ 共 ␶ ,Â,B̂ 兲 ⫽␪
⳵ 具 F̂ 典 a
⳵␣k
⫺␪ 兺
⬁
B 2n iប
n⫽0 共 2n 兲 ! ␪
冉 冊冓 冔 2n
⳵ F̂ 共 2n 兲
⳵␣k
,

再
a
⬁
共 ⫺ ␶ 兲 n⫹1 共 iប 兲 n 共A10兲
⫽Ê ␭ 共 ␶ ,Â,B̂ 兲 ␶ B̂⫺ 兺
n⫽1 共 n⫹1 兲 !
In the above Eqs. 共A8兲–共A10兲, B 2n denotes the Bernoulli

⫻ 关 Â⫹␭B̂ 关 Â⫹␭B̂... 关 Â⫹␭B̂... 兴兴兴 , 冎 共A4兲

number and F̂ (k) is defined by

1
where 关 Â⫹␭B̂,B̂ 兴 ⫽ 兵 共 Â⫹␭B̂ 兲 B̂⫺B̂ 共 Â⫹␭B̂ 兲 其 . 共A11兲
iប

By differentiation of the density matrix ␳ˆ with respect to the Substituting F̂⫽q̂ k into Eq. 共A10兲, one can get the mean-
constant force ␣ i , one can get the relation square atomic displacement from the thermal equilibrium po-
sition, as
1 ⳵⌿ 1
冋 1
兺
⬁
iប
冉 冊 n
册
具 q̂ 共i n 兲 典 a ⫽0,
冉 冊冓 冔
⫹ ⫺ 具 q̂ i 典 a ⫺
␪ ⳵␣i ␪ n⫽1 共 n⫹1 兲 ! ␪ ⳵ 具 q̂ i 典 a
⬁
B 2n iប 2n
⳵ q̂ 共i 2n 兲
共A5兲 具 共 q̂ i ⫺ 具 q̂ i 典 a 兲 2 典 a ⫽ ␪
⳵␣i
⫺␪ 兺
n⫽0 共 2n 兲 ! ␪ ⳵␣i a
.

where 共A12兲

Equation 共A12兲 is used to derive Eq. 共8兲 in the text. The

具 q̂ 共i n 兲 典 ⫽Tr关关 ... 关 q̂ i ,Ĥ 兴 ⫺ Ĥ... 兴 ␳ˆ . 共A6兲 similar formulas can be given for higher order moments as
well.
For equilibrium state, one can show that 具 q̂ (n) i 典 ⫽0 because
⳵␳ˆ / ⳵ t⫽ 关 Ĥ, ␳ˆ 兴 ⫺ ⫽0 is satisfied. One can then derive from APPENDIX B: DERIVATIVES OF COUPLING
共A5兲 the relation PARAMETERS k AND ␥

The second derivatives such as ⳵ 2 ␩ (2) i /⳵x

2
and
⳵⌿
具 q̂ k 典 a ⫽ . 共A7兲 ⳵ ␩ i / ⳵ x , appearing in Eq. 共16兲 in the text are given by the
2 (3) 2
⳵␣k following forms, respectively:

Using the relation 具 F̂ 典 a ⫽Tr关 F̂ ␳ˆ 兴 and Eqs. 共A3兲 and 共A4兲,

one can get the following identities: ⳵ 2 ␩ 共i 2 兲
⳵x2
⫽ 兺j 冋 ⫺3 共 1⫺6l 2i j ⫹5l 4i j 兲 r ⫺3
ij

具 共 F̂⫺ 具 F̂ 典 兲共 q̂⫺ 具 q̂ 典 兲 典 ⫽⫺ ␪ 冉 ⳵ 具 F̂ 典
⳵␣
⫺ 冓 冔冊
⳵ F̂
⳵␣ ⫺2 共 1⫺8l 2i j ⫹7l 4i j 兲 r -2
ij 冉冊 q
r0

冉 冊冓 冔 冉 冊册
⬁
B m iប m
⳵ F̂ 共 m 兲 q 2
⫹␪ 兺
m⫽1
共 ⫺1 兲 m
m! ␪ ⳵␣
, ⫹4l 2i j 共 1⫺l 2i j 兲 r ⫺1
ij
r0
␰ 2i j exp关 ⫺2q 共 r i j /r 0 ⫺1 兲兴 ,

共A8兲 共B1兲

094301-12
THERMODYNAMIC QUANTITIES OF METALS . . . PHYSICAL REVIEW B 67, 094301 共2003兲

⳵ 2 ␩ 共i 3 兲
⳵x2
⫽ 兺j 冋 2 共 1⫺5l 2i j ⫹4l 4i j 兲 r ⫺2 ⫺1
i j ⫹10l i j 共 1⫺l i j 兲 r i j
2 2
冉冊
q
r0
⳵ 2 ␩ 共i 2 兲
⳵y2
⫽⫺ 兺j 冋 共 1⫺3l 2i j ⫺3m 2i j ⫹15l 2i j m 2i j 兲 r ⫺3
ij

⫹4l 4i j 冉 冊册
q
r0
2
␰ 2i j exp关 ⫺2q 共 r i j /r 0 ⫺1 兲兴 . 共B2兲
⫹2 共 1⫺3m 2i j ⫺l 2i j ⫹7l 2i j m 2i j 兲 r ⫺1
ij 冉冊
q
r0

⫺4 共 m 2 ⫺l 2 m 2 兲 r ⫺1
ij 冉 冊册
q
r0
2
␰ 2i j

On the other hand, the first and second derivatives such as ⫻exp关 ⫺2q 共 r i j /r 0 ⫺1 兲兴 , 共B4兲
⳵ ␩ (1)
i / ⳵ y, ⳵ ␩ i / ⳵ y , and ⳵ ␩ i / ⳵ y with respect to the y
2 (2) 2 2 (3) 2
and

冋
variable are given by
⳵ 2 ␩ 共i 3 兲
⳵y2
⫽ 兺j ⫺2 共 l 2i j ⫺4l 2i j m 2i j 兲 r ⫺2
i j ⫺2l i j 共 1
2

⳵ 2 ␩ 共i 1 兲
⳵y
⫽ 兺j 冋 l i j m i j •r ⫺1
i j ⫹2l i j m i j 冉 冊册
q
r0
␰ 2i j ⫺5m 22 兲 r ⫺1
ij ij 冉冊q
r0
⫹4l 2i j m 2i j
q
r0 冉 冊册 2
␰ 2i j

⫻exp关 ⫺2q 共 r i j /r 0 ⫺1 兲兴 , 共B5兲

⫻exp关 ⫺2q 共 r i j /r 0 ⫺1 兲兴 , 共B3兲 respectively.

1 19
P. E. A. Turchi and A. Gonis, Statics and Dynamics of Alloy G. Leibfried and W. Ludwig, Theory of Anharmonic Effects in
Phase Transformations 共Plenum, New York, 1994兲, p. 345. Crystals 共Academic, New York, 1961兲.
2
Foundations and Applications of Cluster Variation Method and 20
J. Rosen and G. Grimvall, Phys. Rev. B 27, 7199 共1983兲.
21
Path Probability Method, edited by T. Morita, M. Suzuki, K. Van de Walle, G. Ceder, and U. V. Waghmare, Phys. Rev. Lett. 80,
Wada, and M. Kaburagi, Suppl. Prog. Theor. Phys. 115, 1 4911 共1998兲.
共1994兲. 22
Van de Walle and G. Ceder, Phys. Rev. B 61, 5972 共2000兲.
3 23
J. L. Morán-López and J. M. Sanchez, Theory and Applications of T. Yokoyama, Y. Yonamoto, T. Ohta, and A. Ugawa, Phys. Rev. B
the Cluster Variation and Path Probability Methods 共Plenum, 54, 6921 共1996兲.
24
New York, 1996兲. T. Yokoyama, K. Kobayashi, T. Ohta, and A. Ugawa, Phys. Rev.
4
R. Ahuja, J. M. Wills, B. Johansson, and O. Eriksson, Phys. Rev. B 53, 6111 共1996兲.
B 48, 16 269 共1993兲. 25
A. I. Frenkel and J. J. Rehr, Phys. Rev. B 48, 585 共1993兲.
5
S. A. Ostanin and V. Yu. Trubitsin, Phys. Rev. B 57, 13 485 26
N. Boccara and G. Sarma, Physics 共N.Y.兲 1, 219 共1965兲.
共1998兲. 27
N. S. Gillis, N. R. Werthamer, and T. R. Koehler, Phys. Rev. B
6
U. Pinsook and G. J. Ackland, Phys. Rev. B 58, 11 252 共1998兲; 165, 951 共1968兲.
59, 13 642 共1999兲. 28
P. F. Choquard, The Anharmonic Crystal 共Benjamin, New York
7
A. Zupan, P. Blaha, K. Schwarz, and J. P. Perdew, Phys. Rev. B 1967兲.
58, 11 266 共1998兲. 29
R. C. Shukla and E. Sternin, Philos. Mag. B 74, 1 共1996兲.
8 30
P. J. Wojtowaz and J. G. Kirkwood, J. Chem. Phys. 33, 1299 R. Hardy, M. A. Day, R. C. Shukla, and E. R. Cowley, Phys. Rev.
共1960兲. B 49, 8732 共1994兲.
9
M. Asta, R. McCormark, and D. de Fontaine, Phys. Rev. B 93, 31
R. C. Shukla, Philos. Mag. B 74, 13 共1996兲, Phys. Status Solidi B
748 共1993兲. 205, 481 共1998兲.
10
V. L. Moruzzi, J. F. Janak, and K. Schwarz, Phys. Rev. B 37, 790 32
F. Ducastelle, J. Phys. 共France兲 31, 1055 共1970兲.
共1988兲. 33
R. P. Gupta, Phys. Rev. B 23, 6265 共1981兲; D. Tomanenk, S.
11
A. P. Sutton, Philos. Mag. A 60, 147 共1989兲. Mukherjee, and K. H. Bennemann, ibid. 28, 665 共1983兲; W.
12
V. Ozolins, C. Wolrerton, and A. Zunger, Phys. Rev. B 57, 4816 Zhong, Y. S. Li, and D. Tomanek, ibid. 44, 13053 共1991兲.
共1998兲; 57, 6427 共1998兲; 58, R5897 共1998兲. 34
M. W. Finnis and J. E. Sinclair, Philos. Mag. A 50, 45 共1984兲.
13 35
K. Persson, M. Ekman, and G. Grimvall, Phys. Rev. B 60, 9999 V. Rosato, M. Guillope, and B. Legrand, Philos. Mag. A 59, 321
共1999兲. 共1989兲.
14
V. Ozolins and M. Asta, Phys. Rev. Lett. 86, 448 共2001兲. 36
F. Cleri and V. Rosato, Phys. Rev. B 48, 22 共1993兲.
15
Van de Walle and G. Ceder, Rev. Mod. Phys. 74, 11 共2002兲. 37
B. von Sydow, J. Hartford, and G. Wahnstrom, Comput. Mater.
16
S. Rubini and P. Ballone, Phys. Rev. B 48, 99 共1993兲. Sci. 15, 367 共1999兲.
17
G. P. Srivastava, The Physics of Phonons 共Hilger, Bristol, 1990兲. 38
C.-H. Chien, E. Blaisten-Barojas, and M. R. Pederson, J. Chem.
18
A. A. Maradudin, P. A. Flinn, and R. A. Coldwell-Horsfall, Ann. Phys. 112„5…, 2301 共2000兲.
Phys. 共N.Y.兲 15, 337 共1961兲; A. A. Maradudin and P. A. Flinn, 39
M. M. Sigalas and D. A. Papaconstantopoulos, Phys. Rev. B 49,
Phys. Rev. 129, 2529 共1963兲. 1574 共1994兲.

094301-13
K. MASUDA-JINDO, VU VAN HUNG, AND PHAM DINH TAM PHYSICAL REVIEW B 67, 094301 共2003兲

40 59
Y. Li, E. B. Barojas, and D. A. Papaconstantopoulos, Phys. Rev. B J. Mei, J. W. Davenport, and G. W. Fernando, Phys. Rev. B 43,
57, 15519 共1998兲. 4653 共1991兲.
41
M. S. Daw and M. I. Baskes, Phys. Rev. B 29, 6443 共1984兲. 60
N. I. Papanicolaou, G. C. Kallinteris, G. A. Evangelakis, and D.
42
S. M. Foiles, M. I. Baskes, and M. S. Daw, Phys. Rev. B 33, 7983 A. Papaconstantopoulos, Comput. Mater. Sci. 17, 224 共2000兲.
共1986兲. 61
R. C. Shukla and C. A. Plint, Phys. Rev. B 40, 10 337 共1989兲, and
43
M. I. Baskes, Phys. Rev. Lett. 59, 2666 共1987兲. references therein.
44
J. M. Holender, Phys. Rev. B 41, 8054 共1990兲. 62
M. Simerska, Acta Crystallogr. 14, 1259 共1961兲; Czech. J. Phys.
45
R. Pasianot, D. Farkas, and E. J. Savino, Phys. Rev. B 43, 6952 B12, 858 共1962兲.
共1991兲. 63
R. Kikuchi, Phys. Rev. 81, 998 共1951兲; T. Mohri, J. M. Sanchez,
46
L. A. Girifalco and V. G. Weizer, Phys. Rev. 114, 687 共1959兲. and D. de Fontaine, Acta Metall. 33, 1171 共1985兲.
47
N. V. Hung and J. J. Rehr, Phys. Rev. B 56, 43 共1997兲. 64
R. Kikuchi and K. Masuda-Jindo, Comput. Mater. Sci. 8, 1
48
U. Hansen, P. Vogl, and V. Fiorentini, Phys. Rev. B 60, 5055 共1997兲.
共1999兲. 65
A. Finel, Suppl. Prog. Theor. Phys. 115, 59 共1994兲; A. Finel and
49
N. Tang and V. V. Hung, Phys. Status Solidi B 149, 511 共1988兲; V. R. Tétot, in Stability of Materials, edited by A. Gonis, P.E.A.
E. Panhim et al., Phase Theory in Alloys 共Nauka, Novosibirk, Turchi, and J. Kudrnovsky, NATO ASI Series, Vol. 355 共Plenum,
1984兲 共in Russian兲. New York, 1996兲, p. 197.
50
N. Tang and V. V. Hung, Phys. Status Solidi B 162, 371 共1990兲. 66
R. A. MacDonald, R. C. Shukla, and D. K. Kahalner, Phys. Rev.
51
Vu Van Hung and K. Masuda-Jindo, J. Phys. Soc. Jpn. 69, 2067 B 29, 6489 共1984兲.
共2000兲; V. V. Hung, H. V. Tich, and K. Masuda-Jindo, ibid. 69, 67
D. Porezag, Th. Kohler, G. Seifert, and R. Kaschner, Phys. Rev. B
2691 共2000兲. 51, 12947 共1995兲; Th. Frauenheim, F. Weich, Th. Kohler, S.
52
D. Frenkel and A. J. C. Ladd, J. Chem. Phys. 81, 3188 共1984兲. Uhlmann, D. Porezag, and G. Seifert, Phys. Rev. B 52, 11492
53
D. Frenkel, Phys. Rev. Lett. 56, 858 共1986兲. 共1995兲.
54
Ya. P. Terletsky and N. Tang, Ann. Phys. 共Leipzig兲 19, 299 68
R. E. Cohen, M. J. Mehl, and D. A. Papaconstantopoulos, Phys.
共1967兲. Rev. B 50, 14 694 共1994兲; E. Wasserman, L. Stixrude, and R. E.
55
S. S. Mitra and S. K. Joshi, J. Chem. Phys. 34, 1462 共1961兲. Cohen, ibid. 53, 8296 共1996兲.
56
K. A. Gschneidner Jr., Solid State Phys. 16, 275 共1964兲. 69
P. Giannozzi, S. de Gironcoli, P. Pavone, and S. Baroni, Phys.
57
American Institute of Physics Handbook, 3rd ed., edited by B. H. Rev. B 43, 7231 共1991兲.
Billings et al. 共McGraw-Hill, New York, 1982兲. 70
D. J. Kouri, Y. Huang, and D. K. Hoffman, J. Phys. Chem. 100,
58
M. W. Zemansky, Heat and Thermodynamics 共McGraw-Hill, 7903 共1996兲; J. M. Thijssen and J. E. Inglesfield, Phys. Rev. B
New York, 1957兲, p. 263. 51, 17 988 共1995兲.

094301-14