j

2 downloads 0 Views 2MB Size Report
Parrinello (CP) [2] pointed out that the ionic and electronic degrees of freedom of a system need not to ... properties of K^Sbi-^ Zintl-phases in the liquid state. 2.
Previous Page

THE VIENNA AB-INITIO SIMULATION PROGRAM VASP: AN EFFICIENT AND VERSATILE TOOL FOR STUDYING THE STRUCTURAL, DYNAMIC, AND ELECTRONIC PROPERTIES OF MATERIALS

J.HAFNER and G. KRESSE

Institut filr Theoretische Physik and Center for Computational Materials Science Technische Universitdt Wien Wiedner Hauptstr. 8-10, A-1040 Wien, Austria

1. INTRODUCTION For more than thirty years the local-density-approximation (LDA) [1] has formed the basis of the progress towards an ab-initio description of complex processes in materials. A remarkable step on this way was realized eleven years ago when Car and Parrinello (CP) [2] pointed out that the ionic and electronic degrees of freedom of a system need not to be treated separately, but may be optimized simultaneously. The work of Car and Parrinello has triggered many fruitful developments and today we are able to treat systems and processes on the basis of fully quantum-mechanical simulations that are far more complex than has been thought to be possible before. However, especially for metals, some serious difficulties remain: (a) The original CP approach is based on a coupled pseudo-Newtonian dynamics for the ionic coordinates and for the electronic wavefunctions and requires a careful control of the adiabaticity of the system. In practice it turns out that adiabaticity can be maintained only if the ionic and electronic subsystems are effectively decoupled because the existence of an energy-gap at the Fermi-level leads to largely different time-scales for electronic and ionic excitations. For metals adiabaticity must be enforced, e.g. by draining the kinetic energy transferred from the ions to the electrons via a Nose thermostat attached to the electrons (the energy lost by the ions must than be restored by a second thermostat) [3]. (b) Most CP-programs are based on plane-wave (PW) basis-sets. This allows to take advantage of the highly efficient fast-Fourier-transform (FFT) technology for calculating the operation of the Hamiltonian on a wavefunction and avoids the appearence of Pulay-terms in the calculation of the Hellmann-Feynman forces acting

on the atoms. PW basis-sets can be used only in conjunction with a pseudopotentialdescription of the electron-ion interaction. Still, with conventional pseudopotentials the convergence of the PW-expansion is slow, especially for transition-metals and for first-row elements, (c) The calculation of the Kohn-Sham orbitals necessitates explicit orthonormalizations - these operations scale with the third power of the number of valence-electrons (C^A^-scaling) and hence rapidly determine the computational effort for large systems. In this paper we describe the Vienna ab-initio simulation program VASP [4, 5, 6] which was developed with the aim of bringing an efficient solution to these problems: (a) The problem of adiabaticity is eliminated altogether by calculating the electronic ground-state exactly after each ionic move. This can be achieved either via a global minimization of the Kohn-Sham (KS) total energy or via an iterative calculation of the occupied and some of the lowest unoccupied eigenstates of the KS Hamiltoninan (combined in an outer iteration-loop with an effective charge-density-mixing algorithm). A brief description and an estimation of the comparative efficiency of different techniques will be given in Sec. 2.3. (b) The electron-ion interaction is described by ultrasoft pseudopotentials which allow to optimize at the same time PW-convergence and transferability (see Sec. 2.1). (c) The impact of the O(N3) orthonormalization operations is minimized by choosing an appropriate strategy for the calculation of the KS-ground-state. We find that this is more easily achieved with iterative diagonalizations of the KS-Hamiltonian than with global minimization methods, (d) To achieve energy-conservation during the MD simulation it is necessary to use finite-temperature density-functional theory [7] where the electronic free energy is the conserved quantity. In addition, the broadening of the one-electron energy-levels leads to a faster convergence of Brillouin-zone integrations. We shall first review the basic principles of VASP and than describe exemplary applications to alloys and compounds: (a) the calculation of the elastic and dynamic properties of a metallic compound (€0812), (b) the surface reconstruction of a semiconducting compound (SiC), and (c) the calculation of the structural and electronic properties of K^Sbi-^ Zintl-phases in the liquid state.

2. TECHNICAL ASPECTS OF AB-INITIO MOLECULAR-DYNAMICS, ESPECIALLY FOR METALS 2.1 Ultrasoft pseudopotentials

Ultrasoft pseudopotentials have been introduced with the aim of optimizing simultaneously PW-convergence and transferability [8, 9]. Fast convergence is achieved by decomposing the charge-density p(r) into a smooth part described in terms of pseudoorbitals Qn (f) and localized augmentation charges P(r) = £ /„ I *n(r) I 2 + ^>n I ft)(ft I S n )^(O n

n

where fn is the occupancy of state n, the /3^s are localized projection states determining the nonlocal part of the pseudopotential according to VNL = ^DpIbWl ij

,

and the Qij(r) are localized augmentation functions. The indices i and j are a compact notation for the quantum numbers and reference energies for the valence states. By using more than one reference energy one can achieve a high accuracy of the pseudopotential over very large energy ranges and hence optimal transferability. The price to pay is that the pseudo-orbitals now satisfy a generalized orthonormality constraint ($n

S I $„) = £„„,

with the overlap matix S denned as

S = l + 5>,-|#)(/3,-| U

where the % = f Qij(r)d3r are the augmentation charges. The $ n 's also obey a generalized KS equation. Ultrasoft pseudopotentials have now been constructed for all elements for all elements from H to Bi [1O]. It has been shown that without any loss of accuracy even for first-row and transition-elements small PW basis-sets comparable in size to those necessary for soft normconserving pseudopotentials for Al or Si can be used. For any further details, see [9, 37]. 2.2 Finite-temperature LDA and &-space integration It has been shown [12, 13] that for metals energy-conservation along the MDtrajectory can be achieved only within the finite-temperature LDA formalism of Mermin [7]. At finite T the conserved quantity is the electronic free energy which is variational with respect to the charge-density />(r), the orbital occupancy / n , and the chemical potential jJ.(f). The FT-LDA formalism also puts the "smearing" methods introduced earlier with the aim of improving the convergence of Brillouin-zone integrations on a firm footing [14, 15]. Within VASP a variety of these techniques best adapted for different applications has been implemented. For dynamical simulations the Methfessel-Paxton [16] technique based on a set of special Appoints is the method of choice because it allows to calculate the Hellmann-Feynman forces directly as the derivatives of the free energy. For a comparative discussion of the relative merits of other approaches we refer to the recent paper by Kresse and Furthmuller [17]. 2.3 Determination of the LDA ground-state For the calculation of the LDA ground-state one can proceed either via the "direct" methods, i.e. via the glocal minimization of the total free energy with respect to the electronic degrees of freedom, or via the the diagonalization (for large PW basis-sets necessarily iterative diagonalization) of the KS Hamiltonian in combination with an iterative update of charge-density and potential. 2.3.1 Direct minimization of the total energy. A variety of minimization techniques have been proposed and implemented for ab-initio -MD simulations: (a) steepest-descent algorithms [18], (b) damped second-order equations-of-motion (EOM) [19], and (c) conjugate-gradient (CG) schemes [20, 21, 22]. The CG approach applied to all bands simultaneously (CGa) is parameter-free and normally as efficient as the second-order EOM, but requires the explicit orthonormalization of the conjugate gradient to all lower-lying bands. Hence this approach involves a large number of O(N3)

operations. A possible way to eliminate this difficulty is to generalize the KS functional to non-orthogonal orbitals [23, 24]. Here we shall not discuss any technical details (we refer to [17]), but we will present below detailed comparisons of the performance of CGa algorithms as compared to iterative diagonalization strategies. 2.3.2 Iterative matrix diagonalization. An important point to realize is that with PW basis-sets the dimension of the Hamiltonian matrix is typically of the order Natom x 100, whereas only 2 x Neiect (i.e. typically 10 percent) eigenstates are occupied by electrons. Hence iterative techniques for the calculation of the lowest eigenstates are infinitely more advantageous than a complete diagonalization. During the developement of VASP [4, 5, 6] a variety of different techniques for the iterative diagonalization of the KS Hamiltonian have been tested. The most efficients methods are (a) the sequential (i.e. band-by-band) CG algorithms originally proposed by Teter et al. [20] for total-energy minimization and later used by Bylander et al. [25] for the iterative calculation of the one-electron energies and (b) a novel variant of the residual-vector minimization method (RMM) first proposed by Wood and Zunger[27] on the basis of the direct inversion in the iterative subspace developed by Pulay [26]. Band-by-band conjugate-gradient (CG) techniques. The sequential CG technique minimizes the Rayleigh quotient

_ ($m I H I Sn.) ~ $m S $m

£app

for an approximate eigenstate | 4> m ). Variation of the Rayleigh quotient with respect to I 4> m ) yields the residual vector

\ R($m)) = (K-tappS)

\$m)

which must be zero for the exact eigenstate. Orthonormality of the eigenstates requires the gradient to be orthogonal to all other eigenvectors, convergence is improved by introducing a preconditioning function K to damp the corrections for the largewavevector components. Hence minimization is performed along the gradient vector

I g($m)) = (1 - Y, I *»(>*» S)K(H-e a p p S) $ m >

,

n

respectively the conjugate gradient vector (for details, see [17, 28]). In addition, for metals a diagonalization in the subspace spanned by the | 3> m ) (subspace rotation) must be performed. The main drawback of the otherwise highly efficient sequential CG algorithm is the need for the explicit orthonormalization of the gradients which not only introduces a large number of O(N3) operations, but also requires a large band-width of the communication between the main memory and the CPU. For both reasons, the efficiency decreases for larger systems. Residual minimization method (RMM-DIIS). Wood and Zunger [27] proposed to minimize the norm of the residual vector instead of the Rayleigh quotient. This is an unconstrained minimization condition. Each minimization step starts with the evaluation of the preconditioned residual vector K | R($lm)) for the approximate eigenstate I ^J n ). Then the norm of the residual vector is minimized in the subspace spanned by I $lm) and the new trial eigenstate is given by

I O = I O + A K I R(Qin))

.

Assuming the linearity of the residual vector, this is equivalent to the minimization of the quantity

EijC=C< bridges between alternating pairs of Si surface-atoms (see Fig. 2(a))[45]. ( b ) Staggered rows of C2 groups (with triple C=C bonds) in Si bridging sites (see Fig. 2(b))[44]. Recently Kackell et al.[46] have studied the stability of both models using VASP. C has been described by an ultrasoft pseudopotential with a cutoff-energy of 240 eV (see Furthmuller et al.[37]), Si by a normconserving pseudopotential. For model (a) the relaxation of the atomic geometry yields a dimer-length of dc-c — 1-38 A (i.e. very close to that of a C-C double bond in a hydrocarbon molecule), and an inwardrelaxation of the C top-layer by 0.29 A. For model (b) the bond-length within the C2 groups is dc-c — 1-23 A (corresponding to a triple-bond), the C top-layer relaxes outward by 0.24 A, Si atoms not connected by C2-bridges form dimers with a bondlength of dsi-si = 2.38 A. Energetically, both surfaces are nearly degenerate, with model (b) only a few meV/surface-atom lower in energy. Surface states within the

Figure 2: Models for the carbon-terminated 3C-SiC(OOl)c(2x2) surface: (a) staggered C=C-dimer model, (b) C 2 groups in Si bridge sites. Cf. text. bulk gap are predicted for the staggered dimer model (a), but not for the CVbriclge model (b). Recent photoemission data have not detected any surface states in the gap, thus confirming the preference for model (b). VASP has also been used to study hydrogenated of the C-terminated surface, as well as Si-terminated surfaces of various stochiometries[43, 46, 47]. 3.3 The structural and electronic properties of K-Sb Zintl phases in the crystalline and molten states The "Zintl"-phases formed by the alkali-metals and B-group metals or semi-metals from groups III to VI occupy a fascinating niche in the physics and chemistry of intermetallic compounds, because they display a wide variety of chemical bonding and electronic properties. Because of the large electronegativity difference one expects that the alkali atoms transfer their single valence electron to the electronegative element and the most stable compound is expected for the stoichiometric composition at which the octet-rule is satisfied. For the I-V compounds this leads to I3V-compounds with a saltlike crystal structure (BiLi3- or AsNa 3 -u/pe). In addition, at the equiatomic composition polyanionic compounds with a covalent character are formed: after charge-transfer, the V~-ions are isoelectronic to the chalcogen-atoms and form, according to the "Zintlprinciple" infinite V^0 spiral chains in analogy to the spiral crystal structures of the chalcogen-elements (NaP- and LiAs-type, see e.g. von Schnering et al. [48]). Note that both types of compounds are non-metallic, with an "ionic" gap in the octet-compounds separating the highest occupied anion-bancl from the lowest empty cation-band and a "covalent" gap in the polyanionic compounds separating anion-states that are bonding and anti-bonding along the spiral chains[49]. Many experiments show that the characteristic physical anomalies do not disappear on melting and suggest that the local order in the melt might be very similar to that in the corresponding crystalline compounds in particular the existence of covalently bonded Sb-chains in liquid alkali-Sb alloys has been suggested on the basis of neutron-scattering experiments[5O]. To confirm such a suggestion by detailed ab-initio molecular-dynamics simulations remains a challenging task, both because of the complexity of the interatomic interactions and the time required to reach local chemical equilibrium. Here we present first results of ab-initio LDA-MD studies of crystalline and liquid K-Sb alloys[51]. The liquid alloys are modelled by 64-atom ensembles of appropriate

DOS [1/eV atom]

Figure 3: (a) Partial pair correlation functions gij(R) in liquid K-Sb alloys, (b) Total, partial, and local electronic densities of states in liquid Ko.soSbo.so- Cf. text. composition, the dynamical simulations are performed in a canonical ensemble with a Nose-thermostat, using a time-increment of 3 fs (with the exact LDA ground-state and Hellmann-Feynman forces calculated after each MD-step). A characteristic property of the Zintl-compounds is the large negative excess-volume which is not known very accurately for the liquid alloys. First, we determined the zero-pressure density at each composition by a series of short const ant-volume MD runs. At the equilibrium density, at each composition MD runs over about 10 ps where performed (for details we refer to the work on pure liquid Sb [52]). Fig. 3(a) summarizes the results for the atomic structure. The partial pair correlation functions show that the dominant features of the liquid structure are: (a) A strong tendency to heterocoordination of the Sb-ions with short K-Sb distances, almos. no direct Sb-Sb neighbours at the octet-composition KsSb, and a Sb-Sb coordination number of Nsb-sb ~ 2 at the equiatomic composition, (b) Together with strong, short Sb-Sb bonds this confirms the existence of broken and entangled Sb-Sb chains in the liquid alloys. Further confirmation of this picture comes from the electronic structure. Fig. 3(b) demonstrates that the characteristic features of the crystalline band-structure are preserved in the liquid phase l : The valence band is split into a low-lying Sb-s band of one-dimensional character and bonding and non-bonding (ppcr) bands just below the Fermi level. The Fermi-energy falls into a "pseudo-gap" between the non-bonding and the anti-bonding (ppcr) bands - the crystalline KSb compound is semiconductor with a, narrow gap of 0.7 eV (we only note in passing that the VASP results for the crystalline compounds are in perfect agreement with the earlier all-electron calculations [49]). The existence of a pseudogap for the liquid alloys explains the electronic transport anomalies observed at this composition[54]. The close similarity of the electronic structure 1

TLe local and partial densities of state have been calculated using projections of the plane-wave components of the eigenstates onto spherical waves centred at the atomic sites[53].

of crystalline and molten KSb with that of crystalline and liquid Se is a very striking illustration of the validity of the "Zintl-principle".

4. CONCLUSIONS In this article we have reviewed very briefly the basic ingredients that make the Vienna ab-initio simulation program VASP a technique of choice for studying the properties of and complex processes in materials: (a) The use of ultrasoft pseudopotentials allows to treat transition-metals and first-row elements at a computational effort that is not much bigger than for simple metals and semiconductors, (b) The exact calculation of the LDA ground-state after each ionic move eliminates all non-adiabaticity problems - this is especially important for systems close to a metal/nonmetal transition, (c) The impact of the cumbersome orthogonalization operations on the scaling of the computational effort with the number of electrons is minimized by using iterative diagonalization via the residual-minimization-method RMM, together with optimized mixing procedures, for the calculation of the ground state, (d) Finite-temperature LDA helps to improve energy-conservation and to speed-up the convergence of Brillouin-zone integrations. The performance of VASP for alloys and compounds has been illustrated at three examples: The calculation of the properties of cobalt dislicide demonstrates that even for a transition-metal compound perfect agreement with all-electron calculations may be achieved at much lower computational effort, and that elastic and dynamic properties may be predicted accurately even for metallic systems with rather long-range interactions. Applications to surface-problems have been described at the example of the 3C-SiC(IOO) surface. Surface physir.s and catalysis will be a particularly important field for the application of VASP, recent work extends to processes as complex as the adsorption of thiopene molecules on the surface of transition-metal sulfides[55]. Finally, the efficiciency of VASP for studying complex melts has been illustrate for crystalline and molten Zintl-phases of alkali-group V alloys.

ACKNOWLEDGEMENTS The contributions of J. Furthimiller, P. Kackell, K. Seifert, R. Stadler, and R. Poclloucky to various parts of the work described in this article is gratefully acknowledged. Part of this work has been supported by the Bundesministerium fur Wissenschaft, Forschung und Kunst through the Center for Computational Materials Science.

References [1] W. Kohn and L.J. Sham, Phys. Rev. 140 (1965) A1133. [2] R. Car and M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471. [3] P.E. Blochl and M. Parrinello, Phys. Rev. B 45 (1992) 9412. [4] G. Kresse and J. Hafner, Phys. Rev. B 47 (1993) RC558. [5] G. Kresse and J. Hafner, Phys. Rev. B 48 (1993) 13115. [6] G. Kresse and J. Hafner, Phys. Rev. B 49 (1994) 14251.

[7] N.D. Mermin, Phys. Rev. 137 (1965) A1441. [8] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892. [9] G. Kresse and J. Hafner, J. Phys.: Condens. Matter 6 (1994) 8245. [10] G. Kresse, unpublished. [11] J. Furthmuller, J. Hafner, and G. Kresse, Phys. Rev. B 50 (1994) 15606. [12] M. Weinert and J. W. Davenport, Phys. Rev. B 45 (1992) 13709. [13] R.M. Wentzcovitch, J.L. Martins, and P.B. Allen, Phys. Rev. B 45 (1992) 11372. [14] C.L. Fu and K.M. Ho, Phys. Rev. B 28 (1983) 5480. [15] A. de Vita and MJ. Gillan, J. Phys.: Condens. Matter 3 (1991) 6225. [16] M. Methfessel and A.T. Paxton, Phys. Rev. B 40 (1989) 3616. [17] G. Kresse and J. Furthmuller, Comput. Mat. Sd. (in print). [18] A.R. Williams and J. Soler, Bull. Arner. Phys. Soc. B 32 (1987) 562; G.X. Qian, M. Weinert, G.W. Fernando, and J.W. Davenport, Phys. Rev. Lett. 64 (1990) 1146. [19] F. Tassone, F. Mauri, and R. Car, Phys. Rev. B *0 (1994) 10561. [20] M.P. Teter, M.C. Payne, and D.C. Allan, Phys. Rev. B 40 (1989) 12255. [21] I. Stich, R. Car, M. Parrinello, and S. Baroni, Phys. Rev. B 39 (1989) 4997. [22] MJ. Gillan, J. Phys.: Condens. Matter 1 (1989) 689. [23] T.A. Arias, M.C. Payne, J.D. Joannopoulos, Phys. Rev. Lett. 69 (1992) 1077. [24] J.M. Holender, MJ. Gillan, M.C.Payne, and A.D.Simpson, Phys. Rev. B 52 (1995) 967. [25] D.M. Bylander, L. Kleinman, and S. Lee, Phys. Rev. B 42 (1990) 1394. [26] P. Pulay, Chem. Phys. Lett. 73 (1980) 393. [27] D. M. Wood and A. Zunger, J. Phys. A: Math. Gen. Phys. 18 (1985) 1343. [28] G. Kresse and J. Furthmuller, Phys. Rev. B. (in print). [29] D. Vanderbilt and S.G. Louie, Phys. Rev. B 30 (1984) 6118. [30] J.F. Annett, Comput. Mat. Sd. 4 (1995) 23. [31] G. Kresse, Proceedings 9th Intern. Conference in Liquid and Amorphous Metals (LAM9), Chicago 1995 (in print). [32] R. Stadler, W. Wolf, R. Podloucky, G. Kresse, J. Furthmuller, and J. Hafner, Phys. Rev. B 54 (15. Aug.1996). [33] L. Hedin and B.I. Lundqvist, J. Phys. C 4 (1971) 2064.

[34] A.D. Becke, Phys. Rev. A 38 (1988) 3098; J.P. Perdew, Phys. Rev. B 33 (1986) 8822. [35] G. Kresse, J. Furthmuller, and J. Hafner, Europhys. Lett, bf 32 (1995) 729. [36] P.H. Giauque, Travail de diplorhe, Institut de Physique Experimental de I'Universite de Lausanne 1992 (unpublished). [37] J. Furthmuller, J. Hafner, and G. Kresse, Europhys. Lett. 28 (1994) 659; Phys. Rev. B 53 (1996) 7334. [38] G. Kern, J. Hafner and G. Kresse, Surf. Sci. 352-4 (1996) 745; 357-8 (1996) 422. [39] G. Kern, J. Hafner and G. Kresse, Surf. Sd. (in print). [40] C. Cheng, R. J. Needs, and V. Heine, J. Phys. C 21 (1988) 1049. [41] P. Kackell, B. Wenzien, and F. Bechstedt, Phys. Rev. B 50 (1994) 17037; ibid. p. 10761. [42] B. Wenzien, P. Kackell, and F. Bechstedt, Surf. Sci. 331-3 (1995) 1105. [43] P. Kackell, J. Furthmuller, and F. Bechstedt, Surf. Sci. 352-4 (1996) 55. [44] J.M. Powers, A. Wander, P.J. Rous, M.A. van Hove, and G.A. Somorjai, Phys. Rev. B 44 (1991) 11159. [45] V.M. Bermudez and R. Kaplan, Phys. Rev. B 44 (1991) 11149. [46] P. Kackell, J. Furthmuller, F. Bechstedt, G. Kresse, and J. Hafner, Phys. Rev. B (in print). [47] P. Kackell, J. Furthmuller, an " F. Bechstedt, Appl. Surf. Sci. (in print). [48] H.G. von Schnering, W. H6nle, and G. Krogull, Z. Naturforsch. (b) 34 (1979) 1678. [49] M. Tegze and J. Hafner, J. Phys.: Condens. Matter 4 (1992) 2449. [50] P. Lamparter, W. Martin, S. Steeb, and W. Freyland, J. Non-Cryst. Solids 61+62 (1984) 279. [51] K. Seifert, Diplomarbeit, Technische Universitdt Wien 1995 (unpublished); K. Seifert, J. Hafner, and G. Kresse, 3rd Liquid Matter Conference, Norwich 1996 . [52] K. Seifert, Proceedings 9th Intern. Conference on Liquid and Amorphous Metals (LAM9), Chicago 1995 (in print). [53] A. Eichler, J. Hafner, and G. Kresse,-Surf. Sci. 346 (1996) 300. [54] H. Redslob, G. Steinleitner, and W. Freyland, Z. Naturforsch. (a) 2 (1982) 587. [55] P. Raybaud, G. Kresse, J. Hafner, and H. Toulhoat, to be published.

Configurational Kinetics Studied by PPM

Tetsuo Mohri Division of Materials Science and Engineering, Graduate School of Engineering, Hokkaido University Kita-13 Nishi-8, Kita-ku, Sapporo 060 Japan

1. Introduction As is well recognized, various macroscopic properties such as mechanical properties are controlled by microstructure, and the stability of a phase which consists of each microstructure is essentially the subject of electronic structure calculation and statistical mechanics of atomic configuration. The main subject focused in this article is configurational thermodynamics and kinetics in the atomic level, but we start with a brief review of the stability of microstructure, which also poses the configurational problem in the different hierarchy of scale. Transmission electron microscopic observaiion reveals various morphologies of precipitates depending upon the constituents and composition of an alloy system, history of heat treatments etc. Typical examples are spherical precipitates found in Ni-Cr-Al system and cuboidal precipitates in Fe-Mo system [I]. The first question raised is what determines the shape of a precipitate. A theoretical framework for this problem was provided by Nabarro [2] nearly half century ago. In his celebrated work, he calculated strain energy, Estrain, of an ellipsoid as a function of c/a ratio where a and c are, respectively, the major and minor axes of an ellipsoid. He demonstrated that a spherical precipitate induces highest strain energy while a minimum energy is reallized by a plate like precipitate. Under a constant volume, however, the total surface area is smallest for a sphere. Eventually a balance between the surface energy, Esurface, and elastic strain energy, Estrain, determines an optimum shape of the precipitate, which may claim the minimization of the total energy of the system Etotal ~ Estrain + Esurface.

(1)

It is reallized that Nabarro's theory deals with a single precipitate in the matrix and interactions among precipitates are not considered. Therefore, the configurational problem i.e. an alignment of precipitates are beyond the scope of the theory. The most striking phenomenon which manifests the interactions among precipitates is the particle splitting observed for 7' precipitates of Ni-Al system by Miyazaki et al [3] about a decade ago. Opposite to well known coarsening phenomenon, a precipitate splits into several pieces. Such is a typical example dictating the importance of the elastic interaction among particles which overcomes the increase of surface energy. Theoretical studies including computer

simulation have been advanced mainly by Khachaturyan and his coworkers [4], and they revealed various roles played by elastic interactions among precipitates. It is emphasized that configuration of microstructure is essentially the subject of energetics. The optimization of various elastic energy contributions dominate a whole realm. The stable configuration and evolution kinetics in the microstructural level are, therefore, controlled and determined through variational process operated on the total energy including elastic interaction energy, Eint, E total = E gtrain + E surf ace + Eint-

(2)

Contrary to the microstructural level, the entropy plays a significant role in the stability and kinetics of atomic level which is the main concern of this article. Among potential theoretical tools, Cluster Variation Method (hereafter CVM) [5] has been recognized as one of the most efficient ways of evaluating and incorporating entropy contribution. In fact, once the analytical expression of the free energy functional based on the CVM entropy formula is obtained, one can derive various thermodynamic quantities such as a phase diagram, instability (Spinodal ordering ) temperature [6], short range order diffuse intensity spectrum [7, 8] from a single free energy formula in a consistent manner [9]. Recent progress of electronic structure calculation even enables us to perform a first-principles calculation of the phase stability. Path Probability Method (hereafter PPM) [10, 11, 12, 13] is recognized as a kinetic version of the CVM, namely the extension of the variational principle into the time-space coordinate. Hence, various fruitful features are inherited from the CVM. Among them is the capability of incorporating wide range of atomic correlations associated with phase transitions, and the ordering and disordering kinetics which claim spatially extended atomic correlations are studied with sufficient accuracy. A noticable feature of the PPM is the fact that state variables (cluster probabilities or correlation functions of the CVM) and their time derivatives are not explicitly dealt, which is in marked difference with other kinetic theories. Also, the fact that relevant quantities derived for the long time limit correctly converge to the equilibriur ones predicted by the CVM ensures the reliability of the PPM. Hence the combination of the CVM and PPM provides invaluable tool to study thermodynamic properties of a given system starting from non-equilibrium to equilibrium states. In the present article, the main subject is the atomic configurational kinetics studied by the PPM. It should be realized that unlike the study of equilibrium thermodynamics for which a model is often mapped onto Ising system, elementary mechanism of atomic motion plays a deterministic role in the kinetic study. In an actual alloy system, diffusion of an atomic species is mainly driven by vacancy mechanism. The incorporation of the vacancy mechanism into PPM formalism, however, is not readily achieved, since the abundant freedom of microscopic path of atomic movement demands intractable number of variational parameters. The present study is, therefore, limited to a simple spin kinetics, known as Glauber dynamics [14] for which flipping events at fixed lattice points drive the phase transition. Hence, the present study for a spin system is regarded as a precursor to an alloy kinetics. The limitation of the model is critically examined and pointed out in the subsequent sections. Since a considerable amount of review articles on both theoretical frameworks and calculated results have been reported[15-25], the main objective of the present study is placed on the comparisons with experimental results. The organization of the present report is as follows: In the next section, for the sake of completeness, a brief theoretical description of the PPM is summarized from the previous articles. In the third section, disorder-L\Q transition is focused and visualized atomic (spin) configuration is compared with recent high resolution electron micrograph. In the fourth section, ordering relaxation

in a single phase field of Ll 0 and Ll 2 phases are studied. The derived long range order parameters are converted into the electrical resistivity change and compared with recent experimental measurements.

2. Path Probability Method In the CVM, the free energy of a given alloy is approximated in terms of probabilities for a selected set of finite clusters. The largest cluster explicitly considered in the free energy functional specifies the level of the approximation. The common practice for an feebased system is the tetrahedron approximation [26] in which nearest neighbor tetrahedron cluster is taken as the largest cluster. Hence, within the tetrahedron approximation, the free energy expression, F,is symbolically expressed as z?

1 77/V J^l L,ss'\ loss's"} T X Z

F =F

V ' \ iJ '

fa

J ' \

^

Lltss's"s"'\\ J » \Wijkl J) ,

/o\ (3)

where z',^', &, / indicate atomic species (up and down spins) , and #f, yfj , z%jk6 and w^ 6 represent, respectively, the cluster probability of atomic configuration specified by subscript (s) on a point, nearest neighbor pair, nearest neighbor triangle and nearest neighbor tetrahedron clusters located on the sublattice designated by superscript (s). Note that all the cluster probabilities are not independent but are interrelated each other through geometric relationships such as ^S X

_ V V ?X -

W

7SS'S" -

~~ 2^ Z^ Zijk S1S" jk

i ~ Z^ Z^ 2/*> 6' j



V

V onSS'S"S'"

2-f Z^ Wijkl S'S"S'" jkl

(A}



(4)

It is remembered that atomic correlation generally exceeds the range of atomic interaction. Hence, in constructing a free energy formula, it is necessary to limit the range of atomic interaction within the spatial extension of the largest cluster assumed in the entropy term S. In the present study, we adopt nearest neighbor pair interaction model for an internal energy term E. The free energy is, therefore, given as F

=

E-TS

^Ve if* -- ^ou6^ V 2^~j~ 2^e^yij SS' ° ij l

-T [- (I) E "I E L (*?) + E < E L (#) - 2 E L ((*?(*)) (ij

S

i

)

+£< (E L№) -EL (#(*))} 66'

-2 I

\ij>kl

£

{^ijkl^mnop

ij

)

L (Wijkl,mnop) -^L (wi]kl(t)) 1 , ijkl

(10)

J

where 0 is the spin flip probability per unit time, Xfj, Y/^i and Wijk^mnop are the path variables for the flipping from one spin configuration to another designated by the subscript^) before (at time t) and after (at time t H- Atf) the comma sign, on a point, pair and tetrahedron clusters over the specified sublattice(s), respectively, and AE* is the change of the internal energy before and after the flipping events.

In the PPF, the first factor PI describes the statistical average of non-correlated spin flip events over entire lattice points, and the second factor P2 is the conventional thermal activation factor. Hence, the product of PI and P2 corresponds to the Boltzmann factor in the free energy and gives the probability that on^ of the paths specified by a set of path variables occurs. The third factor P3 characterizes the PPM. One may see the similarity with the configurational entropy term of the CVM (see eq.(5)), which gives the multiplicity, i.e. the number of equivalent states. In a similar sense, P3 can be viewed as the number of equivalent paths, i.e. the degrees of freedom of the microscopic evolution from one state to another. As was pointed out in the Introduction section, mathematical representation of P3 depends on the mechanism of elementary kinetics. It is noted that eqs.(S)-(IO) are valid only for a spin kinetics. The path variables of the PPM corresponds to the cluster probabilities of the CVM by which the free energy is minimized to obtain the most probable state. Likewise, under a set of constraints, the PPF is maximized with respect to the path variables for each time step, which yields the optimized set of path variables. Since a set of path variables, B0^(£;£ -f At), relates cluster probabilities x^(t) at time t and x(t + At) at time t + A* through A

x+(t + A*) = x+(t) + E s**(