**2. Potential energy function**

According to the embedded atom method, the cohesive energy of an assembly of *N* atoms is given by [16, 17]

$$E\_{htt} = \sum\_{i} F\_i(\rho\_i) + \sum\_{i>j} \phi(r\_{ij}) \tag{1}$$

$$\rho\_i = \sum\_{j(\neq i)} f(r\_{ij}) \,. \tag{2}$$

where *Etot* is the total cohesive energy, ρ*i* is the host electron density at the location of atom *i* due to all other atoms, *f*(*rij*) is the electronic density function of an atom, *rij* is the distance

these approaches: the calculation of the free energy between FCC and HCP structures [6, 7], the investigation of first order phase transition [8], the dependence of the phase diagram on the range of attractive intermolecular forces [9], the investigation of harmonic lattice dynamics and entropy calculations in metal and alloys [10], the calculation of the *P*-*T* diagram of hafnium [11], etc. Recently, the *P-T* diagrams for Ni and Al have been calculated by Gurler and Ozgen [12] by using the MD simulations based on the EAM technique [13]. The reliability of the results obtained from MD simulations depends on the suitable modeling of the interatomic interactions. Interatomic interactions are usually results of fits to various experimental data. However, it is not clear whether simulations performed at other temperatures still reproduce the experimental data accurately. Comparing theoretical and experimental elastic constants and other properties at various temperatures can serve as a measure of reliability and usefulness of potential models [14, 15]. In fact, there are several potential energy functions that can be used for the metallic systems. However, the EAM, originally developed by Daw and Baskes [16, 17] to model the interatomic interactions of face-centered cubic (FCC) metals, has been successfully used to compute the properties of metallic systems such as bulk, surface and interface problems. The reliability of the EAM in

the bulk and its simple form for use in computer simulations make it attractive.

dynamics is applied.

**2. Potential energy function** 

() () *tot i i ij*

given by [16, 17]

When a liquid metal is quenched through the super-cooled region, a phase transition from liquid to glass takes place. Several techniques have been proposed to obtain a disordered state [18-20]. Among them the rapid solidification method is widely used for the amorphous phase. However, due to the demand of a high cooling rate this method is restiricted in most experimental cases. Thus, the computer simulation of molecular

In this study, in order to model Au metallic systems we have used the EAM functions modified by us (Ciftci and Colakoğlu [21]), developed firstly by Cai [22]. In this work, we have carried out MD simulations to obtain the *P*-*V* diagrams at 300 K and the *P*-*T* diagrams of the systems for an ideal FCC lattice with 1372 atoms, by using an anisotropic MD scheme. In addition, the bulk modulus and specific heat of the system in solid phase are determined and results-driven simulations are interpreted by comparing with the values in literature. We have also calculated the pressure derivatives of elastic constants and bulk moduli for Au. The obtained results are compared with the values in the literature. The another purpose of this work is to explore the glass transition and crystallization of Au using EAM .

According to the embedded atom method, the cohesive energy of an assembly of *N* atoms is

*i ij EF r* 

> ( ) ( ) *i ij j i*

where *Etot* is the total cohesive energy, ρ*i* is the host electron density at the location of atom *i* due to all other atoms, *f*(*rij*) is the electronic density function of an atom, *rij* is the distance

 *f r* 

 

(1)

, (2)

between *i* and *j* atoms, *Fi*(ρ*i*) is the embedding energy to embed atom *i* in an electron density ρ*i*, and (*rij*) is the pairwise potential energy function between atoms *i* and *j*.

In this work, we used a modified pairwise potential function in the framework of the Cai version [22] of the EAM. Recently, this potential function has been used by us for predicting several physical properties of some transitional metals [21,23-25]. The present form of the potential makes it more flexible owing to the constants, *m* and *n* in the multiplier forms. Such a factor included in the classical Morse function is treated by Verma and Rathore [26] to compute the phonon frequencies of Th, based on the central pair potential model. The modified parts of the potential and the other terms are as follows:

$$f(r) = f\_c \mathcal{C}^{-a(r - r\_c)},\tag{3}$$

$$F(\rho) = -F\_0 \left[ 1 - \ln \left( \frac{\rho}{\rho\_\varepsilon} \right)^\iota \right] \left( \frac{\rho}{\rho\_\varepsilon} \right)^\iota + D\_2 \left( \frac{\rho}{\rho\_\varepsilon} \right) \tag{4}$$

$$\phi(r) = \frac{D\_1}{(m-1)} \left[ \frac{e^{-m\beta(\frac{r}{r\_c}-1)}}{(\beta \frac{r}{r\_c})^n} - (\beta \frac{r}{r\_c})^n e^{-\beta(\frac{r}{r\_c}-1)} \right],\tag{5}$$

where *α*, *β*, *D1* and *D2* are fitting parameters that are determined by the lattice parameter *a*0, the cohesive energy *E*c, the vacancy formation energy *E*vf, the elastic constants *Cij.* Here ρ*e* is the host electron density at equilibrium state, *re* is the nearest neighbor equilibrium distance, and *F*0=*E*c*E*vf . In this potential model, there are four parameters: *β* and *D*1 are from twobody term, *m* and *n* are adjustable selected constants. The fitting parameters are determined by minimizing the value of exp exp <sup>2</sup> [( ) / ] *W XX X cal* . Here *X* represents the calculated and experimental values of the quantities taken into account in the fitting process. Hence, the potential functions can be fitted very well to the experimental properties of the matter, such as the vacancy formation energy (*E*v), cohesive energy (*E*c), elastic constants (*C*ij), and lattice constants (*a*0) in an equilibrium state. In the fitting process here, the cutoff distance is taken to be *r*cut=1.65*a*0. In the Eq. (3), the *f*e parameter is selected as unity for mono atomic systems because it is used for alloy modeling as an adjustable parameter to constitute suitable electron density. For the selected values of the constants *m* and *n*, the computed potential parameters and experimental input data for Au are given in Table 1.

The cohesive energy changes with the variation of lattice constants of Au calculated from Eq. (1) and from the general expression of the cohesive energy of metals proposed by Rose et al. [32] are compared in Fig.1. The Rose energy is also called as the generalized equation of state of metals and written as

$$E\_R(a^\*) = -E\_0(1+a^\*)e^{-a^\*} \tag{6}$$

$$a^\* = \left(\frac{a}{a\_0} - 1\right) / \left(\frac{E\_\odot}{9B\_m\Omega}\right)^{1/2} \tag{7}$$

A Molecular Dynamics Study on Au 205

where *mi* is mass of particle *i*, **s***i* is the scaled coordinate of atom *i* and is represented by a column vector whose elements are between zero and unity, **h**=(**a**, **b**, **c**); **a**, **b** and **c** vectors are

which represents mass of the computational box, *P*ext is external pressure applied on the cell, *V* is the volume of the MD cell and is obtained from det(**h**). Thus, square of distance

> 1 1 *ii i mi*

> > 1 ( ) *M Pext*

1

restarted with different pressure in each run, to avoid algorithmic errors.

following the procedure given by Karimi et al [14].

*V m*

Also the force on an atom *i* in the system is calculated from the following equation,

1

*j i E FF*

 

*i si i j j i ij*

where the primes denote the first derivatives of the functions with respect to their

In all of the simulation studies, the equation of motion given in Eqs. (9) and (10) were numerically solved by using the velocity version of the Verlet algorithm [35]. The size of integration step was chosen to be 7.87x1015s for Au. Initial structures of the systems were constructed on a lattice with 1372 atoms and an FCC unit cell. It has been observed that, with these initial conditions, the systems were equilibrated in 5000 integration steps. Time averages of the thermodynamic properties of the system in each simulation run were determined by using 30,000 integration steps following the equilibration of the system. The structures of the system in solid phase were examined by using the radial distribution function. Melting temperatures were determined from the plots of the cohesive energy versus temperature. It is possible to classify our simulation runs in two groups as thermal and pressure applications. In the thermal applications, the temperature of the system under zero pressure is raised from 100K to 2400K for Au with an increment of 100K in each run of 35,000 integration step; but near the melting temperatures, the increment is reduced to 20K. The pressure applications are also implemented by repeating the thermal applications under pressure values of 0.5, 1.0, 1.5, 2.5, 5.0, 7.5, 10.0, 15.0 and 20.0 GPa. The simulation is

The temperature dependency of the elastic constants and the bulk moduli are calculated by

where <sup>1</sup> (x x x) ( )*<sup>t</sup> V* **σ b c, c a, a b h** and microscopic stress tensor, Π, is a dyadic given as

1 1

 

*N N <sup>N</sup> ij ii i i i i i j i ij*

. .

*<sup>N</sup>* ˆ*ij*

 

*j ij*

*F*

*r*

**<sup>Π</sup> vv rr** . (11)

*r*

**s F Δ** , (12)

**h**, *M* is an arbitrary constant

*ij ij ij r s Gs* . The classical equations of motion of the

**s F G Gs** (9)

**h ΠΙ σ** , (10)

MD cell axes, the metric tensor **G** is given by matrix product **h**<sup>t</sup>

between particles *i* and *j* is described by <sup>2</sup> *<sup>t</sup>*

system obtained from Eq. (1) become

follows;

arguments.

where *E*0 is a constant to be taken as an equilibrium cohesive energy of solid, *B*m is the bulk modulus, and is the atomic volume in equilibrium. It has been determined that the cohesive energy calculated from Eq. (1) with the parameter given in Table1 for Au is in good agreement with Rose energies in equilibrium.


Table 1. The experimental properties and potential parameters of Au. The experimental lattice parameters (*a*0) at room temperature are from ref. [27]. Bulk modulus (*B*m) and elastic constants (*Cij*) given at zero temperature are from [28], vacancy formation energy (*E*vf) is from ref. [29], melting temperature (*T*m), the coefficient of linear thermal expansion *α* are from [30], and specific heat *C*p is from [31].

Fig. 1. Rose and EAM energies versus lattice constant for Au.

#### **3. Molecular dynamics simulation**

The Lagrange function, written for an anisotropic box, i.e. MD cell, containing *N* particles by Parrinello and Rahman, is given by [33, 34]

$$L\_{PR} = \frac{1}{2} \sum\_{i=1}^{N} m\_i \left( \dot{\mathbf{s}}\_i^\dagger \mathbf{G} \dot{\mathbf{s}}\_i \right) - E\_{tot} + \frac{1}{2} M \text{Tr} \left( \dot{\mathbf{h}}^\dagger \dot{\mathbf{h}} \right) - P\_{\text{ext}} V \tag{8}$$

where *E*0 is a constant to be taken as an equilibrium cohesive energy of solid, *B*m is the bulk

cohesive energy calculated from Eq. (1) with the parameter given in Table1 for Au is in

Au 4.079 2.8842 3.81 0.93 180.32 201.63 169.67 45.44 1337 25.42

*D*<sup>1</sup> (eV)

Table 1. The experimental properties and potential parameters of Au. The experimental lattice parameters (*a*0) at room temperature are from ref. [27]. Bulk modulus (*B*m) and elastic constants (*Cij*) given at zero temperature are from [28], vacancy formation energy (*E*v<sup>f</sup>

from ref. [29], melting temperature (*T*m), the coefficient of linear thermal expansion *α* are

**Au**

3.95 4.00 4.05 4.10 4.15 4.20 4.25

 Rose EAM

a(A)

The Lagrange function, written for an anisotropic box, i.e. MD cell, containing *N* particles by

1 1 ( ) Tr( ) 2 2 *<sup>N</sup> t t*

**s Gs h h**

*L m E M PV*

ext

, (8)

*C*<sup>11</sup> (GPa)

is the atomic volume in equilibrium. It has been determined that the

*D*<sup>2</sup> (eV)

*C*<sup>12</sup> (Gpa)

*C*<sup>44</sup> (GPa)

*T*<sup>m</sup> (K)

*C*<sup>p</sup> (K/mol.K)

) is

modulus, and

*a0*  (Å)

*m n* 

*r0*  (Å)

good agreement with Rose energies in equilibrium.

*Ev* f (eV)

(Å1)

Au 7 0.5 4.3482 3.5361 0.0685 0.3097

*B*<sup>m</sup> (GPa)

*E*<sup>c</sup> (eV)

from [30], and specific heat *C*p is from [31].


**3. Molecular dynamics simulation** 

Parrinello and Rahman, is given by [33, 34]

Fig. 1. Rose and EAM energies versus lattice constant for Au.

1

*i*

*PR i i i tot*




E(eV)




where *mi* is mass of particle *i*, **s***i* is the scaled coordinate of atom *i* and is represented by a column vector whose elements are between zero and unity, **h**=(**a**, **b**, **c**); **a**, **b** and **c** vectors are MD cell axes, the metric tensor **G** is given by matrix product **h**<sup>t</sup> **h**, *M* is an arbitrary constant which represents mass of the computational box, *P*ext is external pressure applied on the cell, *V* is the volume of the MD cell and is obtained from det(**h**). Thus, square of distance between particles *i* and *j* is described by <sup>2</sup> *<sup>t</sup> ij ij ij r s Gs* . The classical equations of motion of the system obtained from Eq. (1) become

$$\ddot{\mathbf{s}}\_i = -\frac{1}{m\_i} \mathbf{F}\_i = \mathbf{G}^{-1} \dot{\mathbf{G}} \dot{\mathbf{s}}\_i \tag{9}$$

$$
\ddot{\mathbf{h}} = M^{-1} (\Pi - \mathbf{I} P\_{ext}) \mathbf{o} \tag{10}
$$

where <sup>1</sup> (x x x) ( )*<sup>t</sup> V* **σ b c, c a, a b h** and microscopic stress tensor, Π, is a dyadic given as follows;

$$\mathbf{III} = V^{-1} \left[ \sum\_{i=1}^{N} m\_i \mathbf{v}\_i, \mathbf{v}\_i - \sum\_{i=1}^{N} \sum\_{j>i}^{N} \frac{F\_{ij}}{r\_{ij}} \mathbf{r}\_i, \mathbf{r}\_i \right]. \tag{11}$$

Also the force on an atom *i* in the system is calculated from the following equation,

$$\mathbf{F}\_{i} = -\Delta\_{s} E\_{i} = -\sum\_{\substack{j=1 \\ j\neq i}}^{N} \left[ F\_{i}^{\prime} \rho\_{j}^{\prime} + F\_{j}^{\prime} \rho\_{i}^{\prime} + \phi^{\prime} i\_{ij} \right] \frac{\hat{\mathbf{s}}\_{ij}}{r\_{ij}} \,, \tag{12}$$

where the primes denote the first derivatives of the functions with respect to their arguments.

In all of the simulation studies, the equation of motion given in Eqs. (9) and (10) were numerically solved by using the velocity version of the Verlet algorithm [35]. The size of integration step was chosen to be 7.87x1015s for Au. Initial structures of the systems were constructed on a lattice with 1372 atoms and an FCC unit cell. It has been observed that, with these initial conditions, the systems were equilibrated in 5000 integration steps. Time averages of the thermodynamic properties of the system in each simulation run were determined by using 30,000 integration steps following the equilibration of the system. The structures of the system in solid phase were examined by using the radial distribution function. Melting temperatures were determined from the plots of the cohesive energy versus temperature. It is possible to classify our simulation runs in two groups as thermal and pressure applications. In the thermal applications, the temperature of the system under zero pressure is raised from 100K to 2400K for Au with an increment of 100K in each run of 35,000 integration step; but near the melting temperatures, the increment is reduced to 20K. The pressure applications are also implemented by repeating the thermal applications under pressure values of 0.5, 1.0, 1.5, 2.5, 5.0, 7.5, 10.0, 15.0 and 20.0 GPa. The simulation is restarted with different pressure in each run, to avoid algorithmic errors.

The temperature dependency of the elastic constants and the bulk moduli are calculated by following the procedure given by Karimi et al [14].

A Molecular Dynamics Study on Au 207

Considering the experimental data in Table 1, it can be seen that the specific heat is

Au P=0GPa Cp=28.2 J/m ol.K

Fig. 3. Variation of the enthalpy with temperature for Au.

5.6%, using the solid-liquid interface technique.

0 200 400 600 800 1000

**T(K)**

There are several methods for determining the melting temperature of a crystal. MD simulations are performed on system at various temperatures, and the cohesive energy is plotted as a function of temperature in one of these methods, as we did here. At the melting point, a discontinuity occurs in the cohesive energy. The other way of determining the melting temperature is to plot caloric curve which is the change of the total energy of crystal versus kinetic energy [36]. Indeed, the melting temperature of metal is obtained as the temperature at which the Gibbs free energy of the solid and liquid phases become equal. The entropy is required to compute the free energy, but it can not be directly calculated from MD simulations. For this reason, some other approaches are required [3]. Another way of determining the melting temperature is to simulate the solid-liquid interface [14]. In this way, the temperature for which the interface velocity goes to zero is determined as the melting temperature and it is reproduced more correctly than the way of caloric curve. Karimi et al [14] estimated the melting temperature for Ni as 1630±50K within an error of -

In the present work, the variations of cohesive energy with temperature for different pressures acted on the system are given in Fig. 4 for Au. We have computed the melting temperatures under zero pressure as 1100±20K for Au. When these values are compared

The radial distribution function (RDF) is used to investigate the structural properties of the solid and liquid phases. The plot of radial distribution functions acquired in solid and liquid phases for Au is given in Fig. 5. First peak location of radial distribution curves represents the distance of the nearest neighbor atoms, *r*0. The second peak location denotes the distances of next nearest neighbors, *a*0. These distances are found to be 2.907Å and 4.144Å, respectively for Au. By comparing with experimental data given in Table1, the calculated

with the experimental ones of 1337K given in Table 1, the error for Au becomes 21%.

calculated with an error of 9.8 % for Au.




**H(kJ/mol)**


For the calculation of glass formation and crystallization, firstly, we run 20 000 time steps to make the system into equilibrium state, then the liquid phase is cooled to 100K at the rate of 1.5833x1013 K/s and 1.5833x1012 K/s , respectively to examine the formation process of amorphization and crystallization.
