The Diffusion Mechanism in Amorphous
Ni
81
B
19
Studied by Molecular
Dynamics Simulations
Ni
81
B
19
Studied by Molecular
The Diffusion Mechanism in Amorphous
Ni
81
B
19
Studied by Molecular
Dynamics Simulations
Proefschrift
ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,
op gezag van de Rector Magnificus prof. ir. K.F. Wakker, in het openbaar te verdedigen ten overstaan van een commissie,
door het College voor Promoties aangewezen, op vrijdag 6 november 1998 te 10.30 uur
door
Léonardus Daniël VAN EE
materiaalkundig ingenieur geboren te Goirle
Prof. dr. ir. A. van den Beukel.
Toegevoegd promotor: Dr. B.J. Thijsse.
Samenstelling promotiecommisie:
Rector Magnificus Voorzitter
Prof. dr. ir. A. van den Beukel Technische Universiteit Delft, promotor
Dr. B.J. Thijsse Technische Universiteit Delft, toegevoegd promotor
Prof. dr. D. Frenkel Universiteit Utrecht en FOM-AMOLF
Prof. dr. S.W. de Leeuw Technische Universiteit Delft
Prof. dr. J.Th.M. de Hosson Rijksuniversiteit Groningen
Prof. dr. ir. E. van der Giessen Technische Universiteit Delft
Dr. ir. J. Sietsma Technische Universiteit Delft
Dr. ir. J. Sietsma heeft als begeleider in belangrijke mate aan de totstandkoming van het proefschrift bijgedragen.
This work is part of the research program of the Stichting voor Fundamenteel Onderzoek der Materie (Foundation for Fundamental Research of Matter), and was made possible by financial support from the Nederlandse Organisatie voor
Wetenschappelijk Onderzoek (Netherlands Organization for Scientific Research).
ISBN 90-9012078-5
Contents
1. General introduction ... 1
1.1 Atomic motion in metallic glasses...1
1.2 Mechanism of diffusion...2
1.3 Molecular dynamics simulations of amorphous Ni81B19...2
1.4 Scope of the investigations ...3
References ...5
2. Computer model of amorphous Ni81B19
... 7
2.1 Introduction ...7
2.2 Experimental data on the structure ...9
2.3 Molecular dynamics algorithm ...11
2.4 Potential energy function...13
2.4.1 Adiabatic approximation ...13
2.4.2 Pair potentials ...14
2.5 Starting structure ...15
2.6 Simulation details ...18
2.6.1 Constants...18
2.6.2 Temperatures and time scales ...19
2.6.3 Equilibrium configurations ...21
2.7 Model characteristics ...21
2.7.1 Glass transition ...21
2.7.2 Radial distribution functions...24
2.7.3 Vibrational density of states ...25
2.8 Conclusions ...26
3. Atomic motion ... 29
3.1 Introduction ...29
3.2 Mean square drift...30
3.2.1 In-cage movement ...30
3.2.2 Diffusion...31
3.3 Van Hove probability function ...36
3.4 Non-Gaussian parameter...40
3.5 Chain jumps ...42
3.6 Conclusions ...47
References ...47
4. Temperature dependence of structural defects
and diffusion... 49
4.1 Introduction ...49
4.2 Diffusion and defects ...50
4.3 Identification of interatomic space ...51
4.4 Void volumes and concentrations...53
4.5 Void shape...56
4.6 Defects ...58
4.7 Diffusion ...62
4.8 Conclusions ...64
References ...64
5. Vibrational modes and localization ... 65
5.1 Introduction ...65 5.2 Vibrational theory ...66 5.2.1 Definitions ...66 5.2.2 Harmonic approximation...67 5.2.3 Dynamical matrix...68 5.2.4 Vibrational modes ...70
5.2.5 Periodicity and wavevectors ...71
5.2.6 Mode and atomic oscillations...72
Contents vii
5.4 Localization ...75
5.4.1 Characterization...75
5.4.2 Distribution of amplitudes ...77
5.4.3 Local structure of low-frequency localized modes ...80
5.4.4 Localization and interatomic space ...82
5.5 Conclusions ...83
References ...84
Appendix A ...85
6. Atomic two-level states and relaxations ... 87
6.1 Introduction ...87
6.2 Launching a vibrational mode ...89
6.3 Atomic motion after the launching ...91
6.4 Reversibility ...97
6.5 Relaxations and two-level states ...98
6.6 Atomic and vibrational mode energies ...99
6.7 Conclusions ...102
References ...103
7. A mechanism for atomic transport ... 105
7.1 Introduction ...105
7.2 Atomic motion and structural changes ...106
7.3 Localization...109
7.4 Correlation between localized modes and jumps...112
7.5 Exciting low-frequency localized vibrational modes ...113
7.6 Conclusions ...115
References ...115
8. General Conclusions... 117
1.
chapter 1
General introduction
1.1 Atomic motion in metallic glasses
This thesis deals with atomic mobility and the nature of the diffusion process in metallic glasses. Metallic glasses are usually obtained by quenching a metallic, multicomponent liquid fast enough to suppress crystallization. The liquid then enters the non-equilibrium amorphous state, which lacks the long-range order of crystals. As will be shown, the diffusion process in a glass is a very complex phenomenon. Although much is known experimentally about structural relaxation [1], diffusion [2,3], and plastic deformation [4,5], the dynamic atomic events responsible for these processes have not yet been captured in a comprehensive classification scheme, let alone in a full theory. In crystalline structures defects (including interstitial atoms, vacancies, and grain boundaries) act as vehicles for the diffusion process. Such a defect can be identified -and can survive- because of the high degree of order in a crystalline structure. In liquids, such defects are not present, but diffusion proceeds through continuous correlated motions of groups of atoms. In metallic glasses neither of these mechanisms can be envisaged to take place. The local energy barriers are too high for hydrodynamic, liquid-like diffusion, but on the other hand, the densely-packed, disordered structure lacks defects similar to vacancies or grain boundaries. Diffusion in glasses is a process with a character of its own. The principal goal of the present study is to identify the essential mechanisms of atomic motion and to find out if there are specific locations in the structure (“defects” in a broader sense), where diffusional dynamic events, or “jumps” are inclined to happen. Experimental methods are only of limited value, because the microscopic information that we are looking for is not always directly accessible. Therefore, in this work molecular dynamics studies have been performed on a realistic computer model of a metallic glass at tempera-tures below and above the glass transition temperature. For reasons to be explaned
1.2 Mechanism of diffusion
The problem of unraveling the nature of the diffusion process in metallic glasses has been approached by various authors, both by means of experiments [1,6] and by means of computer modeling [7-11]. Some experiments have been aimed at the temperature dependence of the diffusivity and its sensitivity to structural relaxation [12], whereas other studies focus specifically on the mechanism of the diffusion process [13]. In the latter studies, however, difficulties are encountered in the interpretation. Two important sources of information on the mechanism of diffusion, namely the pressure dependence and the isotope dependence of the diffusivity, can be interpreted in an approximative way only, based on assumed analogies with crystalline matter [2]. For some alloys the pressure dependence of the diffusivity indicates that the defects involved in the diffusion process are associated with a so-called activation volume of approximately one atomic volume [3], for others the activation volume is found to be zero [2]. A zero activation volume indicates that the diffusion process is direct, i.e. that it does not take place at distinguishable sites in the structure. On the other hand, an activation volume of about one atomic volume indicates that the diffusion process is indirect, taking place at distinct sites in the structure. It is tempting to interpret an activation volume of one atomic volume as indicative for a single-atom diffusion process, involving vacancy-like defects. However, the isotope dependence of the diffusivity
[13] indicates that diffusion in amorphous Co77Fe2Nb14B7 takes place in a
cooperative manner, i.e. it involves a number of atoms moving simultaneously to their new positions. In Ref. [13] it is shown that the characteristics of this process depend on the state of structural relaxation of the glass. Recently, new light is shed on this problem by Teichler [14]. He performed molecular dynamics simulations on
amorphous Ni80Zr20 and Ni50Zr50. Experimental results had suggested that diffusion
is governed by cooperative atomic motion in Ni80Zr20 [15], and by single-atom
motion in Ni50Zr50 [3]. The atomic motion observed in the simulations is in
agreement with the experiments, and it is therefore very likely that, depending on the amorphous material under investigation, either the cooperative or the single-atom mechanism for diffusion prevails.
1.3 Molecular dynamics simulations of amorphous Ni
81B
19Molecular dynamics (MD) simulation is a well-established method of studying the details of atomic motion on a nanosecond and smaller time scale, which can provide valuable complementary information to experimental studies. The choice to
model amorphous Ni81B19 was directed by methodical constraints and the
General introduction 3
simple liquids (argon-like fluids) requires quenching rates that are not feasible, i.e. simple liquids are not “glass-formers”. The existence of a liquid-to-glass transition for such fluids has only been demonstrated by MD simulations [16], in which the quenching rate is typically 6 orders of magnitude faster than the highest experimen-tal quenching rates. These simulated simple glasses are often reported to crysexperimen-tallize easily [8,16]. This unwanted effect can be circumvented by studying a binary or more complex mixture of atoms with different diameters. In this case the crystallization is “frustrated” and will not occur within the time span that is currently available to MD simulations (about 10 ns). In order to simulate a system with k components with MD, the k(k+1)/2 possible types of atomic pair interaction must be known. Studying a three-component system, such as PdNiP, therefore requires that the six pair potentials are known, posing a complex problem. A two-component system is the optimum choice. If the mismatch in atomic sizes is larger than about 20% the phase diagram of the system can be expected to display an eutectic, a feature that greatly enhances the glass-forming ability [7]. The size mismatch between a Ni and a B atom is 32% [17] and the eutectic composition of the NiB system is 81 at. % Ni : 19 at. % B [18].
Research performed by others on amorphous Ni81B19 has led to a number of
results that can be used to assess to what extend the MD-model realistically represents a real glass: The experimental partial reduced radial distribution functions have been determined by Lamparter et al. [17], and the vibrational density of states by Lustig et al. [19]. The atomic interaction potentials have been studied by Hausleitner et al. [20] and Li et al. [21]. The static atomic structure of
amorphous Ni81B19 has been determined by Iparraguirre et al. [22], applying the
reverse Monte Carlo technique.
1.4 Scope of the investigations
The molecular dynamics model of amorphous Ni81B19 is introduced in chapter
2. It consists of 2500 atoms in a cubic box with periodic boundaries. The MD-simulation runs are carried out at ambient pressure and at fourteen constant temperatures between 304 and 1416 K. The equations of motion are solved using the velocity-Verlet algorithm [23] using effective pair potentials of the Weber-Stillinger type [24]. It is shown that the model system undergoes a glass transition
at Tg ≈ 590 K, that it has realistic glassy properties, and that it is in good
The present analysis of the atomic mobility can be separated into three parts: (i) a kinematic analysis of the atomic motion, (ii) a static analysis of amorphous configurations, and (iii) a dynamic analysis of the atomic motion.
The kinematic analysis, presented in chapter 3, shows that the atomic motion in the glassy state can be separated into two contributions: (i) the vibrational movement of an atom around its potential-energy minimum, and (ii) the occasional jump of an atom from its energy minimum into a neighboring potential-energy minimum (diffusion). The diffusion coefficients are determined from the atomic displacements, and the temperature dependence of the diffusivity above the glass-transition temperature is shown to be consistent with both a Vogel-Fulcher relation and a power law. Below the glass-transition temperature an Arrhenian temperature dependence is observed. A special kind of cooperative atomic motion, viz. atoms jumping like a chain, is shown to take place in the model.
In theories on the atomic mobility in metallic glasses (for instance the free-volume theory [25]) a key role is played by the concept of “structural defects”, which are not defined in geometrical terms, but merely denote sites in the structure that can act as catalysts in the diffusion process. The Arrhenian temperature
dependence of the diffusivity below Tg shows that diffusion is a thermally activated
process, and therefore diffusion is assumed to take place in distinct events that occur at these “defects”. However, in an amorphous structure, contrary to crystal-line structures, it is far from straightforward to define a defect. In chapter 4 a method will be used to identify these defects on the basis of a geometrical characterization of the interatomic space. The thus determined defect concentration
at room temperature is on the order of 10–5. The defects have a highly
non-spherical shape and consist of spaces situated between nine or more atoms. Below the glass-transition temperature the defect concentration shows an Arrhenian temperature dependence. Combining the diffusivity data with the defect concentra-tions allows the determination of the jump rate frequency and the migration energies.
In the third approach the dynamic behavior of the amorphous structure, i.e. the actual cause of the atomic motion, is studied. When for an atomic configuration of a disordered solid (quenched to the ground state at 0 K) the vibrational modes are calculated in the harmonic approximation, one often finds low-frequency modes with a low occupation ratio. This means that certain low-frequency vibrations have only a few atoms significantly participating in them, and that they are therefore
localized. Since the low frequency indicates a small curvature in the potential
General introduction 5
void-like defect, is easily made, as described in chapter 5. Consequently, using the harmonic analysis to search for localized low-frequency vibrational modes is synonymous with searching for sites where diffusion events may occur most easily.
In chapter 6 it will been shown that when, starting at 0 K, a single
low-frequency localized mode in amorphous Ni81B19 is excited with sufficient excitation
energy, the system can make a transition from its potential-energy minimum to a neighboring minimum. Such a “jump” is a cooperative process involving some tens of atoms. The jump can be either reversible or irreversible. After an irreversible jump the system contains other locations with similar dynamic properties, which may be the sites for new jump events. The atoms involved in a jump have more than one potential-energy minimum available and the jump location can be called a
two-level state. Two-level states play an important role in the low-temperature
physics of amorphous materials [26].
In chapter 7 it is demonstrated that jumps related to two-level states also take place if all modes are excited (at elevated temperatures), and that they are relevant
for the diffusion process. The model of amorphous Ni81B19 at a temperature just
below the glass-transition temperature is searched for locations where structural changes as a result of diffusion events take place. The vibrational motion of the atoms can be separated from the structural changes by the determination of the atomic equilibrium positions by energy minimization. Close inspection of these configurations shows that structural changes are highly localized, discrete events, and that they are correlated with low-frequency localized vibrational modes. Excitation of such a mode in the way described in chapter 6 can cause a jump resulting in the observed structural change.
References
[1] G. Ruitenberg, P. de Hey, F. Sommer and J. Sietsma, Phys. Rev. Lett. 79,
4830 (1997).
[2] F. Faupel, P.W. Hüppe and K. Rätzke, Phys. Rev. Lett. 65, 1219 (1990).
[3] H.J. Höfler, R.S. Averback, G. Rummel and H. Mehrer, Phil. Mag. Lett. 66,
301 (1992).
[4] G.W. Tecza and J. Hafner, J. Non-Cryst. Solids 202, 1 (1996).
[5] A.I. Taub and F. Spaepen, J. Materials Science 16, 3087 (1981).
[6] S.K. Sharma, M.P. Macht and V. Naundorf, Acta Metall. Mater 40, 2439
[7] Y. Hiwatari, H. Miyagawa, and T. Odagaki, Solid State Ionics 47, 179 (1991).
[8] J.L. Barrat and M.L. Klein, Annu. Rev. Phys. Chem. 42, 23 (1991).
[9] G. Wahnström, Phys. Rev. A 44, 3752 (1991).
[10] D. Thirumalai and R.D. Mountain, Phys. Rev. E 47, 479 (1993). [11] L.J. Lewis, Phys. Rev. B 44, 4245 (1991).
[12] P.A. Duine, J. Sietsma and A. van den Beukel, Phys. Rev. B 48, 6957 (1993). [13] K. Rätzke, P.W. Hüppe and F. Faupel, Phys. Rev. Lett. 68, 2347 (1992). [14] H. Teichler, Defect and Diffusion Forum Vols. 143-147, 717 (1997).
[15] A. Heesemann, K. Rätzke, F. Faupel, J. Hoffmann and K. Heinemann, Europhys. Lett. 29, 221 (1995).
[16] C. Massobrio, V. Pontikis and G. Ciccotti, Phys. Rev. B. 39, 2640 (1989). [17] P. Lamparter, W. Sperl, S. Steeb and J. Blétry, Z. Naturforsch. 37a, 1223
(1982).
[18] Binary Alloy Phase Diagrams, editor T.B. Massulski, American Soc. for Metals (1986).
[19] N. Lustig, J.S. Lannin and R. Hasegawa, Phys. Rev. B 34, 6725 (1986). [20] Ch. Hausleitner and J. Hafner, Phys. Rev. B 47, 5689 (1993).
[21] J.C. Li, N. Cowlam and F. He, J. Non-Cryst. Solids 112, 101 (1989).
[22] E.W. Iparraguirre, J. Sietsma, B.J. Thijsse and L. Pusztai, Comp. Mat. Sci. 1, 110 (1993).
[23] M.P. Allen and D.J. Tildesley, Computer simulation of liquids, 1st ed. (Oxford University Press, 1989).
[24] T.A. Weber and F.H. Stillinger, Phys. Rev. B 31, 1954 (1985). [25] D. Turnbull and M.H. Cohen, J. Chem. Phys. 29, 1049 (1958). [26] W.A. Phillips, Rep. Prog. Phys. 50, 1657 (1987).
1. 2.
chapter 2
Computer model of amorphous Ni
81
B
19
2.1 Introduction
Experimentally, the study of individual atoms is only possible at surfaces and even then with strong limitations. In computational science it is possible to study the behavior of individual atoms in the bulk (see Fig. 2.1). Molecular dynamics (MD) simulations allow atoms to be followed in time, and processes that take place at an atomic scale can be studied, for example diffusion. Furthermore, it is not only possible to keep track of specific atoms, but it is also possible to study the contribution of individual atoms to collective dynamic phenomena, for instance vibrational modes. The results from an MD simulation allow the calculation of statistical and thermodynamical properties and present the opportunity to study the relation between these properties and the structure at an atomic scale. Of course, there are also some drawbacks: (i) simulations can only span a limited time interval (in the present study 10 ns), (ii) simulations have a limitation in the number of atoms that can be studied (in the present study 2500), and (iii) the model system is only an approximation of the real material, especially in its interatomic interac-tions.
The interaction between atoms is very complex. Ab-initio MD simulations allow the determination of all properties of a material from first principles [1,2]. However, this method is very computer-time consuming and the maximum number of atoms is still very limited (about 100). If one does not want to study electronic properties, a classical model suffices and under conditions that are not too far from thermal equilibrium the atomic interaction may even be effectively represented by a pair potential. The validity of this approach can be assessed from a comparison of the macroscopic properties of the model to the experimentally determined proper-ties [3].
In general, MD studies of amorphous materials are performed in order to investigate structural and dynamical properties of the glassy state. Many authors study simple hard-sphere or Lennard-Jones glasses (see for instance [4-8]). Others try to model a specific amorphous material. The latter group of studies use a wide variety of potentials, in which two types can be distinguished: (i) the theory-based potentials, for instance potentials determined with the “nearly free electron” and “tight binding” approximations [9,10,11], and (ii) the empirical potentials, which are especially constructed to have a model with properties that are in good agreement with available experimental data [12,13, this work]. Both types of potentials use experimental information as input, but the latter are not restricted to a theoretical framework. Overall, the general consensus in the literature is that the application of pair potentials to model glasses is adequate and yields reliable results.
Figure 2.1: Example of an amorphous Ni81B19 molecular dynamics configuration at 304
K. The dark spheres represent the Ni atoms and the light spheres represent the B atoms. The cube has an edge length of 29.4 Å.
Amorphous Ni81B19 model 9 In this chapter the simulation method and the characteristics of the model will be discussed and compared with the experimental data. The experimentally
determined structure of amorphous Ni81B19, in the form of partial reduced radial
distribution functions, is given in section 2.2. The properties of the MD model depend on the algorithm (section 2.3), the pair potentials (section 2.4) and the choice of the starting structure (section 2.5). An overview of the performed
simulations is given in section 2.6, and the properties of the Ni81B19 MD model and
its agreement with available experimental data are discussed in section 2.7.
2.2 Experimental data on the structure
Lamparter et al. have performed a neutron diffraction study on three
isotopically different samples of amorphous Ni81B19 [14]. The samples were
produced by melt-spinning, and the nickel in the alloys was respectively composed of naturalNi, 60Ni and a mixture of 60Ni and 62Ni having a zero coherent scattering length for neutrons. The experiments yielded three structure factors S(q), where q is the length of the diffraction vector. Fourier transformation of S(q), using the expression
( )
(
( )
)
( )
G r = q S q − qr q ∞∫
2 1 0 π sin d , (2.1)leads to the reduced radial distribution function G(r), where r is the interatomic distance. Each of the three functions G(r) is a weighted sum of the partial reduced radial distribution functions Gnn'(r),
( )
( )
G r wnnGnn r n n =∑
' ' , ' , (2.2)where the weighting factor wnn' is given by
w c c b b b nn n n n n ' ' ' = 2 . (2.3)
The weighting factors depend on the coherent scattering lengths for neutrons bn of
atomic species n, the molar fraction cn, and <b> = Σn cnbn.
The structural information contained in Gnn'(r) is best understood by expressing
Gnn'(r) in terms of ρnn'(r), the partial number density of n' type atoms at a distance r
from an n atom,
( )
( )
G r r r c nn nn n ' ' ' = − 4π ρ ρ , (2.4)with ρ the average number density. Using this equation the partial reduced radial
distribution functions can be straightforwardly calculated for an atomic model, thus enabling comparison between experiment and computer model. A related function is the partial pair correlation function gnn'(r), given by
( )
( )
g r r c nn nn n ' ' ' = ρ ρ . (2.5)Although both functions contain the same information, the partial reduced radial distribution functions show more structural details than the partial pair correlation functions, especially for the long range structure, and will be used in the following sections to compare the model structure with the experimental structure.
Figure 2.2 shows the experimental partial reduced radial distribution functions determined by Lamparter et al. [14]. The coherent scattering lengths for the isotopic compositions used are given in Table 2.1. In the paper by Lamparter et al. it is mentioned that one of the total structure factors had to be corrected in a rather ad hoc manner in order to remove physically impossible features; this problem will
be discussed in section 2.5. Note that the Ni81B19 structure displays the short-range
ordering typical for metal-metalloid glasses: even though B is the smallest element
in the system, the first peak of the GBB(r) is present at a larger r-value than the first
peak in the GNiB(r). Therefore, no B atom in the structure is a direct neighbor of
another B atom. Lamparter et al. have also measured the density of naturalNi81B19
using the Archimedian method and obtained ρ = 0.103 Å–3.
Table 2.1: Coherent neutron scattering lengths for Ni and B [14].
n bn / 10 –14 [m] B 0.6 natural Ni 1.03 60 Ni (95.4% enriched) 0.309 62 Ni (98.5% enriched) –0.844
Amorphous Ni81B19 model 11
2.3 Molecular dynamics algorithm
The classical description of a system of atoms is fairly simple. The total force fi
acting on an atom i in the system is the result of the interaction of atom i with the other atoms. The acceleration &&xi of atom i (with mass mi) follows from Newton’s second law,
&&ximi= = −∇ ,fi iV (2.6)
where V is the total potential energy of the system and ∇i is the gradient with
respect to the atomic position vector xi. There are a number of algorithms available
to numerically integrate Newton’s equations of motion. The velocity-Verlet
algorithm [15] performs accurately [16] and has the advantage that the positions xi,
the velocities &xi and the accelerations &&xi of the atoms are known simultaneously 5 − 0 5 0 5 0 5 10 0 5 10 15 r [Å] Gnn ' (r ) [ Å − 2 ] B-B Ni-B Ni-Ni
Figure 2.2: Partial reduced radial distribution functions Gnn'(r) for Ni81B19, as obtained
after each time step. Given the quantities xi(t), &xi(t) and &&xi(t) of all atoms at time
t, the quantities xi(t+∆t), &xi(t+∆t) and &&xi(t+∆t) after a time step ∆t are obtained in the following way,
(
)
( )
( )
( )
(
)
( )
( )
(
)
(
)
(
)
(
)
(
)
x x x x x x x x x x x i i i i i i i i i i i i i t t t t t t t t t t t t t t V t t m t t t t t t t + = + + + = + + = −∇ + + = + + + ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ & &&& & &&
&&
& & && .
1 2 2 1 2 1 2 1 2 1 2 (2.7)
Repeating this procedure h times results in a simulation time span equal to h∆t. In
order to avoid surfaces or interactions with container walls, and in order to allow atoms to move outside the simulation box in the study of diffusion, periodic boundary conditions are used [17]. The efficiency of the calculations is increased by the application of the linked-list method [16].
The simulations are performed at constant temperature T and pressure P. As a result of small relaxations in the model and small computational errors (using a finite time step and a finite precision) unavoidable slow drifts in the temperature and pressure will occur, which require correction.
The temperature is obtained from to the total kinetic energy K by
( )
( )
( )
( )
T T t K t k N k N m t t t i i i N t = = 2 =∑
3 1 3 2 B B &x , (2.8)where kB is the Boltzmann constant, N is the number of atoms, and <...>t denotes
averaging over time. During an isothermal simulation the system is weakly coupled to a “temperature bath”, using a method according to Berendsen [18]: the system is made to obey the equation
( )
( )
d d T T t t T t T = − 0 τ , (2.9)where T0 is the desired temperature and τT is the characteristic coupling time. This
is implemented by scaling the atomic velocities at each time step with a factor
( )
( )
f t t T T t T T = + − 1 0 1 1 2 ∆ τ . (2.10)Amorphous Ni81B19 model 13 The pressure P is determined from the virial theorem,
( )
( )
( )
( )
( )
P P t t k T t N t V t t i i i N t = = ρ B −ρ∑
⋅ ∇ 3 x , (2.11)in which ρ(t) is the atomic number density at time t. The pressure is controlled by a
weak coupling to a “pressure bath” [18]. The box length of the simulation box and the atomic coordinates are scaled at each time step by a factor
( )
[
( )
]
f tP T t P P t P = −1 β 0 − τ ∆ , (2.12)in which P0 is the desired pressure, τP is the coupling time and βT is the isothermal
compressibility.
The temperature and pressure control results in slight perturbations from the Newtonian atomic trajectories. Therefore, fluctuation theory can not be used to determine thermodynamic quantities, but time-average thermodynamic quantities are not disturbed [18].
2.4 Potential energy function
2.4.1 Adiabatic approximation
In any material the atomic nuclei and the electrons form a single quantum system. In more simple terms and distinguishing between core electrons and valence electrons, one could say that ionic motion is coupled to the motion of the valence electrons. This is because the electronic arrangement, and hence the contribution of the valence electrons to the total energy of the solid, depends in detail on the particular arrangement of the ion cores. During ionic motion the electronic wave functions will deform in a way that can be quite difficult to deduce with any precision. One copes with this problem by making the so-called adiabatic approximation. This approximation is based on the fact that electronic velocities typically are much greater than ionic velocities. One therefore assumes that at any moment the electrons will be in their ground state for each particular instantaneous ionic configuration, and the way in which this ground state energy (plus the ion-ion energy) depends on the ionic coordinates defines the potential energy V [19].
2.4.2 Pair potentials
In the pair potential approximation the potential energy of the system V is expressed as
( )
( )
( )
V U V nn rii V i N i N = + = + = =∑
∑
c ρ c ρ 1 2 1 1 Φ ' ' ' , (2.13)where Vc(ρ) is the cohesive contribution, which only depends on the density ρ, U is
the pair potential energy, resulting from the sum over all atoms of the effective pair potential interaction Φnn', n is the type of atom i, and r
ii'
is the distance between the atoms i and i'. As will be shown in section 2.7.1, the density dependence of the cohesive contribution is very small and the density fluctuations resulting from the pressure control can be neglected. Therefore, also the cohesive contribution can be neglected, being a constant contribution, which does not affect the atomic movements. For the effective pair potential the following form is used,
( )
( )
( )
(
)
(
(
)
)
Φnn p q r A ar br ar c r c a r c a ' [ ] exp[ ] , = − − < < ≥ − − −1 0 0 (2.14)in which the parameters A, a, b, c, p and q are positive and depend on the types n and n' of the atoms being considered (Weber and Stillinger [20]). This potential has also been used for amorphous NiP [21] and amorphous NiZr [22]. The potential has the advantages that no cutoff distance is needed (this is important when performing low-temperature simulations) and that the potential and its derivatives have no discontinuities in r. Starting from the nearly-free-electron tight-binding-bond potentials presented by Hausleitner and Hafner [11], modified with the data from Li et al. [13] (see Ref. [23]), we have empirically modified the parameters in order to obtain a system that (i) has partial reduced radial distribution functions that show a very good agreement with the experimental data of Fig. 2.2, (ii) has a
density close to the experimentally measured value of 0.103 Å–3 at room
temperature, and (iii) shows no segregation of B atoms, as observed in earlier work [11,23].
Table 2.2: Parameters of the pair potentials used for amorphous Ni81B19 [Eq. (2.14)]. The
quantities rmin and Φmin refer to the bottom of the pair potential-energy well.
n-n' A [eV] a [Å–1] b [Å–1] c p q rmin [Å] Φmin [eV] Ni-Ni 5.164 0.3757 0.3726 2.540 6.376 5.988 2.72 –0.140 Ni-B 2.586 0.3854 0.3726 2.500 4.818 4.422 2.18 –0.273 B-B 0.765 0.3672 0.3726 2.000 10.51 6.995 3.09 –0.026
Amorphous Ni81B19 model 15
The resulting parameters are given in Table 2.2 and the potentials are shown in
Fig. 2.3. With these, the number density of the system becomes 0.0978 Å–3 at
304 K. Other simulation results are compared with the available experimental data in section 2.7. Table 2.2 also reports the interatomic distances and the potential energies of the minima in the pair potentials. As expected in the view of the experimental radial distribution functions, the potentials have a strongly non-additive character; no B atom can be a direct neighbor of another B atom. (The effective diameter of a B atom is about 1.64 Å, as follows from the minimum positions of the Ni-Ni and Ni-B potentials).
2.5 Starting structure
The usual method of obtaining an amorphous material is by quenching a liquid to a low temperature. Experimentally, the maximum achievable cooling rate is in
the order of 106 Ks–1. Computationally, the minimum quench rate is in the order of
1012 Ks–1, on account of the limited available time span of a simulation. From
experiments it is known that the properties of a glass depend on the cooling rate
−0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 0 1 2 3 4 5 6 r [Å] Φnn ' (r ) [eV] Ni-Ni Ni-B B-B
(see for example Ref. [24,25]). Therefore, the large difference between the experimental and the computational cooling rates may result in differences in the glassy properties.
Using an alternative approach we have obtained the starting configuration by applying the reverse Monte Carlo (RMC) method. This is a technique by which the structure of matter, usually disordered matter, can be obtained at the atomic scale [26,27]. RMC generates atomic configurations that correspond to available experimental structural information as closely as possible. Typically this informa-tion is a set of (partial reduced) radial distribuinforma-tion funcinforma-tions or structure factors, obtained from diffraction experiments. In general, the calculations are carried out by iteratively moving a randomly picked atom over a randomly chosen distance and comparing the resulting structural information with the experimental data. If the model shows a better agreement than before, the move is accepted; if not, the move is accepted with a certain probability (which prevents the configuration to remain trapped in a local minimum). The procedure is repeated until the measure of the
goodness-of-fit reaches a maximum, which typically requires 106 atomic moves.
Early RMC calculations on Ni81B19 were performed by Iparraguirre et al. [28],
using the experimental as-quenched partial reduced radial distribution functions
Gnn'(r) of Lamparter et al. [14], the experimental density ρ = 0.103 Å
–3
, and 1920 atoms. In order to account for the uncertainty in the normalization of the measured diffraction intensities (see also section 2.2) the renormalized data
( )
( )
( )
Γnn nn nn r r G r G r ' ' ' = , (2.15)were used as input for the RMC-program, instead of Gnn'(r) itself. Here <|Gnn'(r)|>r
denotes the average over r of the absolute value of Gnn'(r). This approach allows
for a possible multiplicative error in the normalization of the original diffraction intensities.
Table 2.3: Renormalization factors for RMC calculations on the experimental Ni81B19
partial reduced radial distribution functions of Lamparter et al. [14].
ρ [Å–3
] γNiNi γNiB γBB
Iparraguirre et al. [28] 0.103 0.018 0.377 0.077
Amorphous Ni81B19 model 17
Denoting the renormalization factors by 1 + γnn', the values for γnn' obtained by
Iparraguirre et al. are given in Table 2.3. The effect of the normalization errors proved to be considerable, as the renormalization factor for Ni-B distribution function is quite high.
The simulations presented in this work are performed with a different number of
atoms (2500 instead of 1920 atoms) and at a different number density (0.0978 Å–3
instead of 0.103 Å–3). Therefore, we did not use Iparraguirre’s configuration, but
new RMC calculations were performed to obtain the starting configuration for the MD simulations. The RMC method was applied with (i) an equilibrium liquid configuration (1214 K) at high pressure, in order to reach the desired number density, and (ii) the partial reduced radial distribution functions of Lamparter et al. [14], renormalized according to Eq. (2.15), as input data. The RMC method does not use potential energies. Therefore, in order to fine-tune the resulting structure to
5 − 0 5 0 5 0 5 10 0 5 10 15 r [Å] Gnn ' (r ) [ Å − 2] B-B Ni-B Ni-Ni 5 − 0 5 0 5 0 5 10 0 5 10 15 r [Å] Gnn ' (r ) [ Å − 2] B-B Ni-B Ni-Ni
Figure 2.4: Partial reduced radial
distribu-tion funcdistribu-tions Gnn'(r) after the application
of the RMC-method (solid curves), com-pared with the experimental data of Lam-parter (dashed curves) [14] corrected for normalization errors (Table 2.3).
Figure 2.5: Partial reduced radial
distribu-tion funcdistribu-tions Gnn'(r) of two simulations at
304 K. The first simulation started from the RMC configuration (solid curves) and the second simulation started from a quen-ched configuration (dashed curves). The difference in statistics between the two sets of data results from taking the average over a different number of configurations.
the pair potentials, the RMC run was interrupted two times for a short (70 ps) MD
simulation at 304 K. The new renormalization factors γnn' are given in Table 2.3.
The change in density decreases the Ni-B normalization error but increases the B-B normalization error. The resulting partial reduced radial distribution functions are in excellent agreement with the renormalized experimental data, as can be seen in Fig. 2.4.
For comparison a configuration at 1417 K (which in section 2.7.1 will be shown to be in the liquid phase) has been cooled to 304 K with a cooling rate of
1012 Ks–1. The partial reduced radial distribution functions of the resulting
amorphous structure are shown in Fig. 2.5, together with those of a simulation at 304 K that started from the RMC configuration. The differences are very small and
are only visible in GBB(r). Both simulations have equal (within the error) number
densities and equal average potential energies per atom. In early stages of this work, pair potentials were used that were not yet specially adapted to the structure. This resulted in much larger differences, and it was shown that the RMC method yielded a more realistic starting structure than quenching a liquid [23]. For the present case, the differences between the structures obtained by RMC and by quenching a liquid are very small, and it is not possible to determine beforehand which configuration would yield the most realistic result. It was decided to use the RMC configuration as the starting point for all simulations.
2.6 Simulation details
2.6.1 Constants
The amorphous Ni81B19 system consists of N = 2500 atoms (2025 Ni and 475 B
atoms) in a cubic box with periodic boundaries. The equations of motion [Eq.
(2.7)] are solved using a time-step of ∆t = 1.77 fs. All simulations are performed at
constant temperature and ambient pressure. The temperature control [Eq. (2.10)]
with τT = 8.86 fs keeps the standard deviation in the temperature always below
0.6% of the selected temperature. The pressure control [Eq. (2.12)] with βT / τP =
976 Pa–1s–1 results in a standard deviation of the density at 304 K and at 1417 K of
Amorphous Ni81B19 model 19
2.6.2 Temperatures and time scales
In total fourteen simulations are performed at different temperatures. Table 2.4 lists the temperatures and the simulation time spans. Each simulation is started from the RMC configuration. The time span encompassed by a simulation is given
by the total simulation time ts. With decreasing temperature the dynamics of the
system slow down and the time necessary to equilibrate the model increases. This
can be seen in Fig. 2.6, which shows the mean square drift <|r(t)|2>N versus time
for all investigated temperatures, where ri(t) is the drift xi(t) – xi(0) of atom i after
time t, and <...>N denotes the average over all atoms. Note that the mean square
drift is shown for the total simulation time ts. Some of the mean square drift curves
show sudden temporary increases in the slope (a very pronounced one can be seen in the simulation at T = 354 K after t = 2 ns). These periods of increased atomic mobility are thought to be relaxations that take place only when the system has not yet reached equilibrium. For high temperatures these changes are not observed and
these models are already in an equilibrium state after t ≈ 10 ps. With decreasing
temperature the events of increased mobility occur after larger time intervals and the time needed to reach an equilibrium state therefore increases.
Table 2.4: Overview of the simulations: the temperature T, the total simulation time ts,
the equilibration time te, the production time tp, and the number of
configura-tions studied in the production time are given.
T [K] ts [ns] te [ns] tp [ns] number of config. 304 10.0 7.2 2.8 401 354 9.5 5.7 2.8 401 405 9.0 6.2 2.8 401 455 8.5 5.7 2.8 401 506 8.0 5.5 2.5 351 557 5.5 3.4 2.1 301 607 3.0 1.2 1.8 251 658 2.0 0.6 1.4 201 708 1.5 0.4 1.1 151 809 1.0 0.3 0.7 101 911 1.0 0.3 0.7 101 1012 1.0 0.3 0.7 101 1214 1.0 0.3 0.7 101 1417 1.0 0.3 0.7 101
The pronounced change in the model at 354 K (t ≈ 2 ns) causes a small drop (0.04%) in the total potential energy and an increase in the density of about 0.05%.
Relaxations of this kind are also observed by Hörner in a model of Fe90Zr10 [9],
and others [5,8] also report sudden changes in the slope of the mean square drift versus time. Comparison with literature is difficult as many authors only report that the curve of the mean square drift in time contains much noise but do not show these curves. Furthermore, the total simulation times of these simulations are reported to be on the order of 1 ns or less, and, as can be seen in Fig. 2.6, the time intervals between changes usually are considerably larger.
With the exception of the simulation at 354 K, the total simulation time for all
simulations is given by the sum of the equilibration time te, the time needed to reach
equilibrium, and the production time tp, the time used for studying the model. The
model at 354 K did not reach the equilibrium state in the given simulation time. The mean square drift of the model shows a sudden increase at the end of the
10
−110
010
110
210
310
410
−210
−110
010
1t [ns]
<|
r(
t)|
2>
N[
Å
2]
304 K 354 K 405 K 455 K 506 K 557 K 607 K 658 K 708 K 809 K 911 K 1012 K 1214 K 1417 KFigure 2.6: Mean square drift <|r(t)|2>N versus time t, over the total simulation time, at
Amorphous Ni81B19 model 21 simulation time (8.5 - 9.5 ns). The data produced during this period is not used and the production time for this model starts at t = 5.7 ns and ends at t = 8.5 ns. Most likely the model at 304 K has also not reached the equilibrium state, but the mean square displacement does not reveal any relaxations between 7.2 - 10 ns, and this period is used for study.
All time averages presented in this thesis are obtained in the production time, as well as most of the data of events at specific points in time. During the production time configurations are saved every 7.1 ps; the number of configurations that are studied is given in Table 2.4.
2.6.3 Equilibrium configurations
At temperatures above T = 0 K all atoms vibrate in their potential wells. It will be shown in the following chapters that, besides studying an atomic configuration at a certain temperature, it is also useful to study the atomic equilibrium positions, being the positions of the atoms when they are located at the bottom of their potential wells. These positions are obtained by an energy minimization procedure. Starting with a system at temperature T, the set point for the system thermostat is
set at T0 = 50 K and all velocities are set to zero. In a simulation of 3.5 ps the
thermostat is then programmed to cool the system to T = 0 K with dT/dt =
–1.43×1013 K s–1. The resulting temperature initially rises above 50 K and then
decreases towards 0 K. After 3.5 ps the system has a temperature on the order of
10–1 K. In a subsequent simulation of 3.5 ps the temperature control (set to 0 K)
proceeds to cool the system, resulting in a final temperature of about 10–5 K.
2.7 Model characteristics
2.7.1 Glass transition
The model system undergoes a glass transition at Tg ≈ 590 K, as can be clearly
seen in Fig. 2.7, where the time averaged number density <ρ>t, the
Wendt-Abraham ratios WNiNi and WNiB, and the time-averaged pair potential energy per
atom <u>N,t are plotted versus temperature. The Wendt-Abraham ratio
Wnn'=
( )
(
)
g r g r nn nn ' ' min max 1 1 , (2.16)is the ratio of first minimum, at r = rmin1, and the first maximum, at r = rmax1, of the
pair correlation function gnn'(r) [Eq. (2.5)]. The Wendt-Abraham ratio is a measure
for the degree of order in the nearest neighbor shell.
The glass transition temperature of amorphous Ni81B19 is experimentally not
known. However, the value of 590 K obtained from the presented simulations is within the experimentally observed range of glass temperatures for Ni-based glasses. The cubic thermal expansion coefficient is a measure of the anharmonicity
1.6 − 1.5 − 1.4 − 1.3 − 0 0.1 0.2 0.080 0.085 0.090 0.095 0.100 0 500 1000 1500 T [K] < u >N ,t [eV] Wnn ' < ρ >t [ Å − 3 ] Ni-Ni Ni-B Tg
Figure 2.7: Glass transition at Tg = 590 K, as evidenced by a change in slope of the
time-averaged number density <ρ>t (circles), the Wendt-Abraham ratios WNiNi
(diamonds) and WNiB (triangles), and the time-averaged pair potential energy
per atom <u>N,t (squares) as a function of temperature. The lines are
Amorphous Ni81B19 model 23
of the used potentials and has a reasonable value of 8.4×10–5 K–1 below Tg and
16.4×10–5 K–1 above Tg. Omitting the cohesive contribution to the total potential
energy, the specific heat below the glass transition temperature is found to be 25.8
J mol–1K–1, only slightly larger than the classical value 3R (24.9 J mol–1K–1). This
indicates that the density dependence of the cohesive contribution is small.
5
−
0
5
0
5
0
5
10
0
5
10
15
r [Å]
G
nn '(r
) [
Å
− 2]
B-B Ni-B Ni-NiFigure 2.8: Partial reduced radial distribution functions Gnn'(r) of the MD-system at 304
K (solid curves), compared with the experimental data of Lamparter (dashed curves) [14] corrected for normalization errors (Table 2.3).
2.7.2 Radial distribution functions
The configurational average over the production time of the partial reduced
radial distribution functions GNiNi(r), GNiB(r) and GBB(r) of the model system at T =
304 K are shown in Fig. 2.8 together with the experimental data from Lamparter et al. [14]. Only slight differences occur, mainly in the height of the first peaks and in
the splitting of the second peaks. The distribution function GBB(r) for the
5 − 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 0 5 10 15 r [Å] GNiNi (r ) [ Å − 2 ] 304 K 354 K 405 K 455 K 506 K 557 K 607 K 658 K 708 K 809 K 911 K 1012 K 1214 K 1417 K 5 − 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 0 5 10 15 r [Å] GNiB (r ) [ Å − 2 ] 304 K 354 K 405 K 455 K 506 K 557 K 607 K 658 K 708 K 809 K 911 K 1012 K 1214 K 1417 K 5 − 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 0 5 10 15 r [Å] GBB (r ) [ Å − 2 ] 304 K 354 K 405 K 455 K 506 K 557 K 607 K 658 K 708 K 809 K 911 K 1012 K 1214 K 1417 K
Figure 2.9: Partial reduced radial
distribu-tion funcdistribu-tions Gnn'(r) for all
Amorphous Ni81B19 model 25 system shows a small shift to larger r values in the medium and long range order. This is probably caused by the small well depth of the B-B potential (Table 2.2), which was necessary to avoid the segregation of B atoms. Despite these small differences, the overall agreement is very good and corroborates the validity of the constructed pair potentials.
Figure 2.9 shows the partial reduced radial distribution functions for all investigated temperatures. With increasing temperature the order in the system decreases. The splitting in the second peak, which in the past was commonly
thought to be a feature of the glassy structure, does not vanish at Tg, but disappears
approximately 200 K higher in temperature. This was also reported by Lewis [21] and Barrat [29].
2.7.3 Vibrational density of states
The partial vibrational densities of states Jn(ω) for atoms of type n can be
calculated by taking the Fourier transform of the corresponding normalized velocity autocorrelation function Z [30], given by
( ) ( )
( )
Z t i i N i N = & ⋅& & x x x 0 0 2 . (2.17)From the two distributions JNi(ω) and JB(ω) the neutron-weighted total vibrational
density of states J(ω) can be determined, according to
( )
( )
J c b J m n n n n n ω ∝∑
2 ω , (2.18)where the sum is taken over all elements n, mn is the atomic mass, bn is the
coherent scattering length for neutrons (given in Table 2.1), and cn is the molar
concentration [31,32]. In Fig. 2.10 the normalized total and partial density of states
for the model system at T = 304 K are shown together with the experimental data
of Lustig et al. [33]. It can be seen that the B atoms dominate the high frequency part of the spectrum and the Ni atoms perform the slow vibrations. Comparing the total vibrational density of states with the experimental result, we observe that the shoulder on the main peak is more pronounced for the simulation than for the experiment, causing the first observed maximum to be lower and the second to be higher. The reason for this is unknown. Nevertheless, both maxima are situated at the correct frequencies and the overall agreement is satisfactory.
2.8 Conclusions
Molecular dynamics simulations are used to model amorphous Ni81B19. The
starting structure for the simulations is determined by the application of the reverse Monte Carlo method using the experimentally determined partial reduced radial distribution functions. The model consists of 2500 atoms in a cubic box, using periodic boundary conditions. The used pair potentials are of a Weber-Stillinger type with empirically modified parameters to have a model that agrees closely to the available experimental data. Simulations are performed at fourteen constant temperatures and at ambient pressure. The model displays a glass transition at 590 K. The partial reduced radial distribution functions and the vibrational density of states are in good agreement with the experimental data.
0 1 2 3 0 50 100 150 ω [ps−1] J( ω ) [arb. units] JNi JB J
Figure 2.10: Neutron-weighted vibrational density of states. The squares give the experimental data from Lustig et al. [33]. The partial vibrational densities
of states JNi(ω) and JB(ω) are determined from the velocity autocorrelation
functions at 304 K. The total vibrational density of states J(ω) follows from
the partial densities [Eq. (2.18)]. All curves are normalized to the same area.
Amorphous Ni81B19 model 27
References
[1] G. Galli and M. Parrinello, Computer Simulation in Materials Science, 283
(Kluwer, Netherlands, 1991).
[2] J. Hafner, Statics and Dynamics of Alloy Phase Transformations, Edited by
P.E.A. Turchi and A. Gonnis, 269 (Plenum Press, New York, 1994).
[3] V. Vitek, Mat. Res. Soc. Bulletin 21 (2), 20 (February 1996).
[4] D. Thirumalai and R.D. Mountain, Phys. Rev. E 47, 479 (1993).
[5] Y. Hiwatari, H. Miyagawa and T. Odagaki, Sol. Stat. Ion. 47, 179 (1991).
[6] J.M. Delaye and Y. Limoge, J. Non-Cryst. Solids 156, 982 (1993).
[7] G. Wahnström, Phys. Rev. A 44, 3752 (1991).
[8] J. Ullo and S. Yip, Phys. Rev. A 39, 5877 (1989).
[9] A. Hörner, Ph.D. thesis, Fakultät Physik der Universität Stuttgart, Germany,
(1993).
[10] G.W. Tecza and J. Hafner, J. Non-Cryst. Solids 202, 1 (1996).
[11] Ch. Hausleitner and J. Hafner, Phys. Rev. B 47, 5689 (1993); private communication.
[12] U. Buchenau, Yu.M. Galperin, V.L. Gurevich, and H.R. Schober, Phys. Rev. B 43, 5039 (1991).
[13] J.C. Li, N. Cowlam and F. He, J. Non-Cryst. Solids 112, 101 (1989).
[14] P. Lamparter, W. Sperl, S. Steeb and J. Blétry, Z. Naturforsch. 37a, 1223 (1982).
[15] W.C. Swope, H.C. Andersen, P.H. Berens, and K.R. Wilson, J. Chem. Phys. 76, 637 (1982).
[16] M.P. Allen and D.J. Tildesley, Computer simulation of liquids, 1st ed. (Oxford University Press, 1989).
[17] M. Born and Th. Von Karman, Physik. Z. 13, 297 (1912).
[18] H.J.C. Berendsen, J.P.M. Postma, W.F. Van Gunsteren, A. Di Nola, and J.R. Haak, J. Chem. Phys. 81, 3684 (1984).
[19] N.W. Ashcroft and N.D. Mermin, Solid State Physics (Saunders College, Philadelphia, 1976).
[20] T.A. Weber and F.H. Stillinger, Phys. Rev. B 31, 1954 (1985). [21] L.J. Lewis, Phys. Rev. B 44, 4245 (1991).
[22] H. Teichler, Phys. Stat. Sol. b 172, 325 (1992).
[23] B.J. Thijsse, L.D. van Ee and J. Sietsma, Mat. Res. Soc. Symp. Proc. 321, 65 (1994).
[24] H.S. Chen, J.T. Krause, A. Inoue, and T. Masumoto, Scripta Metall. 17, 1413 (1983).
[26] R.L. McGreevy and L. Pusztai, Mol. Sim. 1, 359 (1988). [27] L. Pusztai, Z. Naturforsch. 46a, 69 (1991).
[28] E.W. Iparraguirre, J. Sietsma, B.J. Thijsse and L. Pusztai, Comp. Mat. Sci. 1, 110 (1993).
[29] J.L. Barrat and M.L. Klein, Annu. Rev. Phys. Chem. 42, 23 (1991).
[30] J.P. Hansen and I.R. McDonald, Theory of simple liquids (Acad. Press, 1986).
[31] G.E. Bacon, Neutron Diffraction, 3rd ed. (Clarendon Press, Oxford, 1975). [32] N. Lustig, J.S. Lannin, J.M. Carpenter and R. Hasegawa, Phys. Rev. B 32,
2778 (1985).
1. 2. 3. kk
chapter 3
Atomic motion
3.1 Introduction
The most straightforward approach to study of the dynamics of atoms is simply to follow their motion. However, just generating atomic trajectories is not sufficient. One has to find a way to identify and distinguish certain types of movement. In dense solids there are conceptually two types of atomic motion, and the displacement of the atoms can be separated into two contributions: (i) the vibrational movement of an atom around its potential-energy minimum, and (ii) the occasional jump of an atom (or a group of atoms) from its potential-energy minimum into a neighboring potential-energy minimum. The first type of movement takes place in the “cage” formed by the neighboring atoms, and is therefore called the in-cage movement. This oscillatory movement is the general type of movement performed by all atoms, and the amplitudes of the vibrations are small compared to the interatomic distances. The second type of movement is called diffusion. In liquids, and to a lesser extent in glasses, the local potential-energy minima may change in time and are not fixed in space. Therefore, in the liquid state the two types of motion are strongly coupled and can not be separated, whereas in a crystalline solid the two types of motion can be separated. The glassy state lies somewhere in between, and we can expect to observe both types of motion in a combined manner and possibly even more complicated types. In fact we shall find that the diffusional atomic motion in an amorphous solid is not limited to jumps of single atoms (like vacancy jumps in a crystal), but takes place by a cooperative process, i.e. a number of atoms moving simultaneously to a new potential-energy minimum.
The mean square drift of an atom, which is the most straightforward quantity in the analysis of atomic dynamics, is discussed in section 3.2, where special effort is put in separating the in-cage movement from the diffusion. From the mean square drift the diffusion coefficient can be obtained. In section 3.3 the distributions of atomic displacements are presented in the form of Van Hove probability functions. These functions have a theoretically predicted Gaussian shape in two limiting cases: for purely harmonic in-cage movement (atoms vibrate in local potential wells) and for purely uncorrelated diffusional motion (no correlation of the position of the atoms with their initial positions). The transition from pure in-cage movement to pure diffusion is discussed in section 3.4. A special type of collective atomic motion, in the form of atoms jumping as a chain, is the topic of section 3.5. A more general investigation of collective atomic motion will be presented in the chapters 6 and 7.
3.2 Mean square drift
3.2.1 In-cage movement
If at time t0 an atom i is located at position xi(t0), the displacement of this atom
after time t is defined as
ri(t) = xi(t0 + t) – xi(t0). (3.1)
The mean square drift <|r(t)|2>Nn is defined as the square displacement averaged
over all atoms Nn of type n. At time t0 the atoms have certain displacements from
their equilibrium positions, which results in a non-essential component in the mean
square drift. We have suppressed this influence of the starting positions xi(t0) by
calculating the average mean square drift <|r(t)|2>Nn,t0; all configurations in the
production time are used as starting points (t0), and the resulting mean square drifts
are averaged.
In the case of a purely harmonic potential, the in-cage contribution to the mean square drift is proportional to the temperature. To verify this, the average mean square drift divided by the temperature is shown for short simulation times in Fig. 3.1 for all investigated temperatures. Initially the only contribution to <|r(t)|2>Nn,t0
comes from the in-cage movement. The curves for the simulations performed below
the glass-transition temperature (Tg = 590 K) do not diverge and reach a constant
level after 0.2 ps. The slight temperature dependence of the stationary level is caused by the anharmonic part of the pair potentials and the rare occurrence of diffusion events. However, the effect is only small and the in-cage contribution can
be taken as a constant Cn. For Ni atoms this constant has a value CNi ≈ 2.4×10
Atomic motion 31
Å2K–1 and for B atoms CB ≈ 3.8×10– 4 Å2K–1. Above the glass-transition
tempera-ture the divergence of the curves becomes more pronounced. Again the anharmoni-city of the potentials plays a role but the largest part of the increment is the result of the diffusion process. Both the in-cage contribution and the diffusion contribution to the mean square drift of the B atoms are larger than for the Ni atoms. This indicates that the B atoms have a smaller curvature of the local potential well and have a higher mobility than the Ni atoms.
3.2.2 Diffusion
For isotropic stationary systems, the relation between the self-diffusion
coeffi-cient Dn, with n the atomic type, and the mean square drift is defined by the
Einstein equation for diffusion [1,2], given by
D t t n t Nn = →∞ 1 6 2 lim ( ) d d r . (3.2)
In addition to the Einstein equation, the self-diffusion coefficient is also given by 0 0.0005 0.0010 0 0.2 0.4 0.6 t [ps] <| r( t)| 2 > Nn ,t0 T − 1 [ Å 2 K − 1 ] Ni 1417 K 304 K 0 0.0005 0.0010 0 0.2 0.4 0.6 t [ps] <| r( t)| 2 > Nn ,t0 T − 1 [ Å 2 K − 1 ] B 304 K 1417 K
Figure 3.1: Average mean square drift <|r(t)|2>Nn,t0 divided by temperature T versus time
t for Ni and B atoms. Different curves correspond to different temperatures.
the Kubo equation [3], defined as
(
) ( )
Dn t t t N t n = 1∞∫
+ ⋅ 30 x& 0 x& 0 d . (3.3)For small Dn, the use of the Kubo equation becomes hazardous, because Dn
depends entirely on a cancellation between the positive and the negative parts of the velocity autocorrelation function [Eq. (2.17)], and it consequently requires very long integration times t. The Einstein equation, which involves only the positive contributions to the mean square drift is therefore preferred.
Besides a diffusional contribution, the mean square drift also contains a contribution of the in-cage movement and is influenced by the arbitrary atomic
starting positions xi(t0). In Eq. (3.2) the limit for long times is taken, because both
these effects are time independent. Using a special approach, we determine the mean square drift that only contains information about the diffusion process: (i) we suppress the influence of the starting positions, in a way similar to the procedure in section 3.2.1, by calculating the average mean square drift <|r(t)|2>Nn,t0; (ii) we
determine the mean square drift for equilibrium configurations. In an equilibrium
configuration each atom is located at its equilibrium position xieq(t0 + t), the bottom
of its potential well. The equilibrium configurations are created by applying the energy minimization procedure described section 2.6.3.
10−3 10−2 10−1 100 101 102 103 104 101 102 103 104 t [ps] <| req (t )| 2 > Nn ,t0 [ Å 2 ] Ni 304 K 354 K 405 K 455 K 506 K 557 K 607 K 658 K 708 K 809 K 911 K 1012 K 1214 K 1417 K 10−2 10−1 100 101 102 103 104 101 102 103 104 t [ps] <| req (t )| 2 > Nn ,t0 [ Å 2 ] B 304 K 354 K 405 K 455 K 506 K 557 K 607 K 658 K 708 K 809 K 911 K 1012 K 1214 K 1417 K
Figure 3.2: Average equilibrium mean square drift <|req(t)|2>Nn,t0 versus time t, at
different temperatures, for Ni and B atoms. The heavy line indicates a unit slope.
Atomic motion 33
Figure 3.2 shows the average equilibrium mean square drift <|req(t)|
2
>Nn,t0 versus
time for Ni and B atoms, determined over the production time for all investigated temperatures (see section 2.6.2). In this figure the data is presented on a double log scale. If Eq. (3.2) applies, the mean square drift has a unit slope. For temperatures
above the glass-transition temperature (Tg = 590 K) we find regular diffusive
behavior for long times. At temperatures around and below Tg, the mean square
drift reaches a unit slope only after a period of time. This period increases with decreasing temperature and for the two lowest temperatures (304 K and 354 K) a unit slope is not observed in the investigated time interval. The atomic mobility at these two temperatures is very low and in the investigated time interval the diffusion process has not yet reached a state that is uncorrelated with the atomic
positions at time t0 (this topic will be discussed in more detail in section 3.4).
For the simulations performed around and below the glass-transition tempera-ture one observes, on the basis of a back-extrapolation of the observed unit slope for longer times, that the first part of the mean square drift (t < 100 ps) is higher than expected. In addition, when the initial mean square drift (t = 7.1 ps) is compared for different temperatures, it turns out to contain a large scatter in the low-temperature range. For instance, the initial mean square drift for Ni atoms at 354 K is even larger than the initial mean square drift at 405 K. Apparently, there are atomic displacements taking place that have a relatively larger effect on the drift for short times than for long times. In section 3.5 and chapter 6 and 7 it will be shown that the model contains reversible two-level states, sites at which the system can jump from its potential-energy minimum to a neighboring minimum and back. These reversible two-level states can explain the disproportionally high initial mean square drift. This can be understood as follows: In the calculation of the mean square drift for all available short time intervals, both the forward and backward jumps result in diffusion events observed in many intervals. On the other hand, in the calculation of the mean square drift for all available long time intervals, the atomic motion at the two-level state will not be taken into account in an increasing number of intervals when both the forward and the backward jumps take place within the interval. Therefore, the occurrence of reversible two-level states will result in an increased initial mean square drift as compared to the drift for longer times.
From the equilibrium mean square drifts the self-diffusion coefficients for Ni and B atoms at all investigated temperatures are obtained by a least-squares fit of Eq. (3.2) to the data points. The limit for long times is still applied in order to exclude the previously discussed short-time effect in the mean square drift. In spite
