Molecular Dynamics Simulation of Iron — A Review

Molecular dynamics (MD) is a technique of atomistic simulation which has facilitated scienti¯c discovery of interactions among particles since its advent in the late 1950s. Its merit lies in incorporating statistical mechanics to allow for examination of varying atomic con¯gurations at ¯nite temperatures. Its contributions to materials science from modeling pure metal properties to designing nanowires is also remarkable. This review paper focuses on the progress of MD in understanding the behavior of iron — in pure metal form, in alloys, and in composite nanoma-terials. It also discusses the interatomic potentials and the integration algorithms used for simulating iron in the literature. Furthermore, it reveals the current progress of MD in simulating iron by exhibiting some results in the literature. Finally, the review paper brie°y mentions the development of the hardware and software tools for such large-scale computations.


Introduction
It has already been over 50 years since Alder and Wainwright 1,2 developed molecular dynamics (MD) simulations as a computational tool used for tracing the phase space trajectory of all particles being simulated. Apart from biochemical discipline which This is an open access article published by World Scienti¯c Publishing and distributed under the terms of the Creative Commons employs MD to investigate the properties of biomolecules, materials scientists often employ MD as a step of understanding the mechanisms of physical phenomena caused by metallic atoms. This is achieved by integrating the equations of motion. Then the velocity of the particles follows the Maxwell-Boltzmann distribution, which is temperaturedependent. Accordingly the pressure acting on the particles is determined by the virial theorem. 3 Periodic boundary condition (PBC) 4 has already been employed in the early formalism of MD to avoid the surface e®ect that is common in small simulation samples. Useful physical quantities, such as the di®usion coe±cients, heat capacities and energy changes, can be determined later from the trajectories of the particles saved in the computers.
In the early days, the interatomic potentials were fairly limited to hard sphere approximations, in order to accommodate to the relatively slower calculation capability of the computers at the time MD was just developed. Driven by the demand for more complex materials, numerous interatomic potentials have been devised for a more pertinent representation of the materials as a function of interatomic separation. Initially the potentials focused on pure metals, but later on they could also re°ect the interactions and thermodynamics occurring in alloys. The approaches to formulating the potentials evolved from pure distance dependency to electronic density dependency, followed by bond order dependency.
The advantage of using MD is that one can obtain physical paths of the particles in the course of attaining thermodynamic equilibrium, which is not possible by using Monte-Carlo simulations as it can only return meaningful equilibrium values but random transient states. 2 The method of MD is mainly the solution to particle trajectories derived from the interatomic forces. Numerical integration of atomic motion is performed on the interatomic forces, which results in the particle velocity. The particle position is then obtained by further integrating the velocity. By MD approach, the phase space trajectories of the ensembles can be evaluated.
The success of MD simulation of iron relies on the proper interatomic potentials that address the particular electronic structure of iron accurately. The accuracy of MD simulation of iron is important for nuclear industry, because it can estimate the extent of damage of nuclear power plants. In practical cases, the introduction of impurities in iron potentials is crucial for investigating the e®ect of irradiation which releases a number of impurities that would interact with pure iron. Appropriate potentials of iron are also essential for estimating time evolution of defects that occur in iron, such as vacancies, interstitials, dislocations, and grain boundaries. Besides, MD simulation of iron plays a key role of understanding the e®ect of metal catalyst on the growth of carbon nanotubes.
Because of the application of cuto® distances in atomic force computation, parallel computations of forces can be applied to di®erent portions of a simulation box, with each portion having no e®ect on the other. Recon¯gurable computers and graphics processing units (GPU) can execute parallel computations, so as to accelerate the computation tasks. Science practitioners have to design the algorithms of allocating computing resources for recon¯gurable computers and GPU. The speedup of parallel computation can be over 100 compared to the sequential counterpart.
The organization of this review paper is as follows. The basic principles of MD simulation is discussed, together with a brief introduction to statistical mechanics that is directly relevant to MD formalism. The MD implementation and the corresponding algorithms are then exhibited brie°y. A number of thermostats are mentioned. Some of the techniques applied to MD simulation are provided as a supplement to the conventional MD approach. The history of interatomic potentials for iron in various formats is discussediron without spins, magnetic iron and iron with impurities. A number of categories of MD simulation for iron are exhibited, which demonstrates a wide range of applications of MD in modeling defects and nanotubes. Then the development of computer hardware used in MD simulation is discussed. A summary of the review is presented at the end.

Basics of MD Simulation
A number of references regarding the formalism of MD have been available, such as Refs. 5-10. The remaining portion of this section is a very brief summary of these references, which demonstrates the major points of interest in the MD computation technique. Before brie¯ng the MD technique, important concepts of statistical mechanics that are helpful to the development of MD are stated. The interested reader is referred to Refs. 11 and 12 for much detailed explanations.

Statistical mechanics
Thermodynamic states can be de¯ned by the set of parameters, such as number of atoms N, pressure P and temperature T. These macroscopic quantities can in principle be connected to the microscopic state of the system of interest, and statistics is such a required connection. The study of macroscopic properties via microscopic quantities is known as statistical mechanics.
A microstate of a system of particles is the basis of statistical mechanics. It represents a particular state determined by the set of phase space coordinates with some probability of occurrence. Suppose there are N particles, each with n degrees of freedom. The microstate can then be represented by a point of nN dimensions in the phase space. A particle has 3 position components fr i g and 3 velocity components fv i g, so each particle has 6 degrees of freedom. In this case, the microstate can be represented by a point s ¼ ðfr i g; fp i gÞ of 6N dimensions. The time series of s is known as the phase space trajectory Àðfr i g; fp i gÞ.
The ensemble average of an observable A, based on its probability distribution P ðfr i g; fp i gÞ, is expressed as By ergodicity principle, the ensemble average is equal to the time average as long as every point on the phase space is accessible. The time average has the form where t obs is the observation time. The ergodicity principle is very useful in MD because one can obtain the thermodynamic average from the time evolution of phase space trajectory generated by MD. A microcanonical (NVE) ensemble is commonly used when the system of interest is isolated, so that no energy exchange occurs with the surroundings. Here, the number of particles N, the volume V and the energy E are all kept constant. Each microstate has the same a priori probability. Therefore, the probability of a macrostate depends on the statistical weight ðN; V ; EÞ, which is the number of microstates of that particular macrostate. The entropy of an NVE ensemble is given by where k B is the Boltzmann constant. It is clear from Eq. (3) that the maximum entropy occurs at the maximum statistical weight. Such an equation is vital for MD because one can link the microstates of an ensemble to the thermodynamic states. Another important ensemble is the canonical (NVT) ensemble, in which the temperature T rather than the energy is conserved. In this case, energy transfer to the surroundings is permitted. The probability of occurrence of a macrostate P i follows the Boltzmann distribution: Here, E i is the energy of the macrostate, ¼ 1=k B T is the temperature parameter with Boltzmann constant k B , and Z is the partition function in the form Z ¼ P i e ÀE i , which normalizes the total probability of occurrences to 1. The entropy of an NVT ensemble is given by The average energy of an NVT ensemble, also known as the internal energy, is given by

MD principles
The idea of MD simulation is deduction of the particle motion starting from the interatomic potential. According to Newtonian mechanics, once the potential UðrÞ is given, the time-varying force F i ðtÞ of each particle i can be evaluated as F i ðtÞ ¼ mr :: Here, m is the particle mass, r N is the positions of N particles that are used to de¯ne the potential and r i is the individual particle position. Equation (7) is the basic equation governing particle motion. With atomic forces, one can perform integration to obtain the velocity and then the position after another integration of velocity. For Hamiltonian mechanics, an isolated system of particles with energy E can be expressed in terms of their positions r N and momenta p N : from which one can obtain the equations of motion as The second line of Eq. (9) is simply equal to Eq. (7) in principle. By performing integration on Eq. (9), one can also obtain the velocity and position of individual particles. Regardless of Newtonian or Hamiltonian formalism, the implementation of these integrations in MD involves di®erentiating the potential function numerically, and plugging in the interatomic distances to obtain the interatomic forces. A number of algorithms are available for the numerical integration processes. Here we mention some of them. Theoretically, the interatomic force must be calculated from the interaction of all other atoms. However, this is very time consuming. By employing a cuto® distance from an atom, one can limit the number of nearest atoms within the cuto® that are included in the force evaluation. The suitable cuto® distance should be set according to the interatomic potential, such that atomic interaction beyond the cuto® is negligible.

Verlet algorithm
After the force equations are formulated, the velocity and position of each atom can be obtained by integrating the force equations. In order to allow for numerical integration, the di®erential equations governing the motion have to be discretized in time steps Át. Accordingly, the¯nite di®erence (FD) method is commonly used in MD calculations. One type of the FD methods is the Verlet algorithm, which is derived from the di®erence of two Taylor expansions in position r: Adding them up gives Therefore, one can obtain the particle position for the next time step if one uses the acceleration aðtÞ ¼ r 0 0 ðtÞ derived from the intermolecular forces, the current position and the position for the previous time step. The advantage of using Eq. (11) to determine the position is that we do not need the atomic velocity vðtÞ ¼ r 0 ðtÞ. The velocity of particles is obtained by using¯rst-order central di®erence: The velocity depends on the position for the previous and also the next positions. The merits of the Verlet algorithm are its easy implementation and stability over large time steps.

Velocity Verlet algorithm
The velocity Verlet algorithm helps us to obtain both the velocity and position at t þ Át. The position for the next time step is simply obtained by the Taylor expression: The velocity at the next time step is evaluated by It can be seen that the evaluation of the next velocity step involves the use of the next acceleration step, which is derived from the next position step.

Leapfrog algorithm
In this method, the velocity at half time step is evaluated, which is then used to obtain the position at full time step. After the next position is evaluated, it is used to obtain the velocity for another half time step. This means velocity \leaps" over position, and position \leaps" over velocity in turn. The formulae used are The disadvantage of this algorithm is that the position and velocity cannot be evaluated at the same time step.

Predictor-corrector method
This approach is a three-step process. In the¯rst step, the velocity and position for the next time step are predicted. The acceleration at the next time step is evaluated by the predicted velocity and position. In the¯nal step, the initially predicted velocity and position are corrected with the evaluated acceleration. By modeling particle interaction as produced by harmonic oscillators, we predict the next velocity and position as Here, ! is the angular frequency of an oscillator. Then we evaluate the acceleration for the next time step as After the acceleration is updated, the predicted values are corrected by using those for the next time step:

Gear's predictor-corrector method
This is an improved version of the original predictorcorrector method obtained by employing the¯fthorder Taylor expansion. Therefore, the particle position for the next time step is predicted in terms of ve derivatives: The interatomic forces are evaluated by using the predicted positions. The force is given by where r ij (t) is the interatomic separation, andr ij (t) is the unit vector of the interatomic separation. From the evaluated forces, one can¯nd the di®erence Á r :: between the predicted and evaluated acceleration, such that Ár ð2Þ ¼ ½r where quantities with superscript P is the predicted value for the next time step. The correction would become The values of 's are¯ne-tuned to ensure numerical stability. The 's are determined by the order of the di®erential equations and the order of the predicted Taylor expansion.

Thermostats
Simulation of NVT ensembles requires the application of a thermostat that maintains the ensemble at constant temperature. There are a number of implementations of such thermostats.

Anderson thermostat
The coupling of the Anderson thermostat to an NVT ensemble is achieved by introducing stochastic collision forces that act occasionally on randomly selected particles, such that the particle forces of some atoms are altered for just a short time instant. The frequency of stochastic collision represents the coupling strength to the thermostat, having a Poisson distribution of where P ðt; Þdt is the probability of the next collision at time t þ Át. The motion integration is divided into three steps. First, we initialize the positions r N and momenta p N of N particles, and perform motion integration up to the instant before the¯rst stochastic collision. Second, some particles are randomly chosen to have a collision with the thermostat. Third, the momentum of the particles after the collision is chosen from the Boltzmann distribution at the desired temperature T. All other particles are una®ected.

Nos e-Hoover thermostat
It is an extension to the conventional Lagrangian form by introducing one more coordinate s, such that Here, Q is the e®ective mass associated with s, and g is number of degrees of freedom of the system. The momenta conjugate to r i and s are The Hamiltonian form of the system can be expressed as The Hamiltonian in Eq. (26) leads to the following equations of motion: The extended microcanonical ensemble has 6N þ 2 degrees of freedom, with partition function If we set p 0 ¼ p=s and r 0 ¼ r, as real variables, then If we choose g ¼ 3N þ 1, then where In this condition, the ensemble average of an observable A follows the relation This means that the extended system in real variables can reduce to a canonical ensemble. Also, by letting s 0 ¼ s and t 0 ¼ t=s as other real variables, we can transform the equations of motion in Eq. (27) as

Velocity scaling
This is a very straightforward method by scaling the particle velocity v i by , where Here, T 0 is the target temperature. The disadvantage of this approach is that the result does not correspond to a canonical ensemble. The momentum space generated by this method results in discontinuity.

Berendsen thermostat
Unlike the simple velocity scaling approach that modi¯es the velocity in one step, Berendsen thermostat does the scaling slightly for each time step. The rate of temperature increase is expressed in a di®erential equation where is the coupling strength between the system and the thermostat. The change in temperature ÁT can be expressed as The particle veloctiy v i of each particle can then be scaled to , where Here, T is the time constant characterizing the rate of achieving the target temperature T 0 .

Langevin theromostat
A stochastic approach to maintaining temperature is called the Langevin thermostat, by which a time-varying random force ðtÞ following Gaussian distribution is introduced to the equation of damped motion, such that Here, is the damping constant. The random force satis¯es delta-correlation, such that with being a constant characterizing the strength of the random force. The idea of this thermostat is the choice of that can achieve the target temperature. At this time, the damping force will balance the random force.

Periodic boundary conditions
A simulation box con¯nes the region where particles can be located. Yet, the simulation results generated from this box can fail to represent the bulk condition because of surface e®ect that occurs at the boundary planes con¯ning the box. An approach to correct this problem is the introduction of identical copies of simulation boxes contiguous with the original simulation box. Motion integration has to incorporate the wraparound e®ect when a particle leaves the original simulation box. For example, the x-coordinate of a particle is bounded by ÀL x =2 x L x =2, where L x is the length of the simulation box. If the particle position r x i ! L x =2, then replace the position by r x i À L x . Similarly, if r x i ÀL x =2 then replace the position by r x i þ L x . The same treatment has to be made on interatomic separation when the interatomic potential is updated.

Techniques of Applying MD
In addition to the MD formalism, many techniques have been developed to enrich its applicability to a number of re¯ned situations.

Spin-lattice dynamics
Spin-lattice dynamics (SLD) is a modi¯ed approach of MD, in order to incorporate both spin and lattice degrees of freedom in a Hamiltonian. 13,14 In this formalism, the spin and lattice degrees of freedom are coupled by the exchange integral term, in the sense that the lattice degree of freedom would change the behavior of the spin degree of freedom, and vice versa. The SLD formalism is thus suited for spin-carrying materials, such as iron.
The corresponding Hamiltonian is given by which has four components: lattice kinetic energy, lattice potential energy, magnetic energy as a result of spin-lattice coupling and the magnetic energy due to an external¯eld H ext , respectively. In Eq. (40), m i is the mass of atom i, fp i g is the momentum space, fR i g is the lattice space, fe i g is the classical spin space of unit length, S is the spin vector length. Also, j ij ðR ij Þ represents spin-lattice coupling, which is the product of the exchange integral J ij ðR ij Þ between spins i and j and the norms of the spins, such that j ij ðR ij Þ ¼ S i S j J ij ðR ij Þ, and UðfR i gÞ is the total lattice potential. Physically, e i Á e j signi¯es the spin-spin correlation. The constant g is the gyromagnetic ratio, and B is the Bohr magneton. For the de¯nition here, the direction of the magnetic moment is opposite to that of the classical spin, such that M i ¼ Àg B S i . It is noted that this format of Hamiltonian is isotropic. The equations of motion for the momentum, lattice and spin components can be derived from the time derivatives of Eq. (40), returning In Eq. (43), the e®ective magnetic¯eld H eff i is given by The equations of motion are then implemented using conventional MD approach, except that the spin motion has to be evaluated separately. An application of SLD is exhibited in Ref. 15, which refers to modeling the iron thin¯lm behavior. With the SLD formalism, the thin¯lm magnetization decreases with temperature, having roughly the same temperature dependence as a bulk demonstrates. The thin¯lm temperature dependence is also found to vary with the¯lm thickness. The magnetic transition temperature also decreases with¯lm thickness. The surface magnetization is also di®erent from that inside the bulk. It is also noticed that the introduction of spin-lattice coupling and spin-spin correlation can result in the near-surface relaxation strain that varies with temperature.

Thermodynamic integration
Thermodynamic integration (TI) is a computational approach to evaluating the free energy di®erence between two states. The description below mainly follows Ref. 8. It is known that the Helmholtz free energy is equal to where Q is the partition function in the form Here, Ã ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h 2 =ð2mk B T Þ p is the thermal de Broglie wavelength, with h being the Planck's constant and m being the particle mass, which is not related to the canonical average over the phase space but to phase space volume that is accessible to the system. The Helmholtz free energy change cannot be measured directly from real or computer experiments, because it depends on the partition function which cannot be evaluated numerically.
The idea of TI is the coupling of two thermodynamic states with reference Hamiltonian H I and target H II by a switching parameter . An intermediate state between H I and H II is given by a thermodynamic path The free energy di®erence ÁF between two thermodynamic states characterized by ¼ 0 and ¼ 1 is given by The angle brackets hÁ Á Á i represent taking the ensemble average over . The linear path in Eq. (47) is a convenient choice because So, @F =@ decreases with increasing . The Frenkel-Ladd method 16 of TI is often applied to solid phase governed by a hard-sphere potential. The idea is the construction of a thermodynamic path from the system of interest to a noninteracting Einstein solid having the same structure as the required system. Here, a noninteracting Einstein solid consists of noninteracting atoms coupled to the lattice sites by harmonic springs. Since we cannot switch on the springs while switching o® the hard-sphere potential, the thermodynamic path is modi¯ed to the form Here, H 0 is the unperturbed Hamiltonian, N is the number of atoms, r i0 is the lattice position of atom i and r i is the atom position of atom i. The free energy di®erence between the Einstein solid and the system of interest is then given by The system reduces to an Einstein solid for large value of ! max . Later, the Frenkel-Ladd method is modi¯ed by using an \ideal Einstein molecule", which is a noninteracting Einstein solid except one atom that is not coupled to a harmonic spring. 17 Let the Helmholtz free energy of the ideal Einstein molecule be A ideal . Then we¯x the position of the atom that is not coupled to a harmonic spring (called atom 1). In this situation, the thermodynamic path connecting a hard-sphere potential and a set of harmonic springs is given by where H 0 is the hard-sphere potential, r 1 is the position of the¯xed atom, r i is the position of atom i and r i0 is the lattice position. The Helmholtz free energy due to the harmonic springs is given by whereas the Helmholtz free energy due to the hardsphere potential is given by Here, u 0 ðr ij Þ is the hard-sphere potential depending on the interatomic separation r ij . The free energy of a solid is then found to be By¯xing one atom in¯nding free energy di®erence, the whole lattice needs not be totally¯xed as is necessary for the Frenkel-Ladd method. The implementation of the TI approach is thus more straightforward.
The method of TI has been applied to study the vacancy formation and migration occurring in BCC iron, 18,19 whose Hamiltonian is modeled by the SLD formalism. 13,14 The study uses the SLD formalism to determine the free energy change in the course of vacancy formation and migration by TI. By using both magnetic and nonmagnetic potentials, it is found that a vacancy leads to scattering of magnons, leading to an increase in the total free energy. It is also noted that the magnon state determines the exchange interaction and hence the interatomic force. The temperature dependence of magnon distribution is crucial for the lattice properties as well. The phonon-magnon interaction would lower the energy barrier and increase the entropy of vacancy migration and formation.

Interatomic Potentials
Iron is a commonly found in pure metal form or in alloy form, bringing about various applications in industry. Accordingly, computational tools have to be devised that can return reliable information about iron and its other structures. In fact, it is best to determine the metallic properties by ab initio computations, which can evaluate the interactions among atoms in electronic level. However, the related process is highly computationally demanding and time consuming. Instead of ab initio computations, MD is an e®ective tool of simulation in materials science and engineering, given a potential that describes the atomic interactions pertinently. Then interatomic potentials become to play an important role of determining the time evolution of defects. In order to represent the conditions of iron atoms in various states, a number of interatomic potentials have been devised. Later on, they have been further modi¯ed to suit to practical situations.
In¯tting a potential empirically, many physical quantities are derived from it as a veri¯cation against the experimental values. The remainder of this section reviews several potentials commonly used in modeling iron or iron alloys. Their adapted forms, if any, are also discussed to investigate the improvement in describing more complicated physical phenomena. The veri¯cations against experimental results made by the developers of the potentials are also mentioned.

Pure iron potentials
The embedded-atom method (EAM) 20,21 is an approach based on the density functional theory (DFT), used to model the ground state properties of FCC metals with impurities. This method is an improvement of the previously developed pair potential 22 that requires accurate volume-dependent energy to describe the elastic properties, which can sometimes be ambiguous in situations involving surface defects and cracks. The initial use of EAM is to model the role of hydrogen atoms in steel, which leads to brittle fracture and cracks. It considers the pair potential plus the energy required to \embed" an atom in an electronic density constituted by a host lattice. The energy functional of a system of atoms is then expressed as In Eq. (57), F i is the embedding energy of an atom, to be determined by the local electronic density i ðR i Þ at position R i but without atom i, and V ðR ij Þ is the short-range electrostatic pair potential of repulsive nature due to the neighbors of atom i. The values of F i and V ðR ij Þ are based on the experimental results of the atoms concerned, such as the lattice constants, elastic constants and the migration energy of the impurities. Usually, a monotonically decreasing form is used for i ðR i Þ. Since the electronic density is clearly de¯ned without considering the volume, the ambiguity arising from conventional pair potentials is resolved.
At the time of developing the EAM, this method is not guaranteed to be the universal form of modeling transition metals, especially for BCC metals. Accordingly, many adaptations have been made as an attempt to extend the use of EAM. One can¯nd many methods of implementing the EAM for various transition metals in the literature. This review highlights those relevant to iron. For example, Johnson and Oh 23,24 have extended the application of EAM to determination of the short-range potential for alloys and BCC metals including iron, respectively. The potentials are¯tted to the¯rst and second nearest neighbors for both alloys and pure BCC metals. The idea of¯tting is by means of transformation of an EAM potential to a normalized form, such that the¯rst derivative of the embedding function with respect to the electronic density is zero, i.e., dF i =d i ¼ 0. With this normalization, the potential turns to an e®ective pair potential that is only slightly dependent on the embedding function. The determination of the potential can become conceptually easier, while maintaining the characteristics of an EAM. Wang and Broecker 25 have extended the interatomic distance of the pair potential to the¯fth nearest neighbors by incorporating a Gaussian function as the weighting factor and by using an oscillatory model in terms of sinusoidal functions for the local electronic density. These modi¯cations suggest that the anomaly of the phonon spectra of BCC transition metals can be better reproduced.
Based on the EAM formalism, Mendelev and coworkers 26 have¯tted an EAM potential for crystalline and liquid iron. In fact, this potential forms the basis of some upcoming potentials of Fe alloys, to be discussed later. They established three potentials, each by¯tting to asymmetric crystal defect data by considering atomic force of liquid iron obtained by ab initio calculations, to experimental structure factor of liquid iron at melting point, as well as to symmetric perfect crystal data by the EAM approach, respectively. In this sense the state with small interatomic separation can also be accommodated, and the solid-liquid phase transition of iron can also be described more accurately.
Later, many interatomic potentials of iron have evolved from the work by Mendelev et al. 26 Chamati et al. 27 have followed the EAM approach to develop a potential which is more suitable for BCC iron (-Fe) as well as FCC iron (-Fe). The strength of this potential is that it can reproduce various BCC and FCC parameters such as thermal expansion coe±cient, phonon dispersion relations, meansquare displacements and surface relaxations with-out¯tting with the corresponding experimental results. Other examples of the EAM potential for iron in various states can be found in Refs. 28-30. The modi¯ed embedded-atom method (MEAM) [31][32][33] has been developed for computational simplicity and hence wide applicability to MD simulations. The development responded to the high demand of potentials suitable to model semiconductor physics, in which simple pair potentials are not able to reproduce the elastic constants of covalent structures accurately. In fact, MEAM has also found its values in modeling metallic structures.
The merit of MEAM lies in the addition of angular dependence of the bond-bending forces that constitute the background electronic density, as opposed to linear superposition of radially-dependent atomic densities for the original EAM approach. The main idea of the modi¯cation comes from the host electronic density, which includes angular correction terms in addition to the simple pair potential form. The general expression of the total energy is similar to that of EAM, known as The de¯nition of the symbols follows that of EAM. The change for the modi¯ed EAM lies in the evaluation of F and V ðR ij Þ. The embedding energy F is known as where E 0 is the sublimation energy, and " ¼ n 1 a ðr 1 Þ is the background electronic density due to n 1¯r st nearest neighbor atoms of the reference structure in the form of a monatomic homogeneous solid, with each of the neighbors having equilibrium rst-neighbor distance r 1 . The format of Eq. (59) corresponds to the logarithmic relationship between bond length and number of bonds. The atomic electron density a ðr 1 Þ is exponentially decaying with r 1 : The value of is dependent on the atomic density in units of distance À1 . In Eq. (58), the local electronic density i experienced by atom i exhibits angular dependence by where jik is the included angle formed by atoms j, i and k, a is the¯tting constant determined by shear modulus data. The pair potential V ðR ij Þ in Eq. (58) is given by ÈðrÞ ¼ X s n s n 1 V a s r ð Þ: The notations in Eq. (62) are de¯ned as follows. r is the¯rst-neighbor distance, n s is the number of s-nearest neighbor atoms, a s is the ratio of the s-neighbor distance to the¯rst-neighbor distance, r c is the cuto® distance.
In the early formalism, only the¯rst nearest neighbors were e®ective to calculate the potential, on the basis of DFT computations on BCC tungsten (see Appendix of Ref. 33). The MEAM has been further re¯ned to incorporate the atoms in the second nearest neighbors, 34,35 known as 2NN-MEAM.
The MEAM approach has also been applied by many workers in order to derive more alloy potentials that have been¯tted according to more abundant experimental or ab initio¯ndings. Here we introduce some of the MEAM applied to various situations. Yuan et al. 36 set up a new scheme of determining the embedding energy of BCC transition metals by slightly modifying the embedding energy in Eq. (59) with a¯tting parameter , such that F ðÞ ¼ E 0 ð=" Þ lnð=" À Þ. With , the MEAM potential can be¯tted even to nonbulk systems such as surfaces. The potential returns the crystal elastic sti®ness, the vacancy formation energy and the low-index surface energies, which are close to experimental¯ndings. Jelinek et al. 37 have performed a large-scale formulation of MEAM potentials for Al, Si, Mg, Cu and Fe alloys, which results in improved values of the generalized stacking fault energy curves. The resulting potentials have been validated by comparing them to corresponding DFT results, together with a number of properties such as equilibrium volume, elastic constants and defect formation energy. More examples of the MEAM potentials can be found in Ref. 38. In addition, a detailed review of the MEAM potentials can be found in Ref. 39.
An approach that evolves from EAM to model the potential of transition metals is the Finnis-Sinclair (FS) potential. 40 It aims at rectifying the drawbacks of adopting pure pair potentials in modeling metallic defects such as dislocations and grain boundaries, as a result of the unsatisfactory consideration of elastic constants, Cauchy pressure and vacancy formation energy. The FS potential is an empirical one that solves the problem of pair potentials by considering the metallic cohesion of an atom due to its neighbors according to the tightbinding theory of metals, 41 based on the second moment approximation, as well as the repulsive pair potential. Its initial form focuses on the BCC structure. The general form of the FS potential is expressed as In Eq. (63), the total potential energy U is the sum of the repulsive pair potential U P and the N-body cohesive energy U N . So U P considers the sum of the pair potential V FS ðR ij Þ dependent on the interatomic separation R ij . The cohesive energy is the sum of the square root of the electronic density p i , multiplied by a positive proportionality constant A. The electronic density of each atom is the sum of the cohesion energy ðR ij Þ due to the neighboring atoms. It can be observed that the FS potential simply uses a predetermined square root form of the embedding energy, while the EAM form requires an embedding energy to be ascertained by¯tting.
The parameterized pair potential and cohesive energy of the FS potential for iron is expressed as where c and d are the cuto® distances of V and , respectively, between the second and third nearest neighbors of a BCC structure, and B is a value such that will achieve its maximum within the¯rstnearest neighbor distance. Soon after the FS potential, Ackland and Thetford 42 have improved it by introducing a core term in the pair potential V ðR ij Þ, which increases the short-range force repulsion that is lacking in the original FS potential. With this correction, the atoms at short separation would not fall together. The expression of V ðR ij Þ is recast as Here, V FS ðR ij Þ is the pair potential in the original FS formalism, and b 0 is the nearest-neighbor distance. B and are¯tting parameters. After varying the pair potential term, the pressure-volume relationship of many BCC transition metals becomes more physical. In addition, the altered potential has been veri¯ed by checking the formation energy of an interstitial against those in other theoretical studies. The original FS potential has considered up to the second nearest neighbors, and hence it can be regarded as a short-ranged one. The expressions of the potential are¯tted against some material parameters, such as the lattice constant, elastic constants, bulk modulus and Cauchy pressure of various BCC structures. Initial development of the FS potential has been veri¯ed to give rise to a stable BCC structure, yet its applicability to FCC and HCP structures is questionable due to its short-ranged¯tting strategy.
The FS potential has been adapted since its release in order to extend its range of applications. For example, Dai et al. 43 have attempted to extend the FS potential to FCC transition metals, and corrected the BCC potential to follow the Rose equation of state of metals 44 with increased accuracy. This is achieved by applying a sextic expression to the original quartic form of the repulsive pair potential, and by applying a quartic expression to the original quadratic form of the cohesion energy. So some of the expressions in Eq. (64) have changed to where c and d are the cuto® distances of V and , respectively, between the second and third nearest neighbors of a BCC or FCC structure one is investigating. The extended format can be reduced to the original FS potential when c 3 , c 4 and B are set to zero. This format has found good agreement with the experimental lattice constants, cohesive energies, elastic constants and vacancy formation energies. With this extension, the lattice beyond the equilibrium state can be better represented. Apart from the successful extension to FCC structure, the pressure-volume relationship of BCC structures derived from the extended FS potential becomes more satisfactory after checking with the equation of state. The extended FS potential can even¯nd good agreement with various FCC-BCC cross potential determined by ab initio calculations.
Based on MEAM which considers angular dependence of bonds, Müller et al. 45 have formulated a potential that considers the analytic bond order explicitly, such that the lack of conventional MD in describing electronic degrees of freedom can be remedied. This potential has been successful to model À phase transition and -iron before the melting point.

Magnetic iron potential
Though the EAM and MEAM formalism have been successful in deriving a number of potentials that have found practical values in modeling irradiation damage of steels, the magnetic e®ects of iron has not been considered explicitly in these potentials. In order to¯ll this research gap, potentials speci¯c to magnetic iron have been developed with further modi¯cations after its initial version. Dudarev and Derlet 46 have made a¯rst attempt on -iron potential that can describe both the magnetic and nonmagnetic states of iron. This potential (named as DD potential) applies the EAM formalism with in-depth evaluation of the local electronic density based on Ginzburg-Landau (GL) model and Stoner model up to the second-moment description of the electronic density of states. The GL model describes the second-order phase transition of ferromagnetic iron, whereas the Stoner model describes the correlation e®ect that leads to band magnetism used to characterize the ground state of ferromagnetic 3d transition metals. It is then found that the symmetry-broken solution of the GL model is able to link magnetism and interatomic forces. After applying all these models, the embedding function for both magnetic and nonmagnetic states of Fe: Here, the¯rst term on the right-hand side of Eq. (67) is the FS embedding potential corresponding to the cohesive energy, with A being a constant to be determined. The second term on the right-hand side of Eq. (67) is the magnetic potential term due to Stoner and GL model. In this term, B and are constant, Â is the Heaviside unit step function, and c is the critical electronic density beyond which the magnetic e®ect vanishes. However, there is a problem with this embedding potential: the derivatives are discontinuous at ¼ c .
In view of this shortcoming, Eq. (67) is further modi¯ed as Parameterization takes place on both F ð i Þ and V ðR ij Þ of Eq. (57), with i being written in the Then fðR ij Þ and V ðR ij Þ are written in cubic knot functions: In Eqs. (69) and (70), f n , V n , r f n and r V n are determined by¯tting, given N f and N V terms to include in the knot functions, respectively. The initial two versions of DD potential were¯tted to bulk BCC magnetic properties, vacancy formation energy, isotropic BCC bulk of nonmagnetic iron and interstitial energies. The interested reader is referred to Ref. 47 for an attempt to use the DD potential in large-scale MD simulations containing around one million atoms, which helps to validate the applicability of the potential in evaluating the magnetic moments around an interstitial defect and in describing the migration of self-interstitial and multiple-interstitial con¯gurations in iron. Since this version, some modi¯cations of the DD potential have been achieved by¯tting it to more ab initio data such as the third-order elastic constants and the conditions during the 1=2h111i screw dislocation. 48,49

Potentials for iron with impurities
Practical interatomic potentials should be able to re°ect the interaction with impurities, which is a typical consideration of modeling the materials used for nuclear fusion reactors. For example, appropriate potentials to model the Fe-C alloy are necessary to produce steel. Defects of steel have to be modeled by an appropriate interatomic potential.
Ackland et al. 50 have developed a potential (known as AMS potential) for phosphorus in -iron, through comprehensive analysis of ab initio and experimental results of iron defects, as a tool of investigating irradiation damage caused by phosphorus which shifts the ductile-to-brittle transition temperature of steels. For such a potential to be successful, the P-P interaction, Fe-Fe interaction and Fe-P interaction have to be considered collectively in the evaluation of the pair potentials and the pair electronic densities. The formulation relies on the pure Fe potential in Ref. 26 in view of its applicability to point defect interactions. In addition, large-scale ab initio computations have been performed as the basis of the¯tting process, including Fe monolayers and surfaces, substitutional impurities, vacancies and liquid state. Only those values matching with the experimental results are included in the¯tting process. By using this potential, vacancy and interstitial mechanisms of P atom in Fe matrix can be realized from MD simulations.
As another example, if the e®ect of helium gas on the fusion reactor materials is considered, a potential should be formulated for the Fe-He alloy to better model the irradiation damage which causes void swelling, helium bubbles and blistering. These mechanical defects are often investigated by MD simulations, so a proper interatomic potential is required to model Fe-He materials. Many attempts have been made to develop such a potential. For instance, Seletskaia et al. 51 have formulated a Fe-He potential based on the electronic structure calculations. It consists of a repulsive pair potential and a three-body embedding term. This potential was¯tted with the ab initio computations of the formation and relaxation energy of He defects and He clusters, together with the AMS potential. 50 The detailed implementation of the three-body term is elaborated in Ref. 52. It also considered highly the interstitial properties and applied a three-body potential term to improve the¯tting. Juslin and Nordlund 53 have later developed a Fe-He pair potential based on the C. P. Chui et al. AMS potential 50 to model helium atoms in iron matrices, which was found to be already su±cient to reproduce simple defects of iron due to helium irradiation. Since experimental data for Fe-He clusters were lacking at the time of developing this potential, the e®ectiveness of the potential in modeling migration barriers of helium in iron was realized by verifying it against the DD potential, FS potential and DFT computations. Later, another Fe-He potential based on the multiple lattice inversiontechnique was proposed to solve the problem of tting a potential which requires multiple parameters. 54 Its applicability was established after it was used to reproduce the elastic constant, binding energy and migration barrier of Fe-He crystals obtained from other similar potentials. Gao et al. 55 have developed a Fe-He potential based on the s-band model to describe the many-body interaction, together with the embedding form and repulsive pair potential. To verify it, the binding energy required for an additional He atom to approach a He cluster, together with the migration energy of an He cluster in -iron were reproduced, with fairly good agreement with the ab initio results. Another pair potential for Fe-He materials has been formulated not only by adjusting the method in Ref. 56, but also by¯tting the magnetic potential formalism based on that used in Refs. 46-49. It is found that the values of formation and migration energies of He atom in Fe agree with the ab initio results.
With an abundant choice of Fe-He potentials, the applicability of each potential to a speci¯c physical situation has to be examined carefully. In order to ease this examination process, an interatomic potential design map has been developed for Fe-He potentials, 57 from which one can assess the uncertainties of using a certain potential in modeling a particular defect. It is expected that design maps of this type can be extended to other types of interatomic potentials, thereby facilitating the sci-enti¯c community.
A suitable potential for Fe-Cu binary alloy is important in modeling Cu precipitates, which could lead to embrittlement of reactor pressure vessels. 58 A number of Fe-Cu potentials have been available to investigate this irradiation damage. For example, the Fe-Cu alloy potential has been developed 59 on the basis of 2NN-MEAM, by combining the MEAM potentials for pure Fe and Cu, respectively. Thē tting of this potential is done to reproduce the lattice constants in Fe-rich BCC and Cu-rich FCC phases, enthalpy of the liquid mixture, and the binding energy of a Cu atom in BCC Fe matrix. As another example, Pasianot and Malerba 60 have later developed a binary-alloy potential for Fe-Cu based on EAM, by incorporating the phase diagram data of Fe-Cu systems, such that the thermodynamic functions of the systems can be re°ected in the potential, and that the radiation damage can be modeled with higher accuracy. Other attempts of Fe-Cu potential can be found in Refs. 61-63. The interaction of hydrogen atom with iron is another concern of materials scientists, because it is related to the irradiation damage of steel in nuclear plants and to the physical conditions of containers used to store or transport hydrogen as a source of clean energy. Accordingly, some Fe-H potentials have been designed. For example, Ruda et al. 64 have performed a detailed exposition of the EAM potentials of hydrogen in various metals including iron. Potentials speci¯c for pure metals have been adapted to the determination of the metal-hydrogen pair potentials. Thermodynamic heat of solution of H and lattice expansion in the course of H dissolution form the basis of parameter¯tting. As another example, Lee and Jang 65 have formulated a potential for Fe-H system by means of 2NN-MEAM, with tting parameters coming from experimental parameters such as the dilute heat of solution of H in Fe and the binding energy of H in Fe. With this potential, the role of H atom in vacancies, dislocations and grain boundaries can be predicted. Ramasubramaniam et al. 66 have developed an EAM potential for Fe-H system by adapting the pure Fe potential formulated by Mendelev and coworkers. 26 This form of potential can inherit the property of the Mendelev potential 26 in modeling screw dislocations with comparable accuracy to corresponding DFT computations. The physical quantities derived from this potential can model the di®usion and di®usion of H in -Fe, binding of H to free Fe surfaces, and the trapping of H atoms in defects well.
Carbon is an important impurity of iron, because its introduction increases the tensile strength of iron. The dislocation movement in iron can be impeded by carbon impurities. A number of interatomic potentials for Fe-C alloys have been developed.
Here we illustrate a few of them. An EAM potential of Fe-C alloy has been formulated by¯tting experimental and ab initio data to an e®ective pair interaction. 67 The equilibrium lattice constant, the bulk modulus and cohesive energy in stable and metastable states have been adopted as the¯tting parameters. The potential has been tested by checking against martensite transformation, C interstitials in Fe grain boundaries and C interstitials in a free surface of Fe-C alloys. In order to model the point defects of Fe-C alloys, which is not the strength of the aforementioned potentials, a potential has been made 68 as a remedy. The potential is based on FS formalism, with incorporating C-C potential used to describe defects containing more than 1 C atom. This potential can cater for arbitrary point defect concentration. A number of formation energies are used as the¯tting parameters: carbon interstitials in a perfect BCC Fe lattice, 1C-1V clusters, 2C interstitials in a perfect BCC Fe lattice, 2C-1V clusters and Fe 3 C. Another Fe-C potential used for designing carbon nanotubes from the carbon-saturated metal clusters has been developed in the bond-order formalism, whose quantities used for¯tting are derived from DFT. 69 For example, energies of symmetrical Fe-C clusters for varying bond lengths have been obtained from DFT. The energies of isolated C and Fe atoms of varying bond length have been obtained as well.
Some other Fe potentials with other impurities are brie°y mentioned. Besson and Morillo 70 have developed a potential for B2 and D03 Fe-Al alloys by EAM and pair potential formalism, veri¯ed by checking it against the elastic constants, in order to better study the interfacial properties during segregation of grain boundaries. After that, Lee and Lee 71 have developed a more practical potential for Fe-Al alloys by 2NN-MEAM, such that the Fe-Mn-Al-C system commonly found in steel is considered as a whole. The structural, thermodynamic and elastic properties of Fe-Al binary alloy could be modeled successfully by this approach. Another MEAM potential for Fe-C system has been formulated by using MEAM Fe and C potentials, 72 with intricacy lying in the comprehensive consideration of a carbon atom in various interstitial con¯gurations inside a Fe matrix. After this consideration, the dilute heat of solution of carbon, the vacancy-carbon binding energy and the migration energy of carbon atom in Fe matrix can be reproduced with high degree of experimental agreement. A potential used for modeling high-nitrogen steel has been developed 73 by means of 2NN-MEAM potentials of pure Fe and N, respectively. It is found that nitrogen in iron can result in improved ordering tendency in BCC and FCC iron, compared to carbon in iron, due to the

Simulation Results in the Literature
MD simulation is often regarded as a replacement of real experiments that are technically formidable and less controllable. Some assumptions may have been adopted to improve the computational speed, nevertheless the generated results can su±ciently reveal the physical behavior of atomic interactions. This section exhibits some of the important MD results of iron properties by means of MD calculations, demonstrating the practical value of this widely used technique.

Phase transition
Studies of phase transition or transformation occurring in iron is a hot topic at the time of writing this review, in the sense that martensitic transformation in nuclear power plants and shape-memory materials can be better studied. Determining the temperature at which phase transition occurs is also a major target of the study. Many MD simulations have indicated that the Bain transformation path 80 has been followed in the course of FCC-BCC transformation.
For example, the transformation from FCC to BCC has been demonstrated by MD simulation, 62 showing that FCC-Fe can be transformed perfectly at temperatures of 100-1800 K, under pressures of 0-40 GPa. Such a transformation largely follows the Bain transformation path. Another extensive MD simulation of FCC-BCC transformation has been undertaken on pure iron, in which both the Nishiyama-Wasserman (N-W) and Kurdjumov-Sachs (K-S) orientations of FCC-BCC interfaces have been attempted. 81 Figure 1(a) shows the gradual propagation of FCC structure to BCC structure at 1200 K, near the phase transition temperature of 1185 K between FCC and BCC phase. The growth of FCC structure is mainly planar. The ledge structure is developed at the cross section of the FCC-BCC interface, as can be seen from Fig. 1(b). If K-S interface is adopted, the time evolution of transformation is the one shown on Fig. 1(c). The C. P. Chui et al. arrows in the sub¯gure indicate the needle-like growth of the interface. It is found that the interfacial atoms rearrange themselves, following the Bain transformation path, to reduce the mismatch in the course of phase transformation. In view of the experimental di±culty in capturing the interfacial motion during phase transition, MD simulation has been adopted to explain the corresponding behavior in atomic layer level. 82 It is found that FCC interfaces require some temperature-dependent incubation time to e®ect on a few layers before they undergo very fast transformation to BCC structure. A certain structure similar to screw dislocations has to be established during the incubation time, after which the interface can move quickly. It is realized that the volume to surface area ratio is decisive of the incubation time but not decisive of the transformation rate. The FCC-BCC transformation in iron thin¯lms through the direct and inverse Bain path has been simulated, 83 whose¯ndings are supported by the corresponding variation with the elastic moduli. The correlation between the¯lm thickness and elastic moduli has been identi¯ed. It is found that the change in biaxial strains is responsible for the transformation mechanism. The atomic con¯guration during FCC-BCC transformation has also been investigated by Engin and Urbassek. 84 By using the FS potential, 40  The transition between BCC and HCP phase has also been investigated in the literature. For example, it is found that, by MD simulations in Lagrangian form, 86 the uniaxial tensile stress can induce such a transition. 87 With this approach, the volume can vary with the internal and external stress. The induction of structural change from BCC to HCP due to a uniaxial stress can be realized, together with the intermediate change of structure due to asymmetric shear deformation. The shear deformation becomes more uniform at the end of transformation. The reverse transformation from HCP to BCC can be undertaken by a uniaxial compression, at the expense of a hysteresis loop. The reverse transformation is undertaken by pure homogeneous shear mechanism, which is di®erent from that of the direct transformation. The symmetry breaking mechanism might be responsible for such a di®erence. Morris and Ho 88 have made a step further by analyzing the structure factor in the course of BCC-HCP transformation, suggesting that the Brillouin-zone dependence of scattering is greatest in the course of transformation, indicated by the formation of Bragg peaks that is responsible for the HCP structure. The e®ect of directional loading on transformation between BCC and HCP/FCC is discussed in Refs. 89 and 90, with corresponding time evolution of atomic con¯gurations. For a better representation of both forward and backward phase transition between FCC/HCP and BCC iron states, a speci¯c potential has been formulated. 91 The temperature ranges in which the FCC and HCP phase of iron is unstable have been determined by MD simulations using this potential. The FCC-BCC transformation is found to follow Bain path, 80 while the BCC-FCC transformation is found to follow a Burgers mechanism. 92 Transformation between FCC/HCP and BCC states in iron nanowires has been simulated, 93 which indicates that tensile axial stress can vary the phase transition temperature of iron nanowires and that the transition temperature has an inverse relation with the wire diameter. However, the application of stress beyond a critical value can inhibit transformation from FCC/HCP to BCC phase. Hysteresis e®ect has also been observed in the e®ect of temperature on nanowire length.
MD simulations for characterizing Fe alloys have been attempted, so that the thermodynamic properties of them can be realized. For example, Yang et al. 94 have performed MD simulations on undercooled Fe-Ni alloy in liquid state to¯nd its heat capacity because the time required for measuring the heat capacity already allows the alloy to crystallize, leaving the phase of interest. They have employed the EAM implemented with analytic nearest-neighbor interactions 23,24 to characterize the potentials of both Fe and Ni. Then they have determined the heat capacity of Fe-Ni alloy by di®erentiating the energy with respect to temperature, and have concluded that the composition of an alloy is deterministic of the heat capacity at the undercooled state. Kadau et al. 95 have investigated the phase transition occurring in Fe-Ni nanoparticles, from which the scaling behavior of transition temperature with the inverse of particle diameter is observed. Besides, the N eel temperature of FCC Fe is found to decrease with cluster size.

Interstitials, dumbbells and crowdions
As the atoms in a crystal structure outnumbers the lattice sites, the extra atoms may need to occupy a space that is not reserved for atoms. These extra atoms can be the same as the lattice atoms (known as self-interstitials), or di®erent from the lattice atoms (known as impurity interstitials). Many of the con¯gurations pertaining to interstitials are possible. For example, an atom originally¯xed at a lattice site may be substituted by another atom that has no¯xed site (known as substitutional atom).
The study of interstitials usually involve the time evolution of one self-interstitial atom (SIA) and of a multiple-SIA (n-SIA) cluster. A dumbbell or dimer is another type of interstitials, which has two atoms resting on one lattice site in perpendicular direction to a crystallographic direction. A crowdion is the addition of an interstitial atom along a crystallographic direction, such that the atoms arrange more compactly. The mechanisms of these two interstitials involve an atom displacing another atom for some space to stabilize both atoms. Figure 2 shows the con¯guration of a self-interstitial, a dumbbell and a crowdion.
Regardless of the forms of interstitials, the major concern about interstitials is mainly the interstitial formation energy required for an extra atom to become stable at its occupied space. Another major concern is the mechanism of interstitial movement, i.e., the directions involved in cluster migration. The di®usion coe±cient of the interstitials, usually expressed by using the Arrhenius plot, are sometimes investigated.
Many studies focus on the above-mentioned issues and employ MD as the tool of investigation. Here we study some important¯ndings of interstitial simulations, together with the potentials used in each study. For example, Osetsky and coworkers 96 have attempted the study of Fe clusters evolution by means of MD. By using two potentials suitable for modeling Fe and Cu, 61,62 stable 1=2h111i interstitial loops and glissile h100i loops would be formed on BCC iron. It is found that all SIA clusters would turn to be glissile, even for those clusters that are initially sessile. This implies that the accumulation density of defect in BCC iron could be much lower than that in Cu. As another example, Marian et al. 97,98 have used MD to study the SIA migration in -Fe and in Fe-Cu alloys by using a potential for Fe-Cu systems developed by Ackland and coworkers. 61 It is found that introduction of Cu impurities of 1.0 at.% in Fe results in a lower prefactor of the Arrhenius behavior of SIA di®usion. Also, the migration energy for small clusters that involve three-dimensional motion is found to be larger than that for large clusters that involve linear motion alone. The oversized substitutional Cu solute causes a dilational strain in Fe lattices, which leads to the drop of e®ective migration energy. These¯ndings are attributed to the atomic displacement¯eld interaction which triggers the change in di®usional behavior of atoms with atomic con¯gurations. Based on the variation of the migration energy with the interstitial cluster size, power-law expressions of the prefactor and the migration energy of Arrhenius plots have been extrapolated for a larger cluster size. A more detailed exposition of the mechanisms of self-interstitials of -Fe has been given by Wirth and coworkers, 99 who have calculated the stability conguration of -iron by using the FS potential 40 modi¯ed by Calder and Bacon. 100 They have found that the most stable self-interstitial of -iron is the h110i dumbbell, followed by the h111i dumbbell and h111i crowdion. They have also observed that the migration mechanism of self-interstitials in -iron is composed of two steps. The¯rst step is rotation from the ground-state h110i dumbbell to h111i dumbbell, and the second step is translation along h111i direction through the crowdion saddle point. However, it is noted that the e®ect of angular dependence of bonds in -iron cannot be demonstrated in this work, in view of the adopted FS potential which has not considered this e®ect. Tapasa et al. 101 have investigated the carbon interstitials in -iron by means of two potentials for Fe-Fe interactions 50 and Fe-C interactions, 102 respectively, for use in MD. The study shows that all clusters of SIA of size larger than seven would transform into 1=2h111i con¯guration and migrate along their crowdion axis direction, consistent with the¯ndings in Ref. 96. Furthermore, the study notices that introduction of carbon impurities inhibits cluster mobility. Two C atoms can delay transformation of h111i dislocation loops, and can even make a cluster sessile. Terentyev and coworkers 103 have studied the three-dimensional cluster motion of BCC iron by MD, with the help of the AMS potential that has incorporated DFT computations of SIA of -Fe. 50 Di®usion coe±cients and jump frequencies derived from this potential have been checked against those derived from other empirical potentials too. The study veries the h110i dumbbell as the most stable self-interstitial con¯guration, followed by the h111i crowdion. Their respective formation energies obtained in this study are closer to the DFT results than those by other empirical potentials. The study also suggests that the jump mechanism determined by Johnson 104 is the one that agrees with the DFT results. This observation is di®erent from the established results in Ref. 99. Three major mechanisms have been obtained for -iron: single and diinterstitials in fully three-dimensional motion, 3 to 5 SIA clusters in mixed one-dimensional and threedimensional motion, and SIA clusters of larger size in preferentially one-dimensional motion along h111i directions. Figure 3 compares the major¯ndings of the jump mechanism of -iron clusters.
Studies of Fe-C alloy continue to progress. A more recent study of interstitials of Fe-C alloys can be found in Ref. 105, which involves more recent potentials that considers covalent bonding of carbon explicitly. 77

Vacancies
A vacancy is a point defect which occurs when the lattice sites in a crystal structure outnumbers the atoms, such that the atoms have some space to switch their site. Figure 4 shows the vacancy created by the lack of an atom at a lattice site. Usually, the In some cases, however, the atoms in the vicinity of the vacancy would come close together, making a large vacancy. A vacancy formed by a displaced atom becomes an interstitial atom. Therefore an atom-vacancy pair is formed known as a Frenkel pair. As a displaced atom recombines with a vacancy, the Frenkel pair disappears. The major concern of vacancies is similar to that of interstitials and dumbbells. We care about the vacancy formation energy, vacancy migration energy and the jump and di®usion mechanism of vacancies. Here we introduce some representative papers that focus on Fe. The vacancy binding energy and time evolution of copper precipitates containing vacancies that are found in -Fe have been studied by large-scale MD simulations. 63,106 From this illuminative study, the vacancy binding and migration energies of -Fe have been calculated. The interaction between precipitate and vacancy is found to be anisotropic, preferentially along h011i and h111i directions of the precipitates. The anisotropic interaction suggests the tendency of precipitate phase transformation. One can realize that the di®usion behavior of vacancies within Cu precipitates depends on the vacancy concentration. Besides, the study identi¯es the three stages of time evolution of vacancies within Cu precipitates. Thē rst stage refers to the free migration of vacancies. The second stage is the clustering of vacancies that are free to move initially. The third stage is the di®usion of vacancy clusters. It is also found that the di®usion of vacancy clusters (the third stage) has a larger correlation factor than that of free migration of vacancies, which is of random-walk nature. On the other hand, monovacancy migration within a precipitate results in a smaller correlation factor than a random walk in bulk Fe. The growth of larger Cu precipitates on -Fe has been studied as well. At high vacancy concentration, the time evolution of the precipitates results in partial transformation of their atomic planes from BCC to FCC. The notion of three stages of cluster formation with vacancies has later been challenged by a similar study of Cu di®usion in -Fe, 107 using an Ackland potential for Fe and Cu, 61 from which only stage 1 and 3 can be identi¯ed. Arokiam et al. 108 have performed Cu di®usion in Fe by vacancy mechanism, showing that the di®usion coe±cients are similar for Fe atom in pure Fe in Fe-Cu alloy. The similarity is attributed to the weak interaction between the vacancy and copper atom, and to the short-range behavior of vacancy-Cu binding energy. The study also indicates that single 1=2h111i vacancy jumps dominate the simulation. On the other hand, as the temperature increases to 1500 K, vacancy double jump in h111i direction occurs. Irradiation damage study of nuclear reactors rely heavily on the interatomic potentials of -Fe. Helium is an element that can be generated in fusion reactors, so the Fe-He compound is a major topic for materials scientists. The e®ect of helium clusters on -Fe has been investigated by large-scale MD simulations by changing the number of He atoms and number of vacancies in a cluster. 109 It is found that the binding energies involved in the clusters and the Fe matrix depend on the He-vacancy ratio of clusters, i.e., the number of He atoms to that of vacancies. The binding energies are not dependent on the cluster size. The thermal stability of clusters is also dependent on the He-vacancy ratio, which controls the thermal emission of defects. Another extensive MD simulation of He clusters on -Fe has been performed, 110,111 in which the mechanisms involved in He-vacancy formation and recombination are investigated by using several Fe-He potentials and Fe matrix potentials. The study indicates that the Fe-He potential is a more important factor than the Fe matrix potential to determine the di®usion coe±cient of single He atoms. The additional binding energy required for a He atom to join an interstitial He cluster or a Hevacancy cluster has also been determined. The results show that the speed of Frenkel pair formation and He clustering in -Fe vary quite largely with the potential. The study also shows that He bubbles would expand its radius as more He atoms join the bubbles. The dilation of He bubbles developed after He-vacancy clustering is dependent on the He/V ratio, i.e., the ratio of the number of He atoms to that of vacancies within a bubble.
The vacancy mechanisms of carbon in iron have also been studied. By MD simulation implemented with the AMS potential, 50 Tapasa et al. 112 have determined the jump mechanism of a C atom toward a vacancy and of a vacancy toward a C atom. The corresponding activation energies obtained from MD are similar to those from molecular statics (MS) calculations. Another study of irradiation defects in Fe 113 indicates that vacancy di®usion in Fe-C alloys is faster than soluble carbon di®usion. This means that carbon has an e®ect of retarding microstructure evolution.

Displacement cascade
A neutron colliding with solid or liquid metal can result in displacement cascade of nuclear reactors, in which the atoms with energy greater than a threshold may experience a permanent displacement. The¯rst displaced atom after irradiation is known as the primary knock-on atom (PKA). If the energy of a PKA is su±cient, further displacement of other atoms can occur. The displaced atoms will form point defects and clusters as they migrate to other parts of a reactor. Then the resulting point defects and clusters may continue to migrate and interact, causing an even large cluster and changing the microstructure of the reactors more severely.
The study of displacement cascades typically involves realizing the formation of cascades after bombarding with a PKA with a certain recoil energy, and measuring the production e±ciency of subsequent defects. Calder and Bacon 100 have performed the¯rst MD simulation of displacement cascades of -Fe by using the FS potential, which is modi¯ed to cater for the pressure-volume relation of real metals. It is found from this study that the morphology of cascades would change to collisional phase when the PKA energy is about 1-2 keV. Generally, the number of vacancies and interstitials are greatest in the collisional phase. After a longer time, recombination is prevalent. The relaxation time for vacancy-interstitial recombination is shorter than that of the thermal spike phase. The materials after collisional phase has mainly become a hot solid, instead of a liquid state. Vacancy clustering is not found to occur in them thermal spike phase. Another early attempt has been made by Stoller 114 on displacement cascade in -Fe. The number of surviving defects is found to follow a power law with the cascade energy. The MD simulation can also suggest the presence of threedimensional clusters, which opposes the idea that only planar clusters can be formed. Even a longer simulation time of 100 ps cannot return the threedimensional morphology to the planar one. After that, Soneda and Diaz de la Rubia 115 have performed a large-scale MD simulation on displacement cascade of -iron at 600 K using the Johnson-Oh EAM potential, 23 with the recoil energy of the PKA ranging from 100 eV to 20 keV. They have successfully demonstrated the relations between the recoil energy, cluster size and number of clusters. A larger number of clusters can be formed if the recoil energy is beyond 5 keV. A cluster size of over 10 can be formed as the recoil energy is beyond 10 keV. Figure 5 shows the distribution of these relations.
The MD study of displacement cascade also suggests that a cascade is likely to split into smaller cascades of lower recoil energy. These smaller cascades are also possible to combine back to larger clusters. Similar work has been done on -Fe with carbon in solution form 116 by employing an adapted potential used to account for the short range interactions occurring in Fe-Fe solution form. It shows that for 600 K, the carbon concentration in Fe solution generally has indiscernible e®ect on the number of vacancies formed per cascade. Overlap of cascades has also been studied 117 in which a cascade produced by a larger recoil energy can mask the defects due to that by a smaller recoil energy. Displacement cascade of -Fe by a combination of MD and binary collision approximation (BCA) has been attempted, 118 which is found to be complementary to conventional MD in modeling primary damage. Application of BCA to MD provides similar qualitative results compared to full MD approach, as re°ected in the pair correlation of vacancy-vacancy separation. For a formalism of BCA, which is commonly used for computing atomic trajectories due to elastic and inelastic collisions in a lattice by considering interatomic potentials as well, readers may refer to Refs. 119-121. The BCA is applicable to modeling displacement cascades because it brings about a fairly low variance of statistics. 118 The number of cascades to simulate can be further reduced to save computation time.
Apart from using conventional MD approach, heat propagation in continuum performed by a thermal block surrounding the link cells has been proposed. 122 In this way, one can model the heat dissipated to the surrounding materials that can a®ect the mobility of SIA formed. By this approach, the number of point defects during cascade decrease with the irradiation temperature. On the other hand, the fraction of interstitials in clusters increases with the irradiation temperature. The e®ect of the PKA mass on displacement cascade of -Fe has been studied 123 using the AMS potential. 50 C, Fe or Bi atoms of varying recoil energy are made to knock on -Fe, from which we see that a heavier PKA produces fewer point defects, and that such an e®ect is more pronounced for lower PKA energy. It is the PKA mass instead of the PKA-Fe potential that is crucial for the damage in individual cascades.
A number of studies have been undertaken to examine the choice of an interatomic potential which better models displacement cascades of -Fe. Becquart et al. 124 have performed displacement cascade simulations with some EAM potentials, concluding that the short-range interaction is crucial for studying irradiation damage. Suitable repulsion mechanisms have to be formulated for better description of cascade morphology. Equilibrium properties of the potentials, such as vacancy migration energy and vacancy binding energy, are also important for modeling cascades. A comparative study of the potentials used for displacement cascades can be found in Ref. 125, which shows that MD simulations using DD, 46 AMS 50 and MEA 45 potentials can produce comparable results of Frenkel pair production. Their major di®erence lies in the defect clustered fraction. Malerba 126 has performed an in-depth review of cascades in -Fe simulated by a number of interatomic potentials and their adapted forms. This review¯nds that defect production energies calculated by various potentials exhibit essentially consistent values. It also tells us that, for the potentials attempted, the minimum displacement energy has little e®ect on the¯nal number of Frenkel pairs developed by a given recoil energy. The approach to model the potential at mid-range interatomic distance is found to be more decisive of the cascade behavior. For other critical reviews of the e®ect of interatomic potentials on displacement cascade, readers may refer to Ref. 127 for pure -Fe and Ref. 128 for He-vacancy clusters within -Fe.

Dislocations
Three types of dislocations exist: edge dislocation, screw dislocation and partial dislocation which is the combination of the previous two. An edge dislocation has a Burgers vector perpendicular to the dislocation line, while a screw dislocation has a Burgers vector parallel to the dislocation line. A partial dislocation is an intermediate case where the Burgers vector and the dislocation line intersect at an oblique angle. In many cases, a pure dislocation plane does not exist. Instead, the dislocation plane has to form kinks in order to become stable in a crystal structure. The kink structure aims to stay at the minimum points of the Peierls potential, created by the bulk atoms, that tends to prevent the dislocation plane from gliding. Figure 6 shows the edge dislocation, screw dislocation and a double-kink (DK) dislocation plane that tries to reside on minimum Peierls potential.
By applying MD, trajectories of dislocation glide can be traced out. For example, MD simulation of a=2h111i DK screw dislocation in iron has been performed for evaluating the critical resolved shear stress (CRSS) required to slip a dislocation plane. 129 The resulting glide is found to occur along a (110) plane. The MD critical stress is shown to exhibit temperature dependence as well. The migration of DK dislocation under stress has been simulated by MD approach, 130 showing that the dislocation travels from one Peierls valley to another by the application of moderate stress. As the stress further increases, the travel of dislocation becomes rough and the dislocation begins to rupture into interstitials and vacancies. Pinning points of the dislocation can even be developed, further hindering dislocation migration.
Many of the MD studies of dislocations consider the hardening e®ect of impurities in the course of dislocation glide and climb. For example, the hardening e®ect of Cu precipitates on edge dislocation in BCC iron crystals has been investigated by MD simulation. 131 The corresponding interatomic potential is based on EAM, which superposes all the Fe-Fe, Cu-Cu and Fe-Cu interactions. The MD result con¯rms the notion that the hardening e®ect increases with the diameter of precipitate atoms. In other words, a larger energy is required for an edge dislocation to penetrate through a Cu precipitate of larger radius. This hardening e®ect is caused by the introduction of a dislocation, which leads to phase instability in the particles. Similar studies have been undertaken by another group. 132 The stress-induced interaction between an edge dislocation and voids or Cu precipitates has been studied at¯nite temperature. In the course of dislocation-void interaction, the critical stress decreases with temperature. The dislocation velocity also determines the possibility of passing through the void. In passing through Cu precipitates, the critical stress also decreases with temperature. The dislocation line shape is also found to be dependent on the critical stress. Besides, the approach of passing through Cu precipitates is size dependent. Simple shear displacement occurs on small precipitates, while dislocation climb occurs on large precipitates, accompanied by phase transition of Cu from BCC to FCC. A large-scale MD study, using the AMS potential, 50 has been undertaken 133 to investigate the e®ect of void on hardening of Fe. The MD results are compared to those calculated with one Fe-Cu potential. 61 The simulation result corroborates the fact that voids are strong obstacles of edge dislocation motion because they would deform the dislocation to the screw one at low temperature. Other than this, the dislocation behaviors from the two potentials are basically the same for temperature beyond 100 K. Interested readers may also refer to other studies aimed at investigating the hardening e®ect of Fe due to Cu precipitates. [134][135][136] The e®ect of other impurities around dislocation in iron has been investigated. For example, an MD study of dislocations involves understanding the trajectory of impurity atoms in the vicinity of dislocation cores. For example, the trajectories of He atoms being put close to an edge dislocation core in bulk -Fe have been simulated. 137 It is found that, at 100 K, He atoms on the tension side would migrate to the layer closest to the slip plane as a crowdion atom. The atomic motion is driven by the interaction between He and Fe atoms. On the other hand, a He atom initially on the compression side would travel for a much shorter distance in parallel to the dislocation core, and would become stable at an octahedral site. The much restrained motion on the compression side is due to the higher activation energy required to leave the slip plane.   shows the simulated trajectories. The di®usion of hydrogen atoms in the vicinity of dislocations on a f112g slip plane has also been studied by MD simulation. 138 It is found that H atoms are strongly trapped in the vicinity of the edge dislocation core, so that the corresponding atomic di®usion is very limited. The di®usion of H atom is more mobile as its initial location is 1 or 2 nm beyond the dislocation core.

Grain boundaries
A grain boundary (GB) is the interface between two grains of di®erent orientations. Such a mis-orientation can be quanti¯ed by the mis-orientation axis and angle. One type of GB is the tilt boundary, in which the mis-orientation axis (also called tilt axis) lies in the boundary plane. Another type is the twist boundary, in which the mis-orientation axis (also called twist axis) is perpendicular to the boundary plane. Many of the real situations exist as a combination of both types of GB. Figure 8 shows a simple situation of both tilt and twist boundaries.
An early study of vacancy di®usion in the tilt boundary of -Fe at various temperatures has been undertaken. 139,140 Vacancy jumps have been iden-ti¯ed, with the probability of multiple jumps increasing with temperature. The vacancy jumps are found to be more frequent than in the bulk con¯guration, with a higher tendency of the jump direction along the tilt axis. Short-lived Frenkel pairs would also be established at elevated temperatures. The atomic vibration near the GB has a higher frequency than the bulk counterpart. By tracing the trajectories, it is noticed that vacancies jumps across adjacent layers are found to be more preferential than within the same layer. Also, at 1300 K, the atomic vibration in the GB region is more vigorous than in the bulk region. The GB structure at this temperature is found to be stable, such that the vacancies near the GB can be readily identi¯ed. The simulation can also re°ect the increase in vacancy jumps as the temperature goes to a higher value of 1500 K. Fewer vacancy jumps occur at places far away from the tilt boundary.
Impurity di®usion in GB of -Fe is also a topic of interest. For example, motion of He interstitials in GB of iron has been investigated by Gao and coworkers. 141 A series of MD simulation shows that the maximum binding energy for substitutional/interstitial He atoms and the GB are highly correlated. By¯nding the activation energy of He atoms during di®usion, He interstitials are mobile in the GB. He interstitials primarily di®use in one-dimensional path at low temperature (600 K), and in two-dimensional and three-dimensional paths at higher temperature (800 K and 1200 K). Also, a He atom in a GB would tend to di®use along the GB direction in one-dimensional path. Another MD study focuses on migration of C and H interstitials along GB in -Fe. 142 According to the MD results presented by means of Arrhenius plots, the GB decreases the mobility of H and C atoms in the vicinity of GB, because the activation migration enthalpy across GB is larger than that in bulk -Fe. In other words, GBs in -Fe can trap H and C atoms. Experimental work in this study veri¯es the trapping of C atoms in GB. The penetration of H atoms in GB is found to be more di±cult than in bulk crystal. Figure 9 shows the collection of H atoms in GB, together with the decreased penetration distance for H atoms.
Physical e®ects due to movement of GB have been considered in the literature. For instance, sliding of GB in -Fe has been studied using MD. 143 The e®ect of applying shear stress on tilt GB has been investigated by MD at 1 K, as opposed to MS technique at 0 K, which can indicate the creation of dislocation pairs acting oppositely to each other. The net e®ect of these opposing dislocation pairs results in migration of GB. Besides, the critical shear required to nucleate a dislocation decreases with temperature. In another study of symmetrical tilt GB, minimization of GB energy can be achieved by appropriate translation of grains in GB planes and by speci¯c adjustment of the tilt angle. 144 As another example, the interaction between a dislocation and a tilt GB has been simulated. 145 In this MD study, a tilt GB has been built between a dislocation and a free surface. The result shows that the dislocation glide is determined by the competition between the GB and the free surface. The attraction is strongest when the glide plane is perpendicular to the GB.
Alloy composition ratio is also crucial to the GB energy. Such a study has been performed on a symmetric tilt GB inside the Fe-Cr alloy. 146 The structure of the GB is found to remain stable in the course of thermalization, regardless of the increase in at.% of Cr. It is also noted that the GB energy decreases with higher at.% of Cr. The heterointerface established between Fe-Cr alloy and pure BCC iron has also been simulated. Such a heterointerface structure can be maintained during thermalization, and the GB transition energy is not correlated to increase in at.% of Cr.
The role of GB in fracture of -Fe has been studied by MD simulation. 147 It is noted that a GB stops crack propagation across iron. Besides, the intergranular crack propagation is determined by the angle between a GB and a crack plane. From this study, it is found that the fracture behavior of nanocrystalline materials should be linked to GB accommodation, GB triple-junction activity, grain nucleation and grain rotation. A high volume fraction of GB inside nanocrystalline materials controls crack propagation.
The e®ect of displacement cascades in the vicinity of GB has been investigated extensively. [148][149][150] A PKA of 1 keV, initiated from various directions hits -Fe bulks at varying distances from the tilt and twist GB. The subsequent time evolution of the GB has been recorded. The study shows that the tilt angle has little correlation with the GB energy, such that the GB energy can be regarded as stable. The GB has shown its role as a partial barrier to collision cascades. The PKAs cannot penetrate the GB, and defects (mainly dumbbells) accumulate near the GB in turn. In fact, during collision cascade, the GB su®ers more damage than a bulk lattice does, which can be re°ected by the number of defects formed in these two cases. Most interstitials are formed at the core portion of a GB. In addition, some preferential sites in the (5 3 0) symmetric tilt are discovered for some interstitials after the collision cascade. The largest energy change of collision cascade near a GB occurs at the¯rst half picosecond, after which the energy remains stable in general.

Nanotubes
MD simulation of iron has played its role in understanding the growth of single-walled carbon nanotubes (SWNT), which relies on transition metals such as iron to form metal carbides which act as a catalyst. Carbon atoms are then supplied to the metal carbide leading to the growth of a nanotube. The direction of the nanotube growth is also dependent on the interaction between C atoms and the metal carbides. Such a process is known as the vapor-liquid-solid (VLS) model.
The major concern about nanotubes is its high surface-to-volume ratio. Some adjustments of the MD approach is thus necessary. The MD approach used in simulating nanotubes is basically the same as that used in understanding defects, except that more re¯ned methods are used to¯nd the interatomic and electronic force. Understanding the contribution due to the bonding between C atoms and transition metal atoms is crucial for subsequent time evolution of the nanotube. In view of this criterion, special forms of MD, such as the reactive empirical bond order (REBO-type) MD, 151 quantum mechanical (QM-type) MD and density-functional tight-binding (DFTB-type) MD are adopted to calculate the interatomic potential and electronic forces by quantum mechanics. By this advanced approach, the interaction between the hybridization of C atoms and the orbitals of transition metal atoms can be better evaluated. Once the electronic forces are evaluated, the atomic trajectories are evaluated by conventional force integration. Interested readers may refer to a critical review of the SWNT 152 for an elementary understanding of the various forms of MD approaches to simulating nanotubes.
A number of MD studies involving the nucleation and growth of carbon nanotubes with the use of iron as a catalyst can be found in the literature. For example, the thermodynamics of iron carbide clusters occurring in carbon nanotubes has been investigated by MD. 153 The carbon content dependence of the cluster melting points has been calculated by cooling from much higher temperatures than their melting point. The results show that, for the range of carbon content adopted, the melting point would decrease and then increase due to the formation of stable Fe 3 C phase. The variation of melting point with Fe cluster size has also been obtained showing that the surface e®ect in simulation would result in a lower melting point than the experimental bulk condition. Introduction of carbon also lowers the melting point of Fe clusters because carbon atoms would break the symmetry of Fe clusters thereby destabilizing the structure. The study deduces that below 1200 K, nanotubes might grow from a solid particle or from the molten surface. In another study, the time evolution of growing SWNT on iron nanoparticles have been attempted by means of ab initio MD. 154 It is found that fast carbon di®usion occurs on the metal surface with carbon dimers followed by carbon sp 2 pentagonal and hexagonal bonded caps rooting on top of Fe catalyst without carbon penetration into Fe. The results also indicate that stabilizing a C atom at the Fe cluster surface is more favorable than at the cluster core. The binding energy calculation of Fe on C shows that SWNT growth is possible in iron, where a stronger covalent bond and higher adhesion energy can be achieved. Surface melting of Fe 586 Wul® polyhedral clusters occurring on SWNT has been studied to realize its mechanism. 155 Calculations indicate that the molten surface state can occur at a temperature below the melting point of the clusters. The temperature dependence of Lindemann index (LI) 156 (a measure of atomic disorder) tells us that the melting process of an Fe 586 cluster can be split into three stages. Thē rst stage refer to the slow increase in LI with temperature, while the second stage is the abrupt and nonlinear increase in LI, corresponding to surface melting. This means that the atomic kinetic energy can overcome the binding energy at the cluster surface. The¯nal stage refers to complete melting at high temperature, returning the maximum LI. The graph indicating the temperature dependence is shown in Fig. 10. The radial distribution of LI indicates that surface melting is more prominent at higher temperature. The growth of SWNT on iron clusters has been simulated by MD formalism 157 showing that the growth favors the perpendicular direction that has weaker interaction between the SWNT and the supported substrate of aluminum oxide. The growth angle also increases with simulation time due to the carbon-substrate interaction, which favors the presence of precipitated carbon atoms. In order to reduce the formation along the perpendicular direction, the SWNT-substrate interaction should be increased. MD simulation of coating iron onto SWNT has been studied. 158 By simulating continuous metal evaporation coating, one can see that iron clusters can combine with carbon¯rmly and provide an outward pull on carbon atoms leading to structural deformation of nanotube. The iron atoms tend to cluster together and form second layers.
6. Hardware and Software Development for MD

Recon¯gurable computing
With the advance of¯eld programmable gate arrays (FPGA) as an external device used in cooperation with the conventional CPU machines, one can congure a hardware board that is dedicated to a certain processing stage easily, so that it can accelerate computations that are undertaken ine±ciently in CPU. The approach to using hardware components such as FPGA in operation with appropriate hardware programming is known as recon¯gurable computing. As the phase space trajectory of a particle is practically irrelevant to that of another distant particle, FGPA can come into play for an accelerated performance that a CPU focusing on sequential processing cannot achieve. Accordingly, FGPA is a candidate of high-speed MD computation.
The working principle of recon¯gurable computing is di®erent from CPU computations. In fact, the design principle of recon¯gurable computing shifts to formulating the connection among various memory blocks and logic blocks as well as to devising the force pipeline that results in the interatomic forces used for motion integration. 159 Therefore, developers of applications in FPGA need to program the hardware connections every time a new algorithm is adopted. Despite the hardware-level approach, the idea behind the work°ow remains the same as those applied in CPU and GPU. For example, computational scientists have to organize the data°ow between the host (CPU) and device (FPGA in this case).
Currently in the FPGA market there are two popular brands. Xilinx is the leader of the market, providing the Spartan series for basic computation capability, the Artix and Kintex series for more demanding tasks, and the Virtex series for advanced tasks. Altera, as the major competitor of Xilinx, provides a similar range of products. The Cyclone series targets the Spartan series, whereas the Arria series targets the Artix and Kintex series. The Stratix series is the high-end series o®ering the same level of computational power as the Virtex series.
While recon¯gurable computing which has evolved for over a decade has found its value in MD, its application is still fairly limited to modeling biological and chemical molecules. Modelers in biochemical discipline are concerned about two major types of forcesbonded and nonbonded. The bonded force has a lower computational complexity of O(NÞ which is a®ordable by CPUs, while the calculation of nonbonded forces has a higher complexity of O(N 2 ) and is hence suitable for hardware acceleration. Lennard-Jones (LJ) force 160 is the short-range interaction that is deterministic of the resulting interatomic force. It is derived from the potential which has the form where r ij is the interatomic separation, " is an energy parameter and is a distance parameter. The short-range force is then obtained by numerically di®erentiating this potential with respect to interatomic separation. The velocity and position of individual particles can be calculated by motion integration techniques. Some of the worldwide computation systems that are adaptable to FPGA boards, such as GROMACS, 161 MD-GRAPE, 162 MDGRAPE-2, 163 NAMD 164 and MODEL, 165 have succeeded in processing the computationally complex LJ force and Coulombic force. In order to accelerate the force computation, lookup table storing the LJ potential as a function of r ij might be used instead of direct computation of LJ potential by putting interatomic distances into Eq. (71). Cuto® distance is often used in short-range computations to increase the processing speed, such that the force or interaction is neglected for the interatomic distance beyond a cuto® distance r C . The Coulombic force is a long-range nonbonded force that is often incorporated in the atomic interaction, which is expressed as 166 where in Eq. (72) q i is the charge of particle i, and r ij is the atomic separation vector between particles i and j. Unlike the LJ force, Coulombic force is slowly decaying. It can still have a¯nite amount even after a rather large atomic separation. The approach of using cuto® is thus inapplicable to determine the Coulombic force. 167 A notable method to solve this problem is the Particle mesh Ewald (PME), 168 which relies on computing the force in the reciprocal domain using the three-dimensional fast Fourier transform (FFT), yet the implementation in FPGA is less ine±cient. 169 The multi-grid approach is also a possible approach, whose sequential processing can attain roughly the same computational speed as the PME. 170 A number of attempts have been made to implement the computation of the Coulombic force, see for example Ref. 171.
As large-scale high performance computing is required for many MD tasks, the advantage of FGPA to allow for scaling and parallelism can be utilized. An advantage of FGPA over the conventional computing clusters is the customization of hardware components, which leads to the decreased cost and electric power. 172 The general process of MD calculation on an FGPA can be summarized in few steps, 173 similar to the steps implemented in CPU. First, the cell list is loaded to the FGPA memory, whereas the particle position data are stored in a memory location external to the FGPA. Second, the FGPA adopts the cell list and position data to generate the pairs to be used in interatomic force computation. Newton's third law is often applied, so that the number of force pairs would be reduced by half. After the force pairs are determined, the LJ force can be computed, followed by the update of the acceleration of the particles which are later stored in external memory. Figure 11 illustrates this idea by means of a schematic diagram representing a typical FPGA board.
FPGA designers have to formulate the pipelines dedicated to processing the MD tasks. Readers may refer to some of the implementations in Refs. 166, 169 and 174-176, and this review discusses two of them. Kasap and Benkrid 166 have decomposed the whole MD process into four independent pipelines, and each pipeline handles the nonbonded potentials, the resulting forces and virials due to all other particles in the simulation system. Figure 12 shows the schematic diagram of one pipeline, and the design requires four of them linked to establish the FPGA implementation. Each of these four processors has a dedicated SDRAM allocation for holding the input data. Then the input bu®er of each processor receives the data from the input SDRAM portion and transfers the data to the processor for calculation with the help of the function coe±cients used to interpolate the potentials. The calculation results go to the output bu®er after passing which the results go to the SDRAM portion of the FPGA responsible for storing data. The processors rely on the¯nite state machines (FSM) to coordinate the data transfer.
In view of the advancement of computer networking, Scrofano et al. 175 have developed another pipeline for MD processes performed on a cluster of FGPAs computation nodes. The hardware design for each FGPA node is similar to that in Fig. 11, but the nodes this time are connected to each other to establish a cluster. The schematic diagram of this idea is illustrated in Fig. 13, with parallelization of nonbonded force evaluation relying on spatial decomposition technique. 177 A number of generalpurpose processors (GPP) group themselves to form a GPP element, whereas a number of recon¯gurable hardware (RH) devices group themselves to form an RH element. These two elements are linked together for mutual data transfer. A simulation box is partitioned in a number of simulation cells containing a number of atoms, each of which is assigned an FGPA node. Each node handles the Fig. 11. Block diagram of an FPGA board whose con¯guration targets MD simulations. Reprint of Fig. 1   computations independently, except that cross-cell communication occurs during the atomic movement across the edges of the simulation box and across the cells. They found that a cluster of N accelerated FGPA nodes perform faster than a cluster of 2N computing nodes without applying FGPA acceleration, so that investment in hardware infrastructure can be reduced while maintaining high computation performance.
While recon¯gurable computing packages for MD have been ripe in biochemical discipline, the popularity of recon¯gurable computing in modeling physical phenomena is far lower. Collaborative studies between computer scientists and physicists are therefore anticipated to further extend the implementation of recon¯gurable computing to simulate metals. This can be achievable by formulating a general approach to force computation and integration algorithms tailored for metals. It is expected that metal simulations can gain advantage by using the FPGA accelerator together with appropriate MD acceleration algorithms. The pipelines reviewed in this paper might be treated as some possible guidelines of such a design for MD in metals. Although FGPA has to sacri¯ce the°oating-point precision for a larger number of function units, 176 some approaches implemented in¯xed-point calculations might still be helpful to explore the use of FGPA to materials sciences. After all, the requirement of double precision in MD simulations is questionable. 178

Computing based on GPU
In spite of the progress made by recon¯gurable computing in MD simulation, its development is hindered by the complexity to design custom¯rmware and hardware dedicated to parallel computation tasks. 179 Expertise in electronic hardware and its related programming skills are therefore required for a research team to employ recon¯gurable computing. Furthermore, the programming language for electronic design is not suitable for coding scienti¯c tasks 174 adding di±culties for applying it to MD computations.
With the advance of GPU, large-scale simulation tasks can be performed with more readily accessible hardware components. The speedup factor of using GPU compared to using CPU can often reach 100. Besides, the skills requirement is lowered from understanding hardware architecture of FGPA to realizing parallel programming techniques. In the old days, GPUs were not easy to use for materials modelers because the MD codes had to transform to graphics operations manually by means of appropriate mapping. 180 Fortunately, two streams of programming architecture are now in the market that help users to perform this task. In June 2007, NVIDIA developed the Compute Uni¯ed Device Architecture (CUDA), a proprietary application programming interface (API), which facilitates multicore and parallel computing that is coordinated between CPU and GPU. In 2009, Open Computing Language (OpenCL) framework was established as an open-source and cross-platform counterpart of CUDA, so that parallel processing can also be performed on Intel CPUs, AMD CPUs and AMD GPUs. Nowadays OpenCL is already incorporated in the AMD APP SDK as a tool of computing using ATI GPU cards. Multicore computation units controlled by hardware programming language are applied to complex calculations, with proper code optimization for better parallel computing performance. Supercomputing can hence be performed in software level, without touching on the implementation of FPGA. Moreover, CUDA and AMD APP SDK can be executed on GPU cards initially focused on video gaming, in the sense that more users can experience the high-throughput computations with ease at a more a®ordable cost.
Another advantage of GPU over recon¯gurable computing is the abundance of open-source code libraries for a fairly accessible GPU architecture, as compared to the lack of such programming resources for FPGA. 179 Both GPU card manufacturers release their own code libraries, and some third-party libraries are available for download free of charge. Apart from a series of NVIDIA libraries such as cuBLAS and cuFFT, notable third-party libraries include CULA which acts as an alternative linear Fig. 13. Idea of connecting various FPGA nodes to form a cluster, which is found by the authors of Ref. 175 that N nodes with acceleration perform better than 2N nodes without acceleration. Reprint of Fig. 1  algebra library, JCuda which combines CUDA operations to JAVA libraries, and PyCUDA which allows usage of CUDA codes in Python environment. For AMD, the AMD Core Math Library (ACML) supports the use of AMD APP SDK.
A GPU performs parallel computations by the single-input-multiple-data (SIMD) architecture. The same instructions would be performed on each thread of execution which processes di®erent data values. Since the time to execute each thread can be slightly di®erent, before performing another set of computations, processed threads have to wait until all other threads have¯nished their computations; threads of execution have to be synchronized. Data transfer between the GPU and CPU interweave the computation process, so that data are uploaded to the threads and the results are downloaded to the host machine. MD simulation follows the SIMD architecture, therefore it is suitable for GPU processing.
The general mechanism of CUDA in transferring data and instructions (known as kernels) between the host and the GPU device is depicted in Fig. 14. The host memory containing variables to be evaluated is copied to the device memory, at which the data required are transferred to the CUDA cores for computation to take place. The computed results are stored in the device memory, which are then copied back to the host memory. For computationally intensive processes, one may further raise the speedup by using the shared memory, a fast memory of the order of KB per multiprocessor. 181 Because of its limited capacity, delicate and skillful organization of data transfer to the shared memory is necessary to ensure that only the toughest part of computation is performed by it. At the end of the computation, developers have to free the memory in the GPU, which is initially used to store variables.
Legacy CUDA versions require memory transfer to be performed explicitly, which discourages developers accustomed to pure software programming from employing this API. This shortcoming has been overcome since the release of CUDA 6.0, which allows developers to use the uni¯ed memory that is shared between the CPU and GPU. 182 Explicit transfer actions are no longer necessary, such that CUDA programming is more amenable to software developers whose focus is not on GPU hardware architecture.
AMD APP SDK, on the other hand, operates on another scheme by supporting OpenCL as the language of parallel computations performed on AMD graphics cards. Unlike CUDA which links both CPU and GPU chips, AMD APP SDK is merely the runtime for CPU, so users have to install the Catalyst driver for AMD cards, which includes the runtime for the GPU component. OpenCL employs roughly the same terminology as CUDA to construct the programming architecture. 183 For example, parallel algorithms are performed by kernels. The smallest working unit of OpenCL is known as a work item, which is conceptually equivalent to a thread of execution in CUDA architecture. Memory allocation has to be done on the GPU, and the results generated by the GPU have to be copied back to the host. Also, users have to free the GPU memory after use.
Di®erent collections of NVIDIA GPU cards can help to perform computations of varying complexity. For supercomputing level, the latest Tesla GPU to date has been installed in many famous computing clusters such as Titan, the top cluster till November 2012, 184 which is used for materials science, nuclear engineering and climate research. The relatively cost-friendly Titan series of NVIDIA GPUs can reach the Tesla computation capability of over about 1.3 TFLOPS for double-precision tasks, 185 at the expense of using non-ECC RAM, i.e., memory without error correction code functionality. For single-precision computation, the GTX series primarily designed for video gaming market can already provide a modest speedup as compared to single-core CPU computation.
AMD, as the competitor of NVIDIA, has its series providing comparable computational power. The R series (formerly known as the Radeon series), typically used for video gaming, targets the NVIDIA GTX series. The FirePro series is the collection for professional computation tasks, which can provide over¯ve TFLOPS (W9100 model) for single-precision calculations.
With the rapid growth of the number of transistors assembled in a GPU unit, which is faster than that of a CPU unit, the bottleneck of GPU computation speed lies instead in the bandwidth provided by the PCI Express slots used to transfer data between the host and the device. 186 Worldwide researchers of MD simulations have suggested some algorithms to utilize the GPU architecture and its parallel computation capability that is not the strength of CPU (see  for some of these examples). Besides, open-source GPU software packages such as OpenMM 193 and its later version 194 are also available for public downloading. Some of the aforementioned examples are discussed in this review paper.
Using CUDA, Anderson et al. 188 have performed a comprehensive study of MD simulation on a single GPU not by adapting the CPU code but by completely rewriting a set of code optimized for GPU cards. Particle data are stored in the global memory of the device accessible by all CUDA cores, from which the data is loaded to the texture memory where summing of forces is undertaken. In order to utilize the property of GPU in matrix computation, the neighbor list used to calculate the short-range force is organized in a matrix form instead of the conventional linked list format. Each thread sums the pairwise forces due to the atoms in the neighbor list concurrently. However, the authors did not apply the Newton's third law to¯nd the interatomic force, in fear that it would incur additional memory latency during the read-modify-write process. To reduce the time required to transfer data between the CPU and GPU, integration of the equations of motion is performed in the GPU instead of the CPU, soon after the interatomic forces are determined in the GPU. In simulating a large-scale system of over 1 million particles, the authors demonstrated a speedup factor of about 60 in¯nding the LJ force, and about 30 in generating the neighbor list.
Here we brie°y review another implementation of GPU computing by OpenCL, executable on AMD cards as well as NVIDIA cards. Michael Brown et al. 192 have established an implementation scheme for LAMMPS, 177 a famous object-oriented MD simulation package. In essence, the authors added the OpenCL code simply by adding a derived class of the original code used without GPU acceleration.
The authors were basically doing parallel decomposition of the MD processes including¯nding the neighbor list, calculating the LJ force and the Gay-Berne (GB) force, 195 and integrating the equations of motion. The idea of acceleration mainly lies in the operation of the neighbor lists in the GPU rather than in the CPU. The neighbor list is an improved version, which takes advantage of the linked cell list, and hence reduces the number of particles to check in each time step. The authors also balanced the load by overlapping the short-range computations in GPU with the long-range counterpart in CPU, such that calculations are undertaken on both host and device concurrently. The LJ force requiring low arithmetic intensity, whereas the GB force, a mod-i¯ed LJ force requiring high arithmetic intensity, was employed to demonstrate the system speedup of OpenCL implementation. It was found that the speedup of¯nding the LJ force could reach between 2.9 and 7.8, with a longer cuto® distance returning a higher speedup. The speedup of¯nding the GB force could even reach 11.2, because more arithmetically intensive operations performed in parallel can hide the memory latency between host-and-device transfers during the complex operations involved.
In order to promote the application of OpenCL in scienti¯c computing community, some source code mapping tools have been formulated to facilitate code translation from CUDA to OpenCL. [196][197][198] Though such a translation is trivial in view of the similar architecture of both CUDA and OpenCL, challenges exist for a robust porting between them. 199 For example, separate compilation of CUDA source¯les is possible, yet it is quite di±cult to link the code translated in OpenCL format due to the required reorganization of the initialization code throughout all the source¯les. Besides, source code calling the CUDA libraries, which are not included explicitly in the source¯les, is hard to translate to OpenCL directly.
It has been found that FPGA can still be a highly competitive choice of MD acceleration given that the hardware con¯guration and the pipelines are carefully designed. 200 However, we can realize from this section of the review that GPU is currently more favorable than FPGA for large-scale computations, in view of its accessibility and the skills involved. Also, at the time of preparing this paper, the CUDA is still more popular than OpenCL in MD simulations. It is expected that, with increasingly advanced GPU computation capability, developers can further contribute to both types of programming streams in MD simulations thereby providing more choices of programming tools to the scienti¯c community.

Summary
This review paper starts with a brief discussion of statistical mechanics, which forms the basis of the viability of MD formalism. A number of practical implementations of motion integration then follow. Some common thermostats have been mentioned to maintain the ensemble temperature. SLD for ferromagnetic materials and TI have been investigated, which act as the supplement of the conventional MD approach. The interatomic potentials for iron have evolved from using FS formalism to embedded atom method followed by magnetic iron potential. With this development, a number of iron potentials in pure form and with impurities have been formulated. Examples of applying appropriate interatomic potentials for iron to simulate the time evolution of atoms have been discussed, which are mainly related to the safety of nuclear power plants. They demonstrate that the considerations are further re¯ned when new potentials are adopted, so as to re°ect the increasingly complicated defect conditions. Recongurable computing and GPUs are common hardware components for MD simulation, yet the former is less employed in simulation of metallic materials. OpenCL is more developed than CUDA in terms of MD simulation in view of their current trend of application in scienti¯c community.