**1. Introduction**

24 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

Sato, M.; Yamataka, H.; Komeiji, Y.; Mochizuki, Y. & Nakano, T. (2010). Does Amination of

Sprik, M. & Ciccotti, G. (1998). Free energy from constrained molecular dynamics. *Journal of* 

Tamura, K.; Watanabe, T.; Ishimoto, T. & Nagashima, U. (2008). *Ab Initio* MO-MD

3765

2008), pp. 110–112, ISSN 1348-0634

Formaldehyde Proceed Through a Zwitterionic Intermediate in Water? Fragment Molecular Orbital Molecular Dynamics Simulations by Using Constraint Dynamics. *Chemistry-A European Journal*, Vol. 16, No. 22, (June 2010), pp. 6430-6433, ISSN 1521-

*Chemical Physics*, Vol. 109, No. 18, (November 1998), pp. 7737-7744, ISSN 1089-7690

Simulation Based on the Fragment MO Method: A Case of (-)-Epicatechin Gallate with STO-3G Basis Set. *Bulletin of Chemical Society of Japan,* Vol. 81, No. 1, (January

> Metal clusters and nanoparticles have gained attention in the recent years due to their application as catalysts, antimicrobials, pigments, micro circuits, drug delivery vectors, and many other uses. Many fascinating properties exhibited by nanomaterials are highly size and structure dependent. Therefore, understanding the formation of these nanoparticles is important in order to tailor their properties. The laboratory synthesis and characterization of such clusters and nanoparticles has provided insight into characteristics such as size and shape. However, monitoring the synthesis of such a cluster (or nanoparticle) on the atomic scale is difficult and to date no experimental technique is able to accomplish this. The use of computational methods has been employed to gain insight into the movement and interactions of atoms when a metal cluster or nanoparticle is formed. The most common computational approach has been to use molecular dynamics (MD) simulation which models the movement of atoms using a potential energy surface (PES) often referred to as a force field. The PES is used to describe the interaction of atoms and can be obtained from electronic strcuture calculations, from experimental measurements, or from the combining calculations and measurements.

> Molecular dynamics simulations have been used to study many phenomena associated with nanoparticles. Of particular interests are the geometric structure and energetics of nanoparticles of Au (Erkoc 2000; Shintani et al. 2004; Chui et al. 2007; Pu et al. 2010), Ag (El-Bayyari 1998; Monteil et al. 2010), Al (Yao et al. 2004), Fe (Boyukata et al. 2005), Pb (Hendy & Hall 2001), U (Erkoc et al. 1999) and of alloys such as NaMg (Dhavale et al. 1999), Pt-Ni/Co (Favry et al. 2011), Pt-Au (Mahboobi et al. 2009), Zn-Cd (Amirouche & Erkoc 2003), Cu-Ni/Pd (Kosilov et al. 2008), Co-Sb (Yang et al. 2011) as well as the behavior of nanoparticles during the melting or freezing process such as Au (Wang et al. 2005; Bas et al. 2006; Yildirim et al. 2007; Lin et al. 2010; Shibuta & Suzuki 2010), Na (Liu et al. 2009), Cu (Wang et al. 2003; Zhang et al. 2009), Al (Zhang et al. 2006), Fe (Ding et al. 2004; Shibuta & Suzuki 2008), Ni (Wen et al. 2004; Lyalin et al. 2009; Shibuta & Suzuki 2010), Pd (Miao et al. 2005), Sn (Chuang et al. 2004; Krishnamurty et al. 2006), Na-alloys (Aguado & Lopez 2005), Pt-alloys (Sankaranarayanan et al. 2005; Yang et al. 2008; Yang et al. 2009; Shi et al. 2011), Au-alloys (Yang et al. 2008; Yang et al. 2009; Gonzalez et al. 2011; Shi et al. 2011) and Ag-alloys (Kuntova et al. 2008; Kim et al. 2009). Molecular dynamics simulations have also been applied to study adsorption and desorption of nanoparticles on surfaces, such as Pd/MgO

(Long & Chen 2008) and Mn/Au (Mahboobi et al. 2010), nanoparticle aggregation such as Au (Lal et al. 2011), diffusion processes (Shimizu et al. 2001; Sawada et al. 2003; Yang et al. 2008; Alkis et al. 2009; Chen & Chang 2010), fragmentation of Au and Ag (Henriksson et al. 2005), thermal conductivity of Cu nanoparticles (Kang et al. 2011), and cluster (nanoparticle) formation of Au (Boyukata 2006; Cheng et al. 2009), Ag (Yukna & Wang 2007; Zeng et al. 2007; Hudson et al. 2010), Ir (Pawluk & Wang 2007), Co (Rives et al. 2008), and various alloys (Cheng et al. 2009; Chen & Chang 2010; Chen et al. 2010; Goniakowski & Mottet 2010; Carrillo & Dobrynin 2011).

Formation of metal clusters or nanoparticles can take place in all three phases: in liquid, gas, and on solid surfaces. Different formation mechanisms can be involved in the formation of transition metal nanoparticles. Of particular interest is coalescence, a process by which two droplets or particles collide to form a new daughter droplet or particle. Coalescence is important due to its role in nanoparticle formation and size control. Conventional MD simulations are used to describe coalescence of transition metal nanoparticles and provide information on the dynamics of nanoparticle formation, such as rate constant. However, the change of the electronic properties of the particles can only be probed by performing electronic structure calculations. Therefore, to have a complete picture of the formation of nanoparticles, the coupling of both MD and electronic structure calculations is important and forms the practice of our MD simulations. We denote it as the meta-molecular dynamics (meta-MD) method. In this chapter, we provide a description of the meta-MD method and its application in the study of Fe cluster formations. Before we present the meta-MD method and its application, we provide a general description of conventional MD simulations and the PES that is of ultimate importance in the accuracy of MD simulations.

### **2. Molecular Dynamics (MD) Simulations and Potential Energy Surfaces (PESs)**

In a conventional molecular dynamics simulation, if the motion of atoms in the system is governed by Newton's equations of motion, we numerically solve the position of atom *i* with a mass of *mi* in the Cartesian coordinates *xi*, *yi*, and *zi* by

$$\frac{\mathbf{dx}\_{i}}{\mathbf{dt}} = \frac{\mathbf{p}\_{\mathbf{x}\_{i}}}{\mathbf{m}\_{i}} \; \prime \tag{1a}$$

Advanced Molecular Dynamics Simulations on the Formation of Transition Metal Nanoparticles 27

dp V dt y

dp V dt z

The accuracy of the PES determines the accuracy of the outcome of MD simulations. There are many possible force fields (a.k.a. PESs) (Mazzone 2000; Hendy et al. 2003) but two used most often are the embedded atom method (Daw & Baskes 1984; Zhao et al. 2001; Dong et al. 2004; Lummen & Kraska 2004; Lummen & Kraska 2005; Lummen & Kraska 2005a, 2005b, 2005c; Rozas & Kraska 2007) and the Sutton-Chen potential (Kim et al. 2007; Pawluk &

i

i

i ij i <sup>1</sup> V [ v(i, j) c ] <sup>2</sup>

n

m

ij

ij a r

In the above equations of the Sutton-Chen potential, *rij* is the distance between atom *i* and

stick together. The attraction between atoms is normally described by a 1/r6 potential at long distances, due to van der Waals interaction, and described by an N-body form at short distances. Choosing the value of the parameter *m* to be 6 accomplishes these two things

<sup>n</sup> <sup>f</sup>

j

r

r  , is a cohesive term that describes the tendency for the atoms to

<sup>a</sup> v(i, j)

i j

atom *j*. The parameters; *a, n, m, ε,* and *c* depend on the element that is under study.

Define the lattice sum of a perfect face centered cubic (f.c.c.) crystal to be,

f n j

<sup>a</sup> <sup>S</sup>

The sum is taken over all separations *rj* from an arbitrary atom. *a*f is equal to the f.c.c. lattice

In equilibrium the total energy of the crystal does not change to first order when the lattice

, (2b)

. (2c)

, (3)

, (4)

. (5)

. (6)

yi

i z

Wang 2007; Yukna & Wang 2007; Hudson et al. 2010; Kayhani et al. 2010).

where *ν(i,j)* is an interaction between atoms *i* and *j* given by,

and *ρi* is the local electron density contribution of atom *i* given by,

In the Sutton-Chen PES, *V* is expressed as

The N-body term, i <sup>i</sup>

(Sutton & Chen 1990).

parameter which then defines the unit of length.

parameter is varied. This implies,

$$\frac{\mathbf{d}\mathbf{y}\_i}{\mathbf{d}t} = \frac{\mathbf{p}\_{\mathbf{y}\_i}}{\mathbf{m}\_i} \,\prime\tag{1b}$$

$$\frac{\mathbf{dz}\_i}{\mathbf{dt}} = \frac{\mathbf{p}\_{z\_i}}{\mathbf{m}\_i} \,. \tag{1c}$$

Here xi p , yi p , and zi p are the momentum of the atom *i* in the *x*, *y*, and *z* direction, respectively, and are solved by the gradient of the PES, denoted V:

$$\frac{\text{d}\mathbf{p}\_{\mathbf{x}\_{i}}}{\text{d}\mathbf{t}} = -\frac{\text{\textdegree\prime}}{\text{\textdegree\prime}}\text{\textdegree}\tag{2a}$$

$$\frac{\partial \mathbf{dp}\_{\mathbf{y}\_{i}}}{\partial \mathbf{dt}} = -\frac{\partial \mathbf{V}}{\partial \mathbf{y}\_{i}} \tag{2b}$$

$$\frac{\text{d}\mathbf{p}\_{\mathbf{z}\_{\parallel}}}{\text{d}\mathbf{t}} = -\frac{\text{\textdegree\textquotedblleft}\text{\textquotedblright}}{\text{\textquotedblleft}\mathbf{z}\_{\parallel}}.\tag{2c}$$

The accuracy of the PES determines the accuracy of the outcome of MD simulations. There are many possible force fields (a.k.a. PESs) (Mazzone 2000; Hendy et al. 2003) but two used most often are the embedded atom method (Daw & Baskes 1984; Zhao et al. 2001; Dong et al. 2004; Lummen & Kraska 2004; Lummen & Kraska 2005; Lummen & Kraska 2005a, 2005b, 2005c; Rozas & Kraska 2007) and the Sutton-Chen potential (Kim et al. 2007; Pawluk & Wang 2007; Yukna & Wang 2007; Hudson et al. 2010; Kayhani et al. 2010).

In the Sutton-Chen PES, *V* is expressed as

26 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

(Long & Chen 2008) and Mn/Au (Mahboobi et al. 2010), nanoparticle aggregation such as Au (Lal et al. 2011), diffusion processes (Shimizu et al. 2001; Sawada et al. 2003; Yang et al. 2008; Alkis et al. 2009; Chen & Chang 2010), fragmentation of Au and Ag (Henriksson et al. 2005), thermal conductivity of Cu nanoparticles (Kang et al. 2011), and cluster (nanoparticle) formation of Au (Boyukata 2006; Cheng et al. 2009), Ag (Yukna & Wang 2007; Zeng et al. 2007; Hudson et al. 2010), Ir (Pawluk & Wang 2007), Co (Rives et al. 2008), and various alloys (Cheng et al. 2009; Chen & Chang 2010; Chen et al. 2010; Goniakowski & Mottet 2010;

Formation of metal clusters or nanoparticles can take place in all three phases: in liquid, gas, and on solid surfaces. Different formation mechanisms can be involved in the formation of transition metal nanoparticles. Of particular interest is coalescence, a process by which two droplets or particles collide to form a new daughter droplet or particle. Coalescence is important due to its role in nanoparticle formation and size control. Conventional MD simulations are used to describe coalescence of transition metal nanoparticles and provide information on the dynamics of nanoparticle formation, such as rate constant. However, the change of the electronic properties of the particles can only be probed by performing electronic structure calculations. Therefore, to have a complete picture of the formation of nanoparticles, the coupling of both MD and electronic structure calculations is important and forms the practice of our MD simulations. We denote it as the meta-molecular dynamics (meta-MD) method. In this chapter, we provide a description of the meta-MD method and its application in the study of Fe cluster formations. Before we present the meta-MD method and its application, we provide a general description of conventional MD simulations and

the PES that is of ultimate importance in the accuracy of MD simulations.

with a mass of *mi* in the Cartesian coordinates *xi*, *yi*, and *zi* by

respectively, and are solved by the gradient of the PES, denoted V:

**2. Molecular Dynamics (MD) Simulations and Potential Energy Surfaces** 

In a conventional molecular dynamics simulation, if the motion of atoms in the system is governed by Newton's equations of motion, we numerically solve the position of atom *i*

> <sup>i</sup> i x i

> i yi i

<sup>i</sup> i z i

Here xi p , yi p , and zi p are the momentum of the atom *i* in the *x*, *y*, and *z* direction,

dp V dt x

i

i x

dt m , (1a)

dt m , (1b)

dt m . (1c)

, (2a)

dx p

dy p

dz p

Carrillo & Dobrynin 2011).

**(PESs)** 

$$\mathbf{V} = \varepsilon [\frac{1}{2} \sum\_{\mathbf{i}\uparrow} \mathbf{v}(\mathbf{i}, \mathbf{j}) - \mathbf{c} \sum\_{\mathbf{i}} \sqrt{\rho\_{\mathbf{i}}}],\tag{3}$$

where *ν(i,j)* is an interaction between atoms *i* and *j* given by,

$$\mathbf{v}(\mathbf{i}, \mathbf{j}) = \left(\frac{\mathbf{a}}{\mathbf{r}\_{\parallel}}\right)^{n},\tag{4}$$

and *ρi* is the local electron density contribution of atom *i* given by,

$$
\rho\_{\parallel} = \sum\_{\parallel} \left( \frac{\mathbf{a}}{\mathbf{r\_{\parallel}}} \right)^{\mathbf{m}}.\tag{5}
$$

In the above equations of the Sutton-Chen potential, *rij* is the distance between atom *i* and atom *j*. The parameters; *a, n, m, ε,* and *c* depend on the element that is under study.

The N-body term, i <sup>i</sup> , is a cohesive term that describes the tendency for the atoms to stick together. The attraction between atoms is normally described by a 1/r6 potential at long distances, due to van der Waals interaction, and described by an N-body form at short distances. Choosing the value of the parameter *m* to be 6 accomplishes these two things (Sutton & Chen 1990).

Define the lattice sum of a perfect face centered cubic (f.c.c.) crystal to be,

$$\mathbf{S}\_{n}^{\ell} = \sum\_{\parallel} \left( \frac{\mathbf{a}^{\ell}}{\mathbf{r}\_{\parallel}} \right)^{n}. \tag{6}$$

The sum is taken over all separations *rj* from an arbitrary atom. *a*f is equal to the f.c.c. lattice parameter which then defines the unit of length.

In equilibrium the total energy of the crystal does not change to first order when the lattice parameter is varied. This implies,

$$\mathbf{c} = \frac{\mathbf{n} \mathbf{s}\_n^\ell}{\mathbf{m} \sqrt{\mathbf{s}\_m^\ell}}.\tag{7}$$

Advanced Molecular Dynamics Simulations on the Formation of Transition Metal Nanoparticles 29

Once these theories are chosen, a coupling method has to be employed so that the two types of calculations can be integrated. Appropriate techniques need to be developed in order to integrate the electronic structure calculations seamlessly to the MD simulations. There are several ways to couple MD and DFT calculations. The most straightforward way would be to perform MD simulations first and save the structural information, i.e. the Cartesian coordinates of each atom at each time step. The DFT calculations can be performed using these data. We note that the time step in a typical MD simulation is in the range of 0.01-1 fs, and the simulation can run for ~10ps or longer times. Therefore, the concern with this strategy is that far too many data have to be saved. Additionally, too much computational time is required. An alternative strategy is to carry out MD and DFT calculations simultaneously. One of the advantages of this strategy is to be able to use the electronic wavefunctions generated from the past time step as the initial electronic wavefunction in the subsequent DFT calculation. We are exploring other possibilities to save computer time. Further, the time between two DFT calculations will be set to a longer interval than a simple MD time step. A DFT calculation will be performed between the two DFT calculations only when significant changes have taken place, which will be monitored in the MD simulations. For a particular system of interest, one will need to test extensively which length of time interval will be appropriate to perform DFT calculations in order to find an optimal choice and a common ground in terms of effectiveness

and accuracy. Development of the pragmatic coupling methods is in progress.

our studies of Ag cluster formation previously (Hudson et al. 2010).

McCoy 2003).

In this work, we studied the coalescence of an Fe dimer with an Fe atom and of two 15-atom Fe clusters. The numerical aspects of the MD simulations are similar to the previous MD simulations (Billing & Wang 1992a, 1992b; Ge et al. 1992; Wang & Billing 1992; Wang & Billing 1993; Wang et al. 1994a, 1994b; Wang & Clary 1996; Clary & Wang 1997; Wang & Billing 1997; Wang et al. 1999; Wang et al. 2000; McCoy et al. 2001; Wang &

The simulations for the systems were performed over a 10 ps time period with a time step of 0.01 fs. During the formation of nanoparticles, a small amount of energy was extracted at every single time step from any atom that had a kinetic energy greater than a defined minimum energy. In this work two simulations were run with a minimum of 298 K and 693 K. The subtraction of energy was done to mimic a cooling rate of 1.5625x1011 K/s (1.3464x10-8 eV/fs) and 1.5625x1013 K/s (1.3454x10-8 eV/fs). These cooling rates were used in

DFT calculations were performed in a similar fashion as our previous studies of transition metal clusters (Wang & Ge 2002; Cao et al. 2003; Zhang et al. 2003; Xiao & Wang 2004a, 2004b; Zhang et al. 2004a, 2004b, 2004c) Specifically, spin polarized DFT calculations were carried out via the Vienna Ab-initio Simulation Package (VASP).(Kresse & Hafner 1993; Kresse & Furthmuller 1996a, 1996b) The electron-ion interactions were described by the Projector Augmented Waves method.(Kresse & Joubert 1999) The exchange and correlation energies were calculated using the Perdew-Burke-Ernzerhof (PBE) functional.(Perdew et al. 1996) A plane wave basis set was used with a cutoff energy of 300 eV, which was shown to be sufficient from the convergence test. One k point, the point, was used. In order to eliminate interactions between two neighboring images, we set the nearest distance between images no less than 1.0 nm. The simulation techniques used here are very similar to those in our previous study of Pt clusters.(Xiao & Wang 2004b) Single point calculations were performed based on the structural data from the MD simulations. The binding energy, the

The cohesive energy per atom is given by,

$$\mathbf{E}\_c^\ddagger = \frac{\varepsilon \mathbf{s}\_n^\ddagger}{2\mathbf{m}} (2\mathbf{n} - \mathbf{m}) \tag{8}$$

Finally the bulk modulus, *Bf* , is given by,

$$\mathbf{B}^{\text{f}} = \frac{(\mathbf{2}\mathbf{n} - \mathbf{m})\mathbf{n}\varepsilon\mathbf{s}\_{\mathbf{n}}^{\text{f}}}{36\angle\Omega^{\text{f}}}\_{\text{f}} \tag{9}$$

where *Ω<sup>f</sup> =(af )*3/4 which is the atomic volume. Using the above equations a relation between the cohesive energy, *E*, and the bulk modulus, *B* is given by,

$$\frac{\mathbf{\color{red}{ $\mathcal{Q}$ }}^{\ell}\mathbf{\color{red}{ $\mathcal{B}$ }}^{\ell}}{\mathbf{\color{red}{ $\mathcal{E}\_{\text{c}}$ }}^{\ell}} = \frac{\mathbf{\color{red}{\mathbf{mm}}}}{\mathbf{18}} \,. \tag{10}$$

Using experimental measurements of the cohesive energy, *E*, the bulk modulus, *B*, and the chosen value of *m*=6, an integer value for *n* was to give the value closest in agreement with eq. (10). From the values of *m* and *n* eq. (9) can be used to obtain a value of *ε* and eq. (8) can be used to obtain the value for *c*. (Sutton & Chen 1990) The parameters *ε* and *a*  are defined as units of energy and length, respectively. Thus the values of *ε* and *a* were chosen, for the different metals, to coincide with results obtained from fitting the PES with experimental or computational measurements. The parameters used in the current MD simulations of Fe cluster formations were obtained from Sutton and Chen (Sutton & Chen 1990).

#### **3. Advanced MD simulations: Meta-MD simulations**

In the meta-MD simulations, we couple the conventional MD simulation with the electronic structure calculation to study the formation of transition metal nanoparticles. Such a coupling allows us to record the electronic change of the system during the formation process in addition to the conventional properties in a MD simulation. Furthermore, we will also be able to monitor the accuracy of the PES as well as determine whether the MD simulation on a single PES is valid.

The three ingredients in a meta-MD simulation are electronic structure theory, molecular dynamics theory, and coupling method. In principle, any electronic structure theory can be chosen. Depending on the system of interest, our choice of a particular electronic structure theory is determined by the cost effectiveness and the accuracy of electronic structure calculations. For transition metal systems, the most practical choice of method is density functional theory, where a variety of functionals may be used. The molecular dynamics theory can be quantum scattering, pure classical, mixed quantum-classical, or semi-classical treatment, which also depends on the characteristics of the system to be described. For instance, our current system involves only heavy atoms, we therefore choose classical MD simulations.

ns

m s

f f n c

<sup>s</sup> E (2n m) 2m 

f n f

(2n m)n s <sup>B</sup> 36

f f f c

B nm E 18

Using experimental measurements of the cohesive energy, *E*, the bulk modulus, *B*, and the chosen value of *m*=6, an integer value for *n* was to give the value closest in agreement with eq. (10). From the values of *m* and *n* eq. (9) can be used to obtain a value of *ε* and eq. (8) can be used to obtain the value for *c*. (Sutton & Chen 1990) The parameters *ε* and *a*  are defined as units of energy and length, respectively. Thus the values of *ε* and *a* were chosen, for the different metals, to coincide with results obtained from fitting the PES with experimental or computational measurements. The parameters used in the current MD simulations of Fe cluster formations were obtained from Sutton and Chen (Sutton &

In the meta-MD simulations, we couple the conventional MD simulation with the electronic structure calculation to study the formation of transition metal nanoparticles. Such a coupling allows us to record the electronic change of the system during the formation process in addition to the conventional properties in a MD simulation. Furthermore, we will also be able to monitor the accuracy of the PES as well as determine whether the MD

The three ingredients in a meta-MD simulation are electronic structure theory, molecular dynamics theory, and coupling method. In principle, any electronic structure theory can be chosen. Depending on the system of interest, our choice of a particular electronic structure theory is determined by the cost effectiveness and the accuracy of electronic structure calculations. For transition metal systems, the most practical choice of method is density functional theory, where a variety of functionals may be used. The molecular dynamics theory can be quantum scattering, pure classical, mixed quantum-classical, or semi-classical treatment, which also depends on the characteristics of the system to be described. For instance, our current system involves only heavy atoms, we therefore choose classical MD

f

*)*3/4 which is the atomic volume. Using the above equations a relation between

c

, is given by,

the cohesive energy, *E*, and the bulk modulus, *B* is given by,

**3. Advanced MD simulations: Meta-MD simulations** 

The cohesive energy per atom is given by,

Finally the bulk modulus, *Bf*

*=(af*

where *Ω<sup>f</sup>*

Chen 1990).

simulations.

simulation on a single PES is valid.

f n f m

. (7)

. (10)

. (8)

, (9)

Once these theories are chosen, a coupling method has to be employed so that the two types of calculations can be integrated. Appropriate techniques need to be developed in order to integrate the electronic structure calculations seamlessly to the MD simulations. There are several ways to couple MD and DFT calculations. The most straightforward way would be to perform MD simulations first and save the structural information, i.e. the Cartesian coordinates of each atom at each time step. The DFT calculations can be performed using these data. We note that the time step in a typical MD simulation is in the range of 0.01-1 fs, and the simulation can run for ~10ps or longer times. Therefore, the concern with this strategy is that far too many data have to be saved. Additionally, too much computational time is required. An alternative strategy is to carry out MD and DFT calculations simultaneously. One of the advantages of this strategy is to be able to use the electronic wavefunctions generated from the past time step as the initial electronic wavefunction in the subsequent DFT calculation. We are exploring other possibilities to save computer time. Further, the time between two DFT calculations will be set to a longer interval than a simple MD time step. A DFT calculation will be performed between the two DFT calculations only when significant changes have taken place, which will be monitored in the MD simulations. For a particular system of interest, one will need to test extensively which length of time interval will be appropriate to perform DFT calculations in order to find an optimal choice and a common ground in terms of effectiveness and accuracy. Development of the pragmatic coupling methods is in progress.

In this work, we studied the coalescence of an Fe dimer with an Fe atom and of two 15-atom Fe clusters. The numerical aspects of the MD simulations are similar to the previous MD simulations (Billing & Wang 1992a, 1992b; Ge et al. 1992; Wang & Billing 1992; Wang & Billing 1993; Wang et al. 1994a, 1994b; Wang & Clary 1996; Clary & Wang 1997; Wang & Billing 1997; Wang et al. 1999; Wang et al. 2000; McCoy et al. 2001; Wang & McCoy 2003).

The simulations for the systems were performed over a 10 ps time period with a time step of 0.01 fs. During the formation of nanoparticles, a small amount of energy was extracted at every single time step from any atom that had a kinetic energy greater than a defined minimum energy. In this work two simulations were run with a minimum of 298 K and 693 K. The subtraction of energy was done to mimic a cooling rate of 1.5625x1011 K/s (1.3464x10-8 eV/fs) and 1.5625x1013 K/s (1.3454x10-8 eV/fs). These cooling rates were used in our studies of Ag cluster formation previously (Hudson et al. 2010).

DFT calculations were performed in a similar fashion as our previous studies of transition metal clusters (Wang & Ge 2002; Cao et al. 2003; Zhang et al. 2003; Xiao & Wang 2004a, 2004b; Zhang et al. 2004a, 2004b, 2004c) Specifically, spin polarized DFT calculations were carried out via the Vienna Ab-initio Simulation Package (VASP).(Kresse & Hafner 1993; Kresse & Furthmuller 1996a, 1996b) The electron-ion interactions were described by the Projector Augmented Waves method.(Kresse & Joubert 1999) The exchange and correlation energies were calculated using the Perdew-Burke-Ernzerhof (PBE) functional.(Perdew et al. 1996) A plane wave basis set was used with a cutoff energy of 300 eV, which was shown to be sufficient from the convergence test. One k point, the point, was used. In order to eliminate interactions between two neighboring images, we set the nearest distance between images no less than 1.0 nm. The simulation techniques used here are very similar to those in our previous study of Pt clusters.(Xiao & Wang 2004b) Single point calculations were performed based on the structural data from the MD simulations. The binding energy, the

Advanced Molecular Dynamics Simulations on the Formation of Transition Metal Nanoparticles 31

trigonal structure. The slow cooling rate also shows a wild oscillation of the bond distances between atoms. This result is as expected due to the slower removal of energy from the system. It can also be noticed that the higher minimum energy resulted in a more violent oscillation of bond lengths after cluster formation (Fig. 2 vs. Fig. 1 and Fig. 4 vs. Fig. 3). This is because the greater amount of heat available is expressed as a vibration in the formed cluster. The greater initial energy of the system (0.1 eV vs. 0.5 eV) has an effect only when a

Fig. 2. The interatomic distances of a three-atom coalescence with an initial energy of 0.1 eV at a minimum temperature of 673 K and a cooling rate of 1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right) of atoms 1 & 2 (red/short dash), 2 & 3 (green/solid), and 1 & 3 (blue/long

Fig. 3. The interatomic distances of a three-atom coalescence with an initial energy of 0.5 eV at a minimum temperature of 298 K and a cooling rate of 1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right) of atoms 1 & 2 (red/short dash), 2 & 3 (green/solid), and 1 & 3 (blue/long

slow cooling rate is employed.

dash).

dash).

energy gap between the highest occupied molecular orbital and the lowest unoccupied molecular orbital (HOMO-LUMO), and the magnetic moment of the system were obtained and discussed.

#### **4. Formation of iron clusters**

In this work, we performed MD simulations to study the formation of iron trimers and meta-MD simulations to that of 30-atom iron clusters. The results of these simulations are presented below starting from the formation of iron trimers.

#### **4.1 Iron trimer formation**

The MD results for the Fe atom collinearly colliding with an Fe dimer are summarized in Figs. 1-8.

We first examine the effects of initial kinetic energy, cooling rate, and temperature on the structure of Fe trimers. Figure 1 shows the time evolution of interatomic distances between Fe pairs. These plots show the slower cooling rate producing a faster collision and a longer lived linear trimer than the faster cooling rate, though the final products in both cooling rates are triangular trimers, which are demonstrated by the overlap of all three curves at longer times. When the minimum temperature increases from 298 K to 673 K, the linear trimers at both cooling rates last longer, as clearly shown in Fig. 2. In fact, the linear trimer exists at the end of simulation time when the cooling rate is 1.5625x1011 K/s. When the initial kinetic energy increases from 0.1 eV shown in Figs. 1 and 2 to 0.5 eV shown in Figs. 3 and 4, the time evolution of interatomic distances has trends similar to the lower kinetic energy cases.

Fig. 1. The interatomic distances of a three-atom coalescence with an initial energy of 0.1 eV at a minimum temperature of 298 K and a cooling rate of 1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right) of atoms 1 & 2 (red/short dash), 2 & 3 (green/solid), and 1 & 3 (blue/long dash).

Among the eight collisions depicted in Figs. 1-4, the favored product was the trigonal trimer, which accounts for 6 of them. The other two were the linear trimer at the end of 10 ps. A high minimum energy and a slow cooling rate was the condition that gave the linear trimer regardless of the initial energy given to the system. The slow cooling rate resulted in a longer duration of the linear trimer configuration before the structure converted to the

energy gap between the highest occupied molecular orbital and the lowest unoccupied molecular orbital (HOMO-LUMO), and the magnetic moment of the system were obtained

In this work, we performed MD simulations to study the formation of iron trimers and meta-MD simulations to that of 30-atom iron clusters. The results of these simulations are

The MD results for the Fe atom collinearly colliding with an Fe dimer are summarized in

We first examine the effects of initial kinetic energy, cooling rate, and temperature on the structure of Fe trimers. Figure 1 shows the time evolution of interatomic distances between Fe pairs. These plots show the slower cooling rate producing a faster collision and a longer lived linear trimer than the faster cooling rate, though the final products in both cooling rates are triangular trimers, which are demonstrated by the overlap of all three curves at longer times. When the minimum temperature increases from 298 K to 673 K, the linear trimers at both cooling rates last longer, as clearly shown in Fig. 2. In fact, the linear trimer exists at the end of simulation time when the cooling rate is 1.5625x1011 K/s. When the initial kinetic energy increases from 0.1 eV shown in Figs. 1 and 2 to 0.5 eV shown in Figs. 3 and 4, the time evolution of interatomic distances has trends similar to

Fig. 1. The interatomic distances of a three-atom coalescence with an initial energy of 0.1 eV at a minimum temperature of 298 K and a cooling rate of 1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right) of atoms 1 & 2 (red/short dash), 2 & 3 (green/solid), and 1 & 3 (blue/long dash).

Among the eight collisions depicted in Figs. 1-4, the favored product was the trigonal trimer, which accounts for 6 of them. The other two were the linear trimer at the end of 10 ps. A high minimum energy and a slow cooling rate was the condition that gave the linear trimer regardless of the initial energy given to the system. The slow cooling rate resulted in a longer duration of the linear trimer configuration before the structure converted to the

and discussed.

Figs. 1-8.

**4. Formation of iron clusters** 

**4.1 Iron trimer formation** 

the lower kinetic energy cases.

presented below starting from the formation of iron trimers.

trigonal structure. The slow cooling rate also shows a wild oscillation of the bond distances between atoms. This result is as expected due to the slower removal of energy from the system. It can also be noticed that the higher minimum energy resulted in a more violent oscillation of bond lengths after cluster formation (Fig. 2 vs. Fig. 1 and Fig. 4 vs. Fig. 3). This is because the greater amount of heat available is expressed as a vibration in the formed cluster. The greater initial energy of the system (0.1 eV vs. 0.5 eV) has an effect only when a slow cooling rate is employed.

Fig. 2. The interatomic distances of a three-atom coalescence with an initial energy of 0.1 eV at a minimum temperature of 673 K and a cooling rate of 1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right) of atoms 1 & 2 (red/short dash), 2 & 3 (green/solid), and 1 & 3 (blue/long dash).

Fig. 3. The interatomic distances of a three-atom coalescence with an initial energy of 0.5 eV at a minimum temperature of 298 K and a cooling rate of 1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right) of atoms 1 & 2 (red/short dash), 2 & 3 (green/solid), and 1 & 3 (blue/long dash).

Advanced Molecular Dynamics Simulations on the Formation of Transition Metal Nanoparticles 33

energy plots, the potential energy plots show that the slower cooling rate produces a more oscillatory potential while the faster cooling rate produces a smooth curve with a sudden

When the initial kinetic energy increases to 0.5 eV, the energy distributions at different cooling rate are similar to the case of 0.1 eV. There is a resemblance of the kinetic and potential energy plots between the two minimum temperatures, namely Fig. 7 *vs* Fig. 8. The difference between the temperature lies at the oscillatory part of the potential energy curves. In the case of the 298 K chamber, the oscillations in the potential energy curve occur at the

Fig. 6. The kinetic (red/dash) and potential (green/solid) energy of a three-atom coalescence with an initial energy of 0.1 eV at a minimum temperature of 673 K and a cooling rate of

Fig. 7. The kinetic (red/dash) and potential (green/solid) energy of a three-atom coalescence with an initial energy of 0.5 eV at a minimum temperature of 298 K and a cooling rate of

potential drop.

very last of the MD simulations.

1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right).

1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right).

Fig. 4. The interatomic distances of a three-atom coalescence with an initial energy of 0.5 eV at a minimum temperature of 673 K and a cooling rate of 1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right) of atoms 1 & 2 (red/short dash), 2 & 3 (green/solid), and 1 & 3 (blue/long dash).

We now discuss the energetic aspects of the MD results. Figure 5 shows the changes of the kinetic and potential energy over time at the minimum temperature of 298 K and initial kinetic energy of 0.1 eV. The kinetic energy oscillates during the collision with the slower cooling rate (red curve of the left figure) and is essentially featureless in the case with the faster cooling rate (red curve of right figure). This is due to the slower cooling rate not being able to dissipate the kinetic energy released by the rapid decrease in potential energy.

Fig. 5. The kinetic (red/dash) and potential (green/solid) energy of a three-atom coalescence with an initial energy of 0.1 eV at a minimum temperature of 298 K and a cooling rate of 1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right).

When the minimum energy, or similarly the reaction chamber temperature, was 673 K, similar pictures of the coalescences were obtained. The kinetic energy oscillates in a regular pattern (left of Fig. 6) for the slower cooling rate and has no feature for the faster cooling rate (right of Fig. 6). Again, this is due to the slower cooling rate not being able to remove all of the kinetic energy gained due to a rapid release of potential energy. Similar to the kinetic

Fig. 4. The interatomic distances of a three-atom coalescence with an initial energy of 0.5 eV at a minimum temperature of 673 K and a cooling rate of 1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right) of atoms 1 & 2 (red/short dash), 2 & 3 (green/solid), and 1 & 3 (blue/long

We now discuss the energetic aspects of the MD results. Figure 5 shows the changes of the kinetic and potential energy over time at the minimum temperature of 298 K and initial kinetic energy of 0.1 eV. The kinetic energy oscillates during the collision with the slower cooling rate (red curve of the left figure) and is essentially featureless in the case with the faster cooling rate (red curve of right figure). This is due to the slower cooling rate not being able to dissipate the kinetic energy released by the rapid decrease in potential energy.

Fig. 5. The kinetic (red/dash) and potential (green/solid) energy of a three-atom coalescence with an initial energy of 0.1 eV at a minimum temperature of 298 K and a cooling rate of

When the minimum energy, or similarly the reaction chamber temperature, was 673 K, similar pictures of the coalescences were obtained. The kinetic energy oscillates in a regular pattern (left of Fig. 6) for the slower cooling rate and has no feature for the faster cooling rate (right of Fig. 6). Again, this is due to the slower cooling rate not being able to remove all of the kinetic energy gained due to a rapid release of potential energy. Similar to the kinetic

1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right).

dash).

energy plots, the potential energy plots show that the slower cooling rate produces a more oscillatory potential while the faster cooling rate produces a smooth curve with a sudden potential drop.

When the initial kinetic energy increases to 0.5 eV, the energy distributions at different cooling rate are similar to the case of 0.1 eV. There is a resemblance of the kinetic and potential energy plots between the two minimum temperatures, namely Fig. 7 *vs* Fig. 8. The difference between the temperature lies at the oscillatory part of the potential energy curves. In the case of the 298 K chamber, the oscillations in the potential energy curve occur at the very last of the MD simulations.

Fig. 6. The kinetic (red/dash) and potential (green/solid) energy of a three-atom coalescence with an initial energy of 0.1 eV at a minimum temperature of 673 K and a cooling rate of 1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right).

Fig. 7. The kinetic (red/dash) and potential (green/solid) energy of a three-atom coalescence with an initial energy of 0.5 eV at a minimum temperature of 298 K and a cooling rate of 1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right).

Advanced Molecular Dynamics Simulations on the Formation of Transition Metal Nanoparticles 35

Fig. 9. Energy plots of the coalescence between two 15-atom clusters with an initial energy of

The potential energy plot (Fig. 9 right) shows that the slow cooling rate has a more erratic variation in potential, yet three out of the four simulation conditions gave similar potentials. The one simulation to give a higher potential than the others had a fast cooling rate and low minimum temperature. This is likely due to the cluster being stuck in a higher potential configuration and not being able to overcome an energy barrier to reach a lower energy

Fig. 10. 30-atom clusters formed under 0.5 eV initial energy with 298 K minimum

**(a) (b)** 

temperature and a cooling rate of 1.5625 x 1011 K/s (left, a) and 1.5625 x 1013 K/s (right, b).

0.5 eV and a minimum temperature and cooling rate of 298 K and 1.5625 x 1011 K/s (red/short dash), 298 K and 1.5625 x 1013 K/s (green/solid), 673 K and 1.5625 x 1011 K/s

(blue/long dash), and 673 K and 1.5625 x 1013 (yellow/dotted).

configuration.

Fig. 8. The kinetic (red/dash) and potential (green/solid) energy of a three-atom coalescence with an initial energy of 0.5 eV at a minimum temperature of 673 K and a cooling rate of 1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right).

Figures 5-8 show that the faster cooling rate generates a smoother kinetic and potential energy plot. The minimum energy has little noticeable effect on the kinetic and potential energy plots except that the oscillations of either are more drastic with the greater minimum energy. The greater initial energy causes a spike in the kinetic energy at the beginning of the simulation which occurs later in the simulation when less initial energy is given to the system. Figures 7 and 8 both show a 'flare up' that is noticeably separate from the initial spike of kinetic energy when a slow cooling rate is used. The simulations that resulted in trigonal trimers gave an ending potential energy of around -3.5 eV while the simulations that resulted in the linear trimer (Fig. 6 left and Fig. 8 left) gave an ending potential of around -3.0 eV. This is due to the lower potential when each atom interacts with the other two atoms rather than two of the three atoms only interacting with one other atom.

#### **4.2 Formation of 30-atom iron clusters**

Four MD simulations were carried out with different minimum temperatures and cooling rates in order to investigate how these factors affect the formation process and the structure of the products.

The energy profile evolutions and the final structures of the 30-atom Fe clusters are given in Figs. 9-11.

Figure 9 shows the kinetic energy plot (left) and the potential energy plot (right) of the coalescence of two 15- atom clusters. The slow cooling rate gave a kinetic energy spike after the initial reaction (Fig. 9, left blue and red). The higher minimum temperature results in a higher kinetic energy at the end of the simulation (Fig. 9, left blue and yellow). Figure 9 also shows that the faster cooling rate negates the other parameters. The faster cooling rate plots (yellow and green) have similar kinetic energies while the slower cooling rate plots (red and blue) are very different and the difference is determined by the minimum temperature.

Fig. 8. The kinetic (red/dash) and potential (green/solid) energy of a three-atom coalescence with an initial energy of 0.5 eV at a minimum temperature of 673 K and a cooling rate of

Figures 5-8 show that the faster cooling rate generates a smoother kinetic and potential energy plot. The minimum energy has little noticeable effect on the kinetic and potential energy plots except that the oscillations of either are more drastic with the greater minimum energy. The greater initial energy causes a spike in the kinetic energy at the beginning of the simulation which occurs later in the simulation when less initial energy is given to the system. Figures 7 and 8 both show a 'flare up' that is noticeably separate from the initial spike of kinetic energy when a slow cooling rate is used. The simulations that resulted in trigonal trimers gave an ending potential energy of around -3.5 eV while the simulations that resulted in the linear trimer (Fig. 6 left and Fig. 8 left) gave an ending potential of around -3.0 eV. This is due to the lower potential when each atom interacts with the other

two atoms rather than two of the three atoms only interacting with one other atom.

Four MD simulations were carried out with different minimum temperatures and cooling rates in order to investigate how these factors affect the formation process and the structure

The energy profile evolutions and the final structures of the 30-atom Fe clusters are given in

Figure 9 shows the kinetic energy plot (left) and the potential energy plot (right) of the coalescence of two 15- atom clusters. The slow cooling rate gave a kinetic energy spike after the initial reaction (Fig. 9, left blue and red). The higher minimum temperature results in a higher kinetic energy at the end of the simulation (Fig. 9, left blue and yellow). Figure 9 also shows that the faster cooling rate negates the other parameters. The faster cooling rate plots (yellow and green) have similar kinetic energies while the slower cooling rate plots (red and blue) are very different and the difference is determined by the

1.5625 x 1011 K/s (left) and 1.5625 x 1013 K/s (right).

**4.2 Formation of 30-atom iron clusters** 

of the products.

minimum temperature.

Figs. 9-11.

Fig. 9. Energy plots of the coalescence between two 15-atom clusters with an initial energy of 0.5 eV and a minimum temperature and cooling rate of 298 K and 1.5625 x 1011 K/s (red/short dash), 298 K and 1.5625 x 1013 K/s (green/solid), 673 K and 1.5625 x 1011 K/s (blue/long dash), and 673 K and 1.5625 x 1013 (yellow/dotted).

The potential energy plot (Fig. 9 right) shows that the slow cooling rate has a more erratic variation in potential, yet three out of the four simulation conditions gave similar potentials. The one simulation to give a higher potential than the others had a fast cooling rate and low minimum temperature. This is likely due to the cluster being stuck in a higher potential configuration and not being able to overcome an energy barrier to reach a lower energy configuration.

Fig. 10. 30-atom clusters formed under 0.5 eV initial energy with 298 K minimum temperature and a cooling rate of 1.5625 x 1011 K/s (left, a) and 1.5625 x 1013 K/s (right, b).

Advanced Molecular Dynamics Simulations on the Formation of Transition Metal Nanoparticles 37

Meta-Molecular Dynamics (meta-MD) simulation was developed and described for studying the formation of transition metal nanoparticles. The meta-MD simulation integrates single point electronic structure calculations into the conventional molecular dynamics simulations so that instant changes of the intrinsic electromagnetic properties of the system can be monitored and obtained during the formation of nanoparticles. The results of Fe cluster formation obtained from the meta-MD simulations were presented and discussed. Additionally, the effect of cooling rates was also presented and discussed. Furthermore, using the spin-polarized DFT calculations in meta-MD simulations can also provide indications whether the electronically adiabatic treatment in the MD simulations is

sufficient by monitoring the electronic state changes during the dynamic processes.

The meta-MD technique developed here should also be a good tool in studying heterogeneous catalysis by providing guidance in the design of catalysts. For instance, the detailed picture of local charge distribution may provide insight into the active site and requirement for the catalytic activity. This information will be potentially useful in the preparation of catalysts. The meta-MD simulations described here can also be employed for studying other processes where the changes of the intrinsic electromagnetic properties of the local entities of the system, i.e. subsystems, are important in order to obtain a complete

This work is supported by the National Science Foundation (Grant CBET-0709113) and in

Aguado, A. & Lopez, J. M. (2005). Molecular dynamics simulations of the meltinglike

Alkis, S., Krause, J. L., Fry, J. N., Cheng, H. P. (2009). Dynamics of Ag clusters on complex

Amirouche, L. & Erkoc, S. (2003). Structural features and energetics of Znn-mCdm (n=7,8)

Bas, B. S. d., Ford, M. J., Cortie, M. B. (2006). Melting in small gold clusters: a density functional molecular dynamics study. *J. Phys.: Condens. Matter* 18: 55-74. Billing, G. D. & Wang, L. (1992a). Semiclassical calculations of transport coefficients and rotational relaxation of nitrogen at high temperatures. *J. Phys. Chem.* 96: 2572257-5. Billing, G. D. & Wang, L. (1992b). The use of stratified important sampling for calculating

Boyukata, M. (2006). Molecular-dynamics study of possible packing sequence of medium

Boyukata, M., Borges, E., Braga, J. P., Belchior, J. C. (2005). Size evolution of structures and

Lennard-Jones type potential. *J. Alloys Compounds* 403: 349-356.

energetics of iron clusters (Fe*n*, *n*< 36): Molecular dynamics studies using a

microclusters and Zn50, Cd50, and Zn25Cd25 nanoparticles: Molecular-dynamics

transition in Li13Na42 and Na13Cs42 clusters. *Phys. Rev. B* 71: 075415.

surfaces: Molecular dynamics simulations. *Phys. Rev. B* 79: 121402.

**5. Conclusion** 

picture of the dynamical processes.

part by Illinois Clean Coal Institute (Grant 10/ER16).

simulations. *Phys. Rev. B* 68: 043203.

transport properties. *Chem. Phys. Lett.* 188: 315-9.

size gold clusters: Au2-Au43. *Physica E* 33: 182-190.

**6. Acknowledgement** 

**7. References** 

Fig. 11. 30-atom clusters formed under 0.5 eV initial energy with 673 K minimum temperature and a cooling rate of 1.5625 x 1011 K/s (left, c) and 1.5625 x 1013 K/s (right, d).

Figures 10 and 11 show that all four simulations predict the coalescence product is a 30 atom cluster, though they are structurally different. The faster cooling rate (Fig. 10 right and Fig. 11 right) produce a cluster that is more spreading out than the clusters produced by the slower cooling rate.

DFT calculations were performed for the structures shown in Figs. 10 and 11. The results of these clusters are given in Table 1.


Table 1. The energy difference (eV) between cluster **d** and others, HOMO-LUMO energy gap (eV), and the number of unpaired electrons of 30-atom clusters.

The DFT results in Table 1 show that structure **b** is the least stable isomer of the 30-atom clusters, which agrees with the MD simulations as depicted in Fig. 9 (right). However, the energy differences among the other three clusters are not significant in the MD simulations but are significant in the DFT calculations. More importantly, the number of unpaired electrons of the products is very different, indicating a more complex electronic state of the final product. MD simulations based on a single PES may need to be reexamined for the accuracy of the simulations.

#### **5. Conclusion**

36 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

Fig. 11. 30-atom clusters formed under 0.5 eV initial energy with 673 K minimum

slower cooling rate.

these clusters are given in Table 1.

accuracy of the simulations.

temperature and a cooling rate of 1.5625 x 1011 K/s (left, c) and 1.5625 x 1013 K/s (right, d).

Figures 10 and 11 show that all four simulations predict the coalescence product is a 30 atom cluster, though they are structurally different. The faster cooling rate (Fig. 10 right and Fig. 11 right) produce a cluster that is more spreading out than the clusters produced by the

DFT calculations were performed for the structures shown in Figs. 10 and 11. The results of

Structure Energy difference HOMO-LUMO gap Unpaired electrons

**a** (Fig.10) 0.55 0.15 56 **b** (Fig.10) 9.69 0.13 60 **c** (Fig.11) 2.37 0.11 62 **d** (Fig. 11) 0 0.06 78

Table 1. The energy difference (eV) between cluster **d** and others, HOMO-LUMO energy

The DFT results in Table 1 show that structure **b** is the least stable isomer of the 30-atom clusters, which agrees with the MD simulations as depicted in Fig. 9 (right). However, the energy differences among the other three clusters are not significant in the MD simulations but are significant in the DFT calculations. More importantly, the number of unpaired electrons of the products is very different, indicating a more complex electronic state of the final product. MD simulations based on a single PES may need to be reexamined for the

gap (eV), and the number of unpaired electrons of 30-atom clusters.

**(c) (d)** 

Meta-Molecular Dynamics (meta-MD) simulation was developed and described for studying the formation of transition metal nanoparticles. The meta-MD simulation integrates single point electronic structure calculations into the conventional molecular dynamics simulations so that instant changes of the intrinsic electromagnetic properties of the system can be monitored and obtained during the formation of nanoparticles. The results of Fe cluster formation obtained from the meta-MD simulations were presented and discussed. Additionally, the effect of cooling rates was also presented and discussed. Furthermore, using the spin-polarized DFT calculations in meta-MD simulations can also provide indications whether the electronically adiabatic treatment in the MD simulations is sufficient by monitoring the electronic state changes during the dynamic processes.

The meta-MD technique developed here should also be a good tool in studying heterogeneous catalysis by providing guidance in the design of catalysts. For instance, the detailed picture of local charge distribution may provide insight into the active site and requirement for the catalytic activity. This information will be potentially useful in the preparation of catalysts. The meta-MD simulations described here can also be employed for studying other processes where the changes of the intrinsic electromagnetic properties of the local entities of the system, i.e. subsystems, are important in order to obtain a complete picture of the dynamical processes.
