Advanced Molecular Dynamics Simulations on the Formation of Transition Metal Nanoparticles

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.


Introduction
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 ) 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 ), Cu 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. 2009;Shi et al. 2011), Au-alloys 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;Alkis et al. 2009;, 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 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.

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 m i in the Cartesian coordinates x i , y i , and z i by Here i x p, i y p , and i z 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: 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, 2005cRozas & 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 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 above equations of the Sutton-Chen potential, r ij 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 i   , is a cohesive term that describes the tendency for the atoms to stick together. The attraction between atoms is normally described by a 1/r 6 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, The sum is taken over all separations r j 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, The cohesive energy per atom is given by, Finally the bulk modulus, B f , is given by, where Ω f =(a f ) 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, 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).

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.
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, 1992bGe et al. 1992;Wang & Billing 1993;Wang et al. 1994aWang et al. , 1994bWang & 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.5625x10 11 K/s (1.3464x10 -8 eV/fs) and 1.5625x10 13 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 2004aZhang et al. 2004aZhang et al. , 2004bZhang et al. , 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. ) Single point calculations were performed based on the structural data from the MD simulations. The binding energy, 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 and discussed.

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.

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.5625x10 11 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. 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 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.  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 10 11 K/s (left) and 1.5625 x 10 13 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 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 10 11 K/s (left) and 1.5625 x 10 13 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 10 11 K/s (left) and 1.5625 x 10 13 K/s (right). 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 10 11 K/s (left) and 1.5625 x 10 13 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.

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. 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 10 11 K/s (red/short dash), 298 K and 1.5625 x 10 13 K/s (green/solid), 673 K and 1.5625 x 10 11 K/s (blue/long dash), and 673 K and 1.5625 x 10 13 (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.
(a) (b) 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 10 11 K/s (left, a) and 1.5625 x 10 13 K/s (right, b).
www.intechopen.com (c) (d) 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 10 11 K/s (left, c) and 1.5625 x 10 13 K/s (right, d). Figures 10 and 11 show that all four simulations predict the coalescence product is a 30atom 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 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.

www.intechopen.com
Advanced Molecular Dynamics Simulations on the Formation of Transition Metal Nanoparticles 37

Conclusion
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.

Acknowledgement
This work is supported by the National Science Foundation (Grant CBET-0709113) and in part by Illinois Clean Coal Institute (Grant 10/ER16). Molecular Dynamics is a two-volume compendium of the ever-growing applications of molecular dynamics simulations to solve a wider range of scientific and engineering challenges. The contents illustrate the rapid progress on molecular dynamics simulations in many fields of science and technology, such as nanotechnology, energy research, and biology, due to the advances of new dynamics theories and the extraordinary power of today's computers. This first book begins with a general description of underlying theories of molecular dynamics simulations and provides extensive coverage of molecular dynamics simulations in nanotechnology and energy. Coverage of this book includes: Recent advances of molecular dynamics theory Formation and evolution of nanoparticles of up to 106 atoms Diffusion and dissociation of gas and liquid molecules on silicon, metal, or metal organic frameworks Conductivity of ionic species in solid oxides Ion solvation in liquid mixtures Nuclear structures

How to reference
In order to correctly reference this scholarly work, feel free to copy and paste the following: Lichang Wang and George A. Hudson (2012). Advanced Molecular Dynamics Simulations on the Formation of Transition Metal Nanoparticles, Molecular Dynamics -Theoretical Developments and Applications in Nanotechnology and Energy, Prof. Lichang Wang (Ed.), ISBN: 978-953-51-0443-8, InTech, Available from: http://www.intechopen.com/books/molecular-dynamics-theoretical-developments-and-applications-innanotechnology-and-energy/advanced-molecular-dynamics-simulations-on-the-formation-of-transition-metalnanoparticles