Modelling Of The Structure And Dynamics Of Nanoparticles On Metal

Transcript

Master Science de la matière École Normale Supérieure de Lyon Université Claude Bernard Lyon I Stage 2012 Daniel Förster [email protected] M2 Physique Modelling of the structure and dynamics of nanoparticles on metal substrates Abstract : Epitaxial graphene on metallic (111) surfaces has proven to stabilize metallic nanoparticles placed on top due to a moiré motif. In this work, a simplified system without the graphene layer has been examined. The modelling is based on an empirical many-body potential for the atoms of the nanoparticles and a continuous potential that treats the interaction of the nanoparticle atoms with the substrate. This potential allows to switch on and off the corrugation of the substrate surface. A molecular dynamics code featuring a Nosé-Hoover thermostat has been developed for the exploitation of these potentials. With this code, simulations were carried out on iridium nanoparticles of three different sizes: 13, 55 and 147 atoms. This lead to equilibrium structures of the nanoparticles on the substrate. Furthermore, results concerning the melting behaviour of iridium nanoparticles as a function of their size as well as the presence and corrugation of the substrate could be obtained. Key words: “Iridium nanoparticles”, “Molecular dynamics”, “Sutton-Chen potential”, “Melting” Supervising professors: Florent Calvo and Franck Rabilloud [email protected] / phone: (+33) 4 72 44 83 14 Laboratoire de Spectrométrie Ionique et Moléculaire (LASIM) LASIM-UMR 5579 CNRS / UCBL Domaine Scientifique de la Doua Université Claude Bernard Lyon 1 Bâtiment Alfred Kastler 43 Boulevard du 11 novembre 1918 69622 Villeurbanne cedex http://www-lasim.univ-lyon1.fr August 15, 2012 Acknowledgements Tout d’abord je remercie Florent Calvo qui m’a donné l’opportunité d’effectuer ce stage au sein du groupe physico-chimie théorique du LASIM. Son expertise et sa disponibilité pour répondre à toutes mes questions m’ont encouragé à me surpasser tout au long du stage. Je tiens également à remercier Franck Rabilloud pour ses conseils pertinents. Merci enfin à Ma Quianli et Bai Xueshi qui par leur bonne humeur ont rendu le travail quotidien au bureau très agréable. Contents 1 Introduction 1 2 Simulation methods 2.1 Empirical potentials . . . . . . . . . . . . . . . . . . 2.1.1 Sutton-Chen potential . . . . . . . . . . . . . 2.1.2 Implicit interaction potential with a substrate 2.2 Molecular dynamics . . . . . . . . . . . . . . . . . . 2.2.1 Basic principles . . . . . . . . . . . . . . . . . 2.2.2 Velocity Verlet integrator . . . . . . . . . . . 2.2.3 Nosé-Hoover thermostat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 2 3 4 4 6 6 3 Implementation details and observables 3.1 Adjustment of simulation parameters . . . . . . . . . . . . . . . . . 3.1.1 Interaction with the substrate and coordination dependence 3.1.2 Stability of trajectories . . . . . . . . . . . . . . . . . . . . . 3.1.3 Choosing the mass of the Nosé-Hoover thermostat . . . . . 3.2 Observables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Stable structures . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Thermodynamical properties and Lindeman index . . . . . 3.2.3 Diffusion constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 7 7 9 10 12 12 12 13 4 Results 4.1 Relaxed configuration of Irn clusters in presence 4.2 Dynamical relaxation to equilibrium geometry . 4.3 Thermodynamics of Irn clusters and Lindemann 4.4 Dynamics and diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 14 15 16 18 of a substrate . . . . . . . . . . indices . . . . . . . . . . . . . . . . . . . . . 5 Conclusion and perspectives 19 References 20 A Force calculation from Sutton-Chen potential 21 B Force calculation from continuous substrate potential 22 C Nosé-Hoover thermostat and canonical ensemble 24 D Parameters for continous substrate potential 26 Modelling of the structure and dynamics of nanoparticles on metal substrates 1 Daniel Förster Introduction Nanoparticle assembly on substrates Nanoparticles may have strongly different properties as compared to the same material in bulk form. Metallic nanoparticles such as of iron, nickel or cobalt feature large magnetic moments which makes them interesting for the use in high-density magnetic data storage devices. This technology would require not only the production of monodisperse nanoparticles, but also the capability of depositing them regularly and stably on a substrate [1]. One possibility of arranging nanoparticles is the use of metal substrates with an epitaxial layer of graphene on top. The metal substrate may consist of iridium presenting a (111) surface; a scheme of such a system is shown in Figure 1. The lattice constants of the graphene layer and the underlying bulk iridium are incommensurate. This gives rise to a moiré motif with a rhombic unit cell of a lattice repeat vector of (25.3 ± 0.5) Å in length which corresponds to about nine times the atomic nearest neighbour distance of bulk iridium [2]. Experimentally it has been observed that this moiré may act as a template for well-ordered iridium cluster patterns. In this context, also DFT calculations suggesting that the moiré corresponds to the corrugation of the graphene layer, have been carried out [2]. Iridium nanoparticles 1111111111111111111111111111111111111 0000000000000000000000000000000000000 0000000000000000000000000000000000000 1111111111111111111111111111111111111 0000000000000000000000000000000000000 1111111111111111111111111111111111111 Graphene layer 0000000000000000000000000000000000000 1111111111111111111111111111111111111 Iridium substrate Figure 1: Scheme of the system: Epitaxial graphene on iridium acts as a template for a Irn cluster array on top Present work The long-term objective of the type of research described in the following is the modelling of an assembly of nano-particles on the above cited substrate. However, in the present work the model contains only a semi-infinite block of bulk iridium presenting a (111) surface on which one nanoparticle can be placed. The modelling of the nanoparticle is based on a semi-empirical N -body potential and the interaction of the atoms within the nanoparticle with the substrate is treated by an additional continuous potential describing implicitly the periodicity and corrugation of the substrate. Further work would involve a potential that also describes the effect of an epitaxial graphene layer on the metallic surface so as to obtain corresponding results for the model system of Figure 1. The thermodynamical and dynamical properties of iridium nanoparticles of three different sizes (13, 55 and 147 atoms) are explored with the help of molecular dynamics. The effect of the presence of the substrate and its corrugation on these properties is also examined. For this purpose a molecular dynamics code featuring a Nosé-Hoover thermostat that allows to carry out simulations in the canonical ensemble has been developed and evaluated. Simulations run with this code form the basis for the latter analysis of the behaviour of iridium nanoparticles in different configurations. In section 2 a description of the empirical potentials can be found, which model the atomic interactions and the interactions with the substrate surface as well as an explanation of the working principle and the features of the molecular dynamics code. Section 3 explains the choice of the parameters for the mentioned potentials as well as the details of the simulations, such as their parameter settings and the calculation of the observables. The results that could be obtained from the simulations are presented in section 4 accompanied by a discussion. The focus is there on how the size of the nanoparticles and the presence of a substrate influences their melting and dynamical properties. Finally, this report closes with section 5 by giving an outlook on how the modelling can be extended in order to obtain more telling results in the future. 1 Modelling of the structure and dynamics of nanoparticles on metal substrates 2 Daniel Förster Simulation methods 2.1 2.1.1 Empirical potentials Sutton-Chen potential Molecular dynamics simulations depend crucially on the potential which is used for the force computation. In this work, the semi-empirical many-body Sutton-Chen potential has been employed to model the interactions between the iridium atoms of the nanoparticles. Bonding in metals is due to delocalized electrons, which would suggest to use electronic structure theory to compute the interactions. This quantum approach, however, is very time consuming making the Sutton-Chen potential the model of choice at the expense of accuracy. The potential treats the atomic interactions classically and depends on the nuclear positions only, neglecting the electronic degrees of freedom. When compared to pair potentials like the Lennard-Jones potential, it is more realistic as it takes into account the repartition of all other atoms of the system when considering the interaction of a pair of atoms. Another advantage is that it is efficient to evaluate numerically and that it allows analytic expressions for various quantities as a function of its parameters. The Sutton-Chen potential has the following form (see equation (1)). Like Lennard-Jones potentials, it has a repulsive and an attractive part, that correspond to the first and second term of equation (1) respectively. The attractive term is long-ranged and has a genuine many-body character. At small separations the potential is dominated by a repulsive term, which accounts for the Pauli repulsion:   X  a n X√ 1 VS−C = ε  ̺i  (1) −C 2 rij i i,j|i6=j X  a m (2) ̺i = rij j6=i where ̺ corresponds to the local density of atoms, which in this approximation is a good estimate for the q electron density. The separation of atoms i and j is designated by rij = (xj − xi )2 + (yj − yi )2 + (zj − zi )2 with xi , yi and zi the Cartesian coordinates of atom i. The potential contains two parameters ε, a, that set the scales of energy and length, and a parameter C that is determined by the equilibrium conditions of the particular lattice structure. The two exponents m and n characterize the potential and are restricted to integer values with n ≥ m. This means that if two metals having the same crystal structure share the same pair of exponents, all results obtained with this potential can be translated from one metal to the other by simply rescaling the units of energy and length according to ε and a [3, 4]. The set of parameters for iridium that is used for the simulations is taken from [3] and is shown in Table 1. m 6 n 14 ε [eV] 2.4489 · 10−3 a [Å] 3.84 C 334.94 Table 1: Parameters of the Sutton-Chen potential for iridium, as given in [3]. As explained in section 2.2, it is necessary to compute the forces on the atoms that derive from this potential. The derivation is given in appendix A and leads to the following expression for the force f on atom k:   m   n  1 X  r~kj a Cm a − 21 −2 + (3) fk = ε −n ̺k + ̺ j 2 rkj 2 rkj rkj j6=k In this expression r~kj = r~j − r~k designates the vector connecting k to j. 2 Modelling of the structure and dynamics of nanoparticles on metal substrates 2.1.2 Daniel Förster Implicit interaction potential with a substrate The model considered in this work contains not only the nanoparticles, but also a semi-infinite iridium substrate. However, it does not include the epitaxial layer of graphene that was present in the corresponding experiments and DFT calculations. The effect of this layer would be subject to further investigation and would require another potential similar to the one described here, when treated implicitly. As the here considered system contains only iridium atoms, the nanoparticles as well as the substrate, it could already be entirely described by the Sutton-Chen potential presented in the previous section. Such calculations, however, are quite time-consuming when taking into account a reasonable number of atoms of the substrate. So, the idea here is to treat the interactions between the atoms belonging to the nanoparticle and those of the substrate implicitly. A continuous potential that contains the periodicity and corrugation of the substrate is therefore introduced. Bulk iridium has a face-centered cubic lattice structure of which a (111) surface is considered here. This surface has a triangular lattice, which imposes the form of this potential. The potential is approximated by the three first terms of a series expansion in terms of the surface corrugation: X VSubstrate = Vi,0 + Vi,1 fi,1 + Vi,2 fi,2 (4) i The term Vi,0 represents the substrate potential in the lowest order approximation, i.e. a planar surface. It depends therefore only on the z coordinate of atom i and its coordination ci . The term "coordination" and the calculation of its numerical value is explained later on. The two other terms account for the corrugation of the surface, where fi,1 and fi,2 are the first two harmonics of the triangular lattice that depend on the x and y coordinates only and Vi,1 , Vi,2 correspond to the coefficients of this development that depend on the distance from the surface and the coordination of the considered atom. The functions fi,1 and fi,2 have the following form for a triangular lattice with a lattice constant b:         2π yi yi 2π 4π yi √ fi,1 = cos xi + √ xi − √ + cos + cos (5) b b b 3 3 3        √  √  2π 2π  4πxi xi + 3yi + cos xi − 3yi + cos (6) fi,2 = cos b b b This lattice constant b is related to the one a of the face-centered cubic lattice of bulk iridium via √ b = a/ 2. The dependence of the coefficients Vi,m on the z coordinate is similar to the form of the Sutton-Chen potential: αi,m,1 αi,m,2 Vi,m = αi,m,3 − αi,m,4 (7) zi zi The coordination of atom i is taken into account via the αi,m,l parameters. This fact ensures that the many-body character of the Sutton-Chen potential is not lost in this implicit description of the interaction with the substrate. The functional form of this dependence has been chosen empirically. This choice is illustrated in section 3.1.1: αi,m,l = βm,l [1 + dm,l exp (−γm,l ci )] (8) The numerical values for the parameters βm,l , dm,l and γm,l have been determined via a curve fitting procedure, where the parameters of this implicit potential are adjusted according to the explicit SuttonChen potential for the iridium substrate with a (111) surface. The details of this procedure are also discussed in section 3.1.1. Finally, ci stands for the coordination of atom i for which the potential is calculated. Roughly speaking, it corresponds to the number of its nearest neighbouring atoms, where it is necessary to specify the maximal distance of two atoms that are considered being neighbours. It is, however, not possible to use a sharp cut-off at a certain distance since this would lead to jumps in the potential which is undesirable for use in geometry relaxation and molecular dynamics simulations. 3 Modelling of the structure and dynamics of nanoparticles on metal substrates Instead of a sharp cut-off, a smooth transition of the following form is employed:  1 rij ≤ rmin  h  i X rij −rmin 1 1 + cos π rmax −rmin rmin ≤ rij ≤ rmax ci = csurf + 2 j6=i  0 rij ≥ rmax Daniel Förster (9) If the distance of another atom to the one in question is smaller than rmin it is counted entirely, if its distance is between rmin and rmax then it is counted partly and for any larger distances not at all. Another contribution csurf to the coordination stems from the fact that the number of neighbouring atoms is not only due to atoms of the nanoparticle that are treated explicitly by the model, but also due to atoms belonging to the substrate. These atoms do not appear explicitly in the model, and it is therefore necessary to estimate their number as a function of the distance of the atom in question to the substrate. The estimation is done on the basis of the following seventh order polynomial whose coefficients are adjusted in a way that csurf reproduces well the real physical additional coordination number due to the substrate atoms: (
zi < rmax a5 zi5 + a4 zi4 + a3 zi3 + a2 zi2 + a1 zi + a0 (zi − rmax )2 (10) csurf = 0 zi ≥ rmax If an atom is farther from the substrate surface than rmax there is no contribution of the substrate atoms to its coordination. It should be noted that the polynomial is continuously differentiable at rmax . In section 3.1.1 the curve fitting procedure for the adjustment of the parameters ai is explained. For molecular dynamics simulations demand an expression for the force on the atoms, also this potential with all its dependencies needs to be differenciated. The derivation, along with the expression for the forces can be found in appendix B. In the simulation code, there is a switch that allows to choose whether one wants to activate this potential and its corrugation. This potential is switched off for simulations of a free nanoparticle. In simulations of the nanoparticle in interaction with the substrate without corrugation, only the first term of the above mentioned series expansion is used. Finally, simulations including the the entire contribution of this potential correspond to the situation of a nanoparticle on a corrugated substrate. 2.2 2.2.1 Molecular dynamics Basic principles In order to derive physical properties from these interaction potentials, the method of molecular dynamics is employed. The idea of classical molecular dynamics is to integrate numerically the Newton equations of motion for each atom in the system, neglecting quantum nuclear effects. The forces that appear in these equations are derived from potentials like the ones described in section 2.1. This way, trajectories of the system through phase space may be generated, allowing to compute macroscopic properties as a function of the atomic coordinates. These trajectories are generally chaotic, i.e. due to the Lyapunov instability, trajectories of two systems with nearby initial conditions diverge exponentially in the course of time. This is not necessarily a problem since the aim is not to predict the precise trajectory for a given initial condition, but to sample a certain statistical ensemble. It can even be seen as an advantage, because due to this chaotic behaviour the system may explore more of the phase space, which is a demanded by the ergodic hypothesis [5, 6]. Initialisation The molecular dynamics procedure usually begins with an initialisation process, where a configuration at equilibrium at zero temperature is prescribed and random initial velocities are assigned to each atom. The velocities are to be rescaled according to the equipartition theorem so that one can already impose at the beginning an approximate temperature to the system. However, one should use a temperature about twice as high as desired, because about half of the kinetic energy will be transformed into potential energy in the equilibration part of the simulation. This is not strictly 4 Modelling of the structure and dynamics of nanoparticles on metal substrates Daniel Förster exact if the interaction potential is anharmonic, which is usually the case, but this method can put the system at least in the vicinity of the target temperature. Later on, in the course of the simulation the temperature can be determined also via the equipartition theorem, which assigns to a given total kinetic energy Ekin a temperature value T . Ekin =
f 1X 2 2 2 = kB T + viy + viz mi vix 2 2 (11) i The mass of atom i is designated with mi and is throughout this work always equal to the atomic mass of iridium. The viα are the components of the velocity vector of atom i, kB is the Boltzmann constant, and f the number of degrees of freedom of the entire system. It should be noted that f depends on the presence of the substrate. For a system without interaction to the substrate there are 3N − 6 degrees of freedom, 3N for the atoms minus six that are lost to the conservation of linear and angular momentum. In the case of a system in presence of a surface without corrugation we would obtain 3N − 3, where the minus three comes from the rotational symmetry with respect to the z axis and the translational invariance along the x and y coordinate axis. Finally 3N degrees of freedom are found for the system in interaction with a corrugated substrate, since in this case all spatial symmetries are broken. Force evaluation After the initialisation the algorithm enters the main part of the simulation which consists of a loop over the force evaluation and the numerical integration of the equations of motion. The forces on every atom are computed with the help of the potentials presented in section 2.1. There are in fact two ways of achieving this: The first method is to differentiate the expression for the total potential analytically: ~ i Vtot f~i = −∇ (12) The expressions are derived in appendices A and B. The analytically determined forces are exact and fast to evaluate. However, this involves quite lengthy expressions, so the second method is used as a valuable test to check whether the implementation of the analytic derivation is correct. This consists of differentiating these potentials numerically. This can be achieved by evaluating the total potential energy Vtot , which depends on all the position vectors of the system of N atoms, at two virtual configurations for each component of the force vector of every atom: one where the α-component of the position vector of atom i in question has been increased by δ, which represents a small distance, and the other when it has been decreased by δ. In theory, the computed force becomes more accurate as δ approaches zero, but this is limited by the numerical precision: fiα = − Vtot (r1x , r1y ...ri,α + δ, ...rN,z ) − Vtot (r1x , r1y ...ri,α − δ, ...rN,z ) 2δ (13) Integration The next step is to integrate Newtons equation of motion with respect to the time t: d2 r~i f~i = mi · 2 dt (14) For this purpose the velocity Verlet integrator is employed, which is discussed in section 2.2.2. In the case of a simulation in the canonical ensemble, the situation is a bit more complicated since a thermostat comes into play. For the simulations presented in this report the choice fell on the NoséHoover thermostat which is deterministic and easier to check. The corresponding extended equations of motion as well as the fact that its usage corresponds to sampling the canonical ensemble is shown in section 2.2.3 and appendix C. Throughout the simulation a variety of quantities of interest can be calculated. One great advantage of molecular dynamics simulations is that it is possible to calculate every quantity that can be expressed as a function of the microscopic variables. When interested in equilibrium properties that do not depend on time, the statistics for the quantity in question becomes better as the simulation time increases. A 5 Modelling of the structure and dynamics of nanoparticles on metal substrates Daniel Förster major change in the average of such a quantity in the course of the simulation would indicate that the system is not in equilibrium. The computation of dynamical properties is more complicated, statistics on the considered quantity can be obtained by carrying out many simulations of the same type. Another restriction is that the presence of a thermostat affects the dynamics of the system in an uncontrolled way, therefore it is somewhat unsafe to compute dynamical properties in the canonical ensemble. 2.2.2 Velocity Verlet integrator The integration of Newton’s equations of motion in the simulations is achieved by means of the so-called velocity Verlet scheme. It is not only simple and fast, but it also displays long-term energy conservation. The latter is important, as this property ensures that the system stays in the intended microcanonical ensemble. Higher-order schemes may guarantee a higher short time precision, but often they produce a long-term energy drift, which makes them less interesting for molecular dynamics. Another interesting property is that the generated trajectories are time-reversible (symplectic). This holds of course only within the numerical precision, but it is independent of the length of the time step [5]. The velocity Verlet integrator is given by the following set of equations, where v~i designates the velocity vector of atom i and dt the timestep. r~i (t + δt) = r~i (t) + v~i (t)δt + v~i (t + δt) = v~i (t) + f~i (t) 2 δt 2mi f~i (t + δt) + f~i (t) δt 2mi (15) (16) The error is O(δt4 ) for both the positions and the velocities. 2.2.3 Nosé-Hoover thermostat The Nosé-Hoover thermostat algorithm allows to drive the system to an imposed temperature and maintain it there. More importantly, this type of thermostat generates configurations that follow the canonical distribution [5]. In contrast to the Andersen thermostat that also may reproduce correctly canonical statistics, the Nosé-Hoover thermostat is time-reversible and deterministic [5], i.e. it does not rely on random changes to the system. Instead, an extra variable with a coordinate s and momentum ps model the heat bath. Its interaction with the physical system is computed with the help of a specially designed potential that allows the system to stay within canonical statistics. The NoséHoover thermostat features one parameter which corresponds to the its inertia Q. This parameter can also be seen as the coupling strength to the heat bath. The numerical value of this coupling parameter is important for the efficiency of the thermostat and must therefore be chosen with care which is discussed in section 3.1.3. The entire system including the thermostat is again found to sample a microcanonical ensemble. This fact allows the definition of a conserved quantity, namely the total energy including the thermostat. This quantity can therefore be used to validate the simulation code. It should be noted that also in the canonical ensemble temperature fluctuations occur for finite systems. Therefore, in contrast to simple velocity rescaling thermostats, the instantaneous temperature of a system in contact with the Nosé-Hoover thermostat fluctuates. These fluctuations follow, however, the real canonical fluctuations. This property of the thermostat is shown in appendix C. The equations of motion for the system in contact with the Nosé-Hoover thermostat read [6]: 6 Modelling of the structure and dynamics of nanoparticles on metal substrates d~ ri p~i = dt mi d~ pi ∇i Vtot ps p~i =− − dt mi Q 2 X p~i dps = − f kB T dt mi Daniel Förster (17) (18) (19) i ds ps = dt Q (20) The target temperature of the thermostat is designated with T , which appears together with the Boltzmann constant kB and the number of degrees of freedom of the system as a whole f . The coordinates of atom i of mass mi are r~i and their corresponding momenta p~i , while Vtot corresponds to the total potential for the physical system (without the thermostat). The actual scheme for this thermostat that has been used in the simulation code is taken from [7], which gives the velocity Verlet integrator of the equations of motion shown above. 3 Implementation details and observables 3.1 3.1.1 Adjustment of simulation parameters Interaction with the substrate and coordination dependence Continuous substrate potential The previous sections explained the fundamentals on which the simulations of this work are based. Before being able to actually run simulations, it is however necessary to choose the parameters for the potentials, as discussed here. The parameters of the Sutton-Chen potential for iridium did not have to be adjusted, since they could be adopted from [3]. In principle, the system is already described by this potential, since the nanoparticles as well as the substrate contain iridium atoms only. However, as mentioned in section 2.1.2, the interaction with the substrate has been treated with a continuous potential containing the periodicity of the substrate. This potential has a number of parameters that have been adjusted to the potential as is given by the Sutton-Chen potential in a fully atomistic approach. The first step towards the determination of the parameters is the computation of the Sutton-Chen potential of a single atom (coordination zero) as a function of its distance from the substrate. This has been carried out at different positions of the triangular unit cell of the surface lattice. These positions and the corresponding values of fi,1 and fi,2 are as follows: √ • "Centre": Centre of the equilateral triangle (e.g. x = b/2, y = b/(2 3)), fi,1 = −1.5, fi,2 = 3 • "Bridge": Middle of an edge of the triangular unit cell (e.g. x = b/2, y = 0), fi,1 = −1, fi,2 = −1 • "Top": Edge of the triangular unit cell (e.g. x = 0, y = 0), fi,1 = 3, fi,2 = 3 • "Average": Average value over the entire unit cell, fi,1 = 0, fi,2 = 0 which is also illustrated in Figure 2. 7 Modelling of the structure and dynamics of nanoparticles on metal substrates Daniel Förster Bridge Centre Top Average b y x Figure 2: Considered position in the unit cell The potential has also been computed for a dimer (coordination one), a trimer (coordination two), a partial layer (coordination three) and a complete layer (coordination six) of the (111) surface. The data from these calculations have then been used to adjust the parameters of the potential. The potential energy data from the Sutton-Chen potential and the potential curves with the adjusted parameters are shown in Figure 3. The panel on the left shows the potential energy for a single atom (coordination zero) in dependence of the distance from the substrate at the four considered positions of the surface lattice unit cell. The right panel shows the same thing, but for an atom with coordination six. In appendix D the corresponding curves for coordination one, two and three can be found. The fitting procedure relies on two steps: First, with the help of the potential curves at the "average position", the parameters for the planar substrate can be adjusted, that are the 12 parameters that occur in the zeroth order term of VSubstrate , where l = 0. In a second step, the numerical values for the remaining 24 parameters are obtained by a global adjustment of these parameters for the full potential (containing the substrate periodicity) leaving the parameters determined in the first step unaltered. This way a total of 36 parameters has been computed; their numerical values can be found in appendix D. This fitting procedure was subject to previous work, it was not carried out in the scope of this work. (a) Single atom (b) Monolayer 2 centre bridge top average 0 Binding energy [eV/atom] Binding energy [eV/atom] 2 -2 -4 -6 1 2 3 4 Distance to surface [Å] 5 0 -2 -4 -6 1 6 centre bridge top average 2 3 4 Distance to surface [Å] 5 6 Figure 3: Sutton-Chen potential energy as a function of the distance to substrate at different positions within the triangular unit cell of the surface lattice. The continuous lines show the result of the curve fitting procedure. The right panel shows the potential energy for a single atom (c = 0) whereas the left panel shows the potential energy for an atom with coordination c = 6. The corresponding curves for the other coordination numbers are shown in appendix D. Determination of the coordination The formula that gives the numerical value of the coordination c of a certain atom has two variables rmin and rmax , as described in section 2.1.2. Since the coordination should give the number of neighbours of an atom, the choice was made in a way that the nearest neighbours of the atoms within the nanoparticles are counted entirely and that the second neighbours are not counted anymore. This means rmin needs to be larger than the nearest neighbour distance and rmax needs to be smaller than the distance to the second neighbours. Additionally, it is 8 Modelling of the structure and dynamics of nanoparticles on metal substrates Daniel Förster preferable to have a low difference between rmin and rmax . It should however be large enough to ensure that the potential does not become too steep, which would increase numerical errors. The final choice for the parameter pair was: rmin = 2.9 Å and rmax = 3.8 Å. During the evaluation of the forces in the molecular dynamics code, the contribution of the substrate atoms to the coordination is also taken into account, as mentioned in section 2.1.2. The substrate atoms however do not appear explicitly in the model which makes it necessary to estimate their contribution to the coordination of an atom within the nanoparticle. An atomistic description of the substrate is used in order to be able to calculate the real physical coordination of an atom as a function of the distance from the substrate. Then a seventh order polynomial is adjusted to these data in order to give an approximation of the additional coordination of an atom due to the presence of the substrate. For distances that exceed rmax there is no contribution of the substrate to the coordination number. It is imposed at rmax that the polynomial is zero as well as its first derivative. In order to fix the remaining five degrees of freedom, a curve fitting procedure is used. In reality, the coordination depends, however, not only on the z coordinate but also on the horizontal position within the unit cell of the surface lattice. Therefore, data from 500 random positions have been used to adjust the free coefficients of the polynomial. The data as well as the adjusted polynomial are shown in Figure 4. 10 computed coordination adjusted polynomial Coordination 8 6 4 2 0 0 0.5 1 1.5 2 2.5 3 3.5 4 Distance from surface [Å] 4.5 5 Figure 4: Coordination of a single atom due to the presence of explicit substrate atoms as a function of distance to surface. The Figure shows mean values over 500 random (x,y) positions, while the error bars correspond to the resulting standard deviation. A seventh degree polynomial has been adjusted to the data (not all points shown) imposing the polynomial as well as its first derivative to vanish at a distance of rmax . 3.1.2 Stability of trajectories The quality of energy conservation during the simulation depends on the choice of the time step that is used for the numerical integration. It should be chosen as high as possible in order to save computation time, but as low as necessary for an adequate energy conservation. The stability of the total energy, or more precisely its standard deviation during a simulation of 128 ps has been studied. The simulations have been carried out in the microcanonical as well as the canonical ensemble. For the simulations in the canonical ensemble the total energy includes also the energy of the thermostat. The results of this study are shown in Figure 5. In the simulations in the microcanonical ensemble the standard deviation of the total energy becomes smaller then the numerical precision for simulations with a time step smaller than 1 fs. The same happens in the canonical ensemble, but only for time steps smaller then 0.25 fs. It can be noted that the energy conservation is slightly better in all simulations when 9 Modelling of the structure and dynamics of nanoparticles on metal substrates Daniel Förster Standard deviation total energy [%] carried out in the microcanonical ensemble. For all the latter simulations a time step of 1 fs has been chosen. microcanonical ensemble canonical ensemble 0.1 0.01 0.001 0.0001 1e-05 1e-06 0.1 1 10 Time step [fs] 100 Figure 5: Dependence of the standard deviation of the total energy on the length of the time step dt in microcanonical and canonical simulations 3.1.3 Choosing the mass of the Nosé-Hoover thermostat The Q parameter of the Nosé-Hoover thermostat must be chosen with care. If the value of this pseudomass of the thermostat is chosen too big, the convergence to the target temperature is slow. In the limit of an infinite Q the system behaves as without the thermostat i.e. follows the microcanonical ensemble. On the other hand, however, for low values of Q, thermal energy may flow in and out of the thermostat periodically, which leads also to a slow convergence to the target temperature. It is recommended to choose a numerical value according to: Q≈ f kB T ω2 (21) where ω is a typical system frequency [6]. The typical frequency of the nanoparticles is determined in dedicated simulations, where the initial velocities are not distributed randomly, but all the velocity vectors pointing to the centre of the nanoparticle. This kind of excitation leads to periodic breathing motion of the nanoparticle of which the frequency has been recorded, as shown in Figure 6. 10 Modelling of the structure and dynamics of nanoparticles on metal substrates Frequency [THz] 12 Daniel Förster 53-atoms nanoparticle 13-atoms nanoparticle dimer 10 8 6 4 2 0 0.0001 0.001 0.01 0.1 1 10 Kinetic energy per atom [eV] 100 Figure 6: Determination of a typical system frequency for Ir nanoparticles of different sizes as a function of the imposed total energy One observes that the frequency of this breathing oscillation decreases with growing energy, which is due to the anharmonicity of the potential. Furthermore this study shows that this type of oscillation has a lower frequency for the larger nanoparticles. More interesting for the latter simulations is the determination of a typical system frequency. The choice fell on f ≈ 4.5 THz, which corresponds to a wavenumber of ω = 150 cm−1 . This choice was partly motivated also by a subsequent series of tests. The stability of the system temperature when the thermostat is applied was therefore examined. Simulations at different values of the parameter ω were run for the Ir55 nanoparticle over a duration of 100 ps at 400 K. The result of this series of tests is shown in Figure 7. It appears that the choice of ω = 150 cm−1 is not the value with the lowest temperature fluctuations. Small temperature fluctuations, p however, are expected also in the canonical ensemble. They are in the order of 2/(3N ), where N is the number of atoms. For 55 atoms at 400 K this would lead to fluctuations with a standard deviation of about 44 K, which is close to the value observed in this study for the chosen value of ω. Standard deviation [K] 1000 temperature stability 100 10 1 10 100 1000 10000 100000 -1 ω [cm ] Figure 7: Stability of temperature in presence of the Nosé-Hoover thermostat as a function of its parameter ω 11 Modelling of the structure and dynamics of nanoparticles on metal substrates 3.2 3.2.1 Daniel Förster Observables Stable structures All simulations presented in this report are carried out in equilibrium with the exception of the one about deformation dynamics, shown in section 4.2. Therefore, the starting point has to be an equilibrium structure of the nanoparticle. In the case of simulations with free nanoparticles, the equilibrium structure had already been determined. For the simulations with interaction with the substrate, the equilibrium structures had first to be found. In order to do so, a slightly changed version of the simulation code has been used for determining a local energy minimum by simply following down the energy gradient along the so-called steepest descent path. The structure in the local minimum has then been used in simulations like the one discussed in 4.2, where this local minimum structure is heated up and then from time to time a configuration of the nanoparticle is taken out and relaxed to its local energy minimum structure. Among a total set of 25 structures obtained this way, the one with the lowest potential energy has been used as initial configuration for equilibrium simulations. This energy minimization procedure has been followed for the three nanoparticles (Ir13 , Ir55 , Ir147 ) in the case of nanoparticles in the presence of the non-corrugated and corrugated substrates. The resulting structures are presented in section 4.1. 3.2.2 Thermodynamical properties and Lindeman index Along the molecular dynamics simulations it is possible to extract quantities that are based on the microscopic variables of the system. The thermodynamical variables considered in this work are the internal energy U and its derivative, the heat capacity Cv . The precision increases for longer simulation runs, because these quantities are based on the calculation of averages. The internal energy U is just the time average over the potential energy E of the system by virtue of the ergodicity principle. In the following, the time average is denoted by h·i: U = hEi (22) The heat capacity Cv can be determined by the following relation, which necessitates to sum up also E 2 during the simulation runs [5]:   hE 2 i − hEi2 ∂U Cv = = (23) kB T 2 ∂T During the initialization process, the atoms have been assigned random velocities. It is, therefore, necessary to let run the simulation for a certain amount of time until the velocity distribution follows the Boltzmann statistics. In the following, it is assumed that the system has reached equilibrium after 1 ps or a thousand time steps. This is long enough as it has been verified by observing the evolution of the velocity distribution. After equilibration, it is then possible to start summing up the variables in order to obtain the mentioned averages. Another quantity of interest is the Lindemann index δ, which is used to quantify the degree of rigidity of the atoms within a nanoparticle of N atoms: q 2 i − hr i2 X hrij ij 2 δ= (24) N (N − 1) hrij i i