DE Original scientific paper Received: February 18, 2016 Accepted: March 24, 2016 DOI: 10.1515/rmzmag-2016-0002 Molecular dynamics simulation of aluminium melting Simulacija taljenja aluminija z uporabo metode molekularne dinamike Jakob Novak Faculty of Mathematics and Physics, University of Ljubljana, Jadranska ul. 19, 1000 Ljubljana, Slovenia novak.jakob@yahoo.com Abstract Solid-liquid phase transition has been simulated by the molecular dynamics method, using isobaric-iso-enthalpic ensemble. For interatomic potential, glue potential has been selected. The original algorithm for bookkeeping of the information on neighbouring relationships of the atoms has been developed and used in this research. Time consumption for calculation of interatomic forces has been reduced from o(N2) to o(N) by the use of this algorithm. Calculations show that phase transition from solid to liquid occurs between 1,000 K and 1,300 K. The simulated temperature of phase transition is higher than the experimental value due to the absence of crystal defects. If constant heat flux is supplied, temperature decreases during melting because the superheated state becomes unstable. During the cooling process, no significant changes of the observed variables were detected due to the high cooling rate, which prevents crystallisation. Key words: molecular dynamics simulation, isobaric-isoenthalpic ensemble, glue potential, bulk aluminium melting Izvleček Fazni prehod trdno-tekoče je bil simuliran z metodo molekularne dinamike, z uporabo izobarnega-izoental-pnega ansambla. Za medatomski potencial je bil uporabljen večdelčni potencial, ki ga je razvil F. Ercolessi (ang. Glue potential) [1]. Originalen algoritem, za shranjevanje informacije o sosedskih odnosih med atomi, je bil razvit in uporabljen v pričujoči raziskavi. Časovna poraba izračuna medatomskih sil je z njegovo pomočjo zmanjšana iz o(N2) na o(N). Izračuni kažejo, da dobimo fazni prehod iz trdnega v tekoče, med 1000 K in 1300 K. Simulirana temperatura faznega prehoda je višja od eksperimentalne vrednosti, zaradi odsotnosti defektov v kristalu. Če dovajamo konstanten toplotni tok, med taljenjem temperatura vzorca pade, saj postane pregreto stanje nestabilno. Pri ohlajanju taline, ni bilo zaznanih signifikantnih sprememb v opazovanih spremenljivkah, zaradi visokih hitrosti ohlajanja, ki preprečijo kristalizacijo. Ključne besede: metoda molekularne dinamike, izo-baren-izoentalpen ansambel, večdelčni medatomski potencial, taljenje razsežnega aluminija 9 na iimi-am © 2016 Jakob Novak, published by De Gruyter Open. 10 Introduction Melting and solidification constitute an important segment of metal alloy production. The process of phase transition is influenced by many parameters. One of them is the interatomic potential energy, which is the consequence of previous processing of the material. This study is an introduction to research into the influence of crystal defects on the characteristics of solid-liquid phase transition of metals. Our special point of interest is the influence of defects on the enthalpy of fusion and short-range order in the melt. To observe the influence of crystal defects on the behaviour of a sample, simulation, at the atomic level, is required. Two different methods are commonly used for simulating the behaviour of materials, on the basis of their atomic configuration: molecular dynamics and Monte Carlo simulation [2]. Both methods extract macroscopic quantities from the time evolution of particles on the microscopic level. The difference is that in molecular dynamics, time evolution is deterministic, in contrast to Monte Carlo simulations, wherein atomic states are chosen randomly for every time step, according to suitable distributions. An example of kinetic Monte Carlo simulation considering asymmetric diffusion in Fe-Cu alloy is presented in the study by Bombac and Kugler [3]. In this research, the molecular dynamics method has been chosen, because it automatically includes the influence of sample history on the behaviour of the material, which is much harder to incorporate in Monte Carlo simulations. This is very important for the study of the influence of crystal defects on the phase transition. In literature, molecular dynamics simulations are used for the examination of different processes in materials, e.g., for calculation of transport and equilibrium diffusivity in Lennard-Jones fluids [4] and for investigation of surface diffusion [5]. Molecular dynamics simulation has also been widely used for examination of the behaviour of crystal defects. In previous studies [6, 7], migration of grain boundaries has been simulated; some studies [8, 9] discuss the dislocation processes as a consequence of mechanical deformation in aluminium. In this study, simulation of aluminium phase transition is presented. Crystal defects have not been applied yet. The sample is an ideal, infinite, face-centred cubic (FCC] crystal. The infinite size of the crystal has been accomplished by periodic boundary conditions. The reason for application of periodic boundary conditions is that the current computational capabilities typically enable molecular dynamics simulations of samples of sizes up to 10 nm. This type of size is satisfactory for simulating the behaviour of nanoparticles, but it does not satisfy for observation of macroscopic sample, so called bulk aluminium. Even though we are dealing with the same element in the case of nanoparticle and bulk aluminium, there are important differences in the properties of the two samples [10]. The latter has no free surface. This results in different physical properties of phase transition because the free surface acts as a nucleation site for phase transition [11]. Theoretical introduction In this study, the isobaric-isoenthalpic ensemble (NPH) has been selected because it allows energy and volume to fluctuate. This is not the case in a microcanonical ensemble (NVE], which uses Newtonian equations of motion and periodic boundary conditions. Energy and volume fluctuations are very important for the description of phase transition because energy fluctuations and thermal expansion affect the melting point. The isobaric-isoenthalpic ensemble uses volume as the additional degree of freedom to describe the coupling of the observed sample with its surroundings. The equations of motion for N atoms are expressed as follows [12]: i 2 i \ MV=-P+— W and > (12) where ^=(0,0,0), a2=(0,a,0) and a3 = (0,0, fl) are primitive vectors of the cubic Bravais lattice and nx, n2 and n3 can be any integer. For the FCC lattice, at every point of the Bravais lattice, a basis consisting of four atoms is placed. The positions of atoms are therefore determined by the Bravais lattice vector and one of the four basis vectors: r0 («J,w2,n3,k) =rB (n1,n2,n3)+bk k = 1,2,3,4 (13] 5 =(0,0,0), b2=^(a,a, 0), Vector R can be expressed as small fluctuations E = —,eN) around the minimum: R — Rq + E. After the Taylor expansion, Eq. (11] gets the following form: dw = A' exp(-(Zij -Ji-^. (*o) h + + I,ilmivf)/kBT)dEdV, (14] Velocities of particles are the independent degrees of freedom because the probability for a particle to have certain velocity is independent of the velocities and positions of other particles. But that is not the case for positions. For that reason, transition to normal coordinates, q., is necessary. Normal coordinates are the normal modes of the crystal. They are obtained by diagonalisation of the potential U. Amplitude of the normal mode, qi = becomes the independent degree of freedom and probability for it can be chosen randomly, according to probability distribution [14]: dw{qi) = ,4'exp (-mrfqf / 2kj) dqt , (15] where fti». is the eigenfrequency of the crystal, calculated from the eigenvalues of the matrix of second derivatives of the potential energy. Velocities can be determined in the standard coordinates dw(v,.) = ^"exp (-mtf / 2kBT) dvt. (16] Eqs. (15] and (16] are useful only if we know the minimum vector R„. In the case of ideal FCC lattice, we know the global minimum from experiments; however, if crystal defects are present in the solid, it is much harder to determine the exact vector of the local minimum. That is why this research has been carried out to investigate how fast atoms equilibrate. The conclusion is that there is no need to set correct initial conditions, but an arbitrary initial condition can be selected, after which atoms spontaneously transit into the equilibrium state on a time scale shorter in comparison to the simulation running time. Computational framework All the results presented in this study were computed with the code written especially for the needs of this research, and no commercial program for molecular dynamics has been used. In this section, the algorithms constituting the simulation are presented, with emphasis on the algorithm for bookkeeping of neighbouring relationships of the atoms, which has been developed for the needs of this study and differs from the commonly used Verlet list or Linked list. The sample observed was a cube, with initial side length of 3.24 nm, consisting of 2,048 atoms. Vectors of Bravais lattice were parallel to the cube sides. To simulate the melting of bulk aluminium, periodic boundary conditions were applied. When an atom passes through the face of a cube, it enters through the opposite face of the cube. Atoms located at the surface of the cube also feel the potential of the atoms from the opposite side. With this procedure, an effectively infinite crystal is generated. The system of second-order differential equations (1] and (2] has been solved numerically by the Gear predictor corrector algorithm [15]. This algorithm is especially appropriate for systems with large numbers of particles because it is not time consuming and is a very good energy conserver. The time step has been chosen to be 1 fs, considering that the time scale for vibration of atoms is of the same order [11]. During the heating process, constant heat flux is supplied to the sample. At every time step, the average energy per particle is increased by 2.6x10-6eV, by scaling the velocities. The same heat flux, with a negative sign, is set for the cooling process. Table 1 shows the list of all parameters used to obtain the results presented in this paper. The most time-consuming part of the simulation is the computation of forces because the total force on the selected particle is the sum of forces due to all the particles in the simulation. In order to shorten the computational time, only the forces due to the surrounding atoms Table 1: List of all parameters used in the simulation, with corresponding values. Parameter name Parameter value Number of atoms (N) 2,048 Aluminium lattice parameter (a) 0.4046 nm Mass of aluminium atom (m) 25.15 GeV/c2 Fictitious mass (M) 610 keV/c2 Integration time step 1 fs Potential cutoff distance r c 0.555 nm Pressure (P) 0 Pa Heat flux 2.6x10-6 eV/fs are calculated; other forces are assumed to be negligible. To reduce time consumption from o(N2) to o(N), the algorithm that retains the information of neighbouring relationships has been applied, to avoid calculation of the distances due to all atoms. In literature, the most commonly used algorithms for the bookkeeping of neighbouring relationships are the Verlet list and the Linked list. The first one reduces the time consumption to o(N3/2) [2] and is the most efficient for smaller samples. The second reduces the time consumption to o(N) [2], but it is less efficient than the algorithm presented in this study for most of the situations. The basic principle of the algorithm is the following. Before the simulation starts, every atom is given a unique number, which serves as its identity. Two lists are constructed for every atom. The first one contains the numbers of all atoms situated in a sphere of radius rc (herein referred to as 'near neighbours'], where rc is an arbitrary selected value that represents the range of interatomic interaction. Because we are dealing with a numerically computed potential, rc coincides with the potential cutoff distance. In this study, the cutoff has been set to 0.555 nm because at this distance, the functions