The iron melting point determination by 2D simulation

In the past research, molecular simulation has been widely used to simulate various properties of materials and understand atomic theory. The effect of the material’s surface tension and different closed-packed structures can also be successfully simulated by using the Lennard-Jones potential. However, the two-dimensional simulation of the iron’s melting point has not been studied deeply. Thus, by MATLAB, this project used the Lennard-Jones potential and 2D simulation to simulate the melting point of the iron (698.82 K-725.08 K) under certain pressure. This pair coefficient equation of the iron atom from 300K to 900K was found successfully. The pair coefficient equation is generally in line with expectations, but the calculated melting point is far away from the experimental result (1811 K). This is mainly because of the limitation of 2D simulation. (the vertical force is not contained) The advice for the future study is to use the 3D simulation to do the more accurate iron’s melting point calculation.


Introduction
Molecular interaction simulation is a wildly used simulation process to explain and determine some chemical and physical properties of materials, and molecular basis theories, such as permeation, binding, folding, conformational changes [1,2].For example, by using molecular dynamic (MD) simulation, Urooj Fatima successfully predicted the excess volume (V E ) and viscosity deviations (∆) of the ionic liquid: 3-methylimidazolium bis(trifluoromethylsufonyl)imide + 1-butanol/1-propanol mixtures [1].The results of experiment and MD simulation are basically consistent.During the calculation, molecular simulation considers the effect of entropy [3].The Lennard-Jones potential is used for describing how the intermolecular interaction changes with distance [4].It can be used to study the effect of surface tension and different closed-packed structure including: sc, bcc, and fcc [5,6].The two-dimensional simulation of the metal's melting point has not been discussed.
This project uses MATLAB to simulate the melting point of iron atoms in two dimensions.The intermolecular force is gained from the Lennard-Jones potential.The motion of iron atoms is simulated by the modified forward Euler method.Ar atom is introduced to the system to simulate atmospheric pressure.
The purpose of this project is to find out the melting point of the iron and pair correlation function at each temperature by using 2D simulation.

Theory
As one of the most popular intermolecular interactions between particles, the Lennard-Jones potential is used in this research project for simulating interaction among iron and argon particles.The Lennard-Jones potential function has the general form of: (1) Where,  is the distance when the potential energy is zero. corresponds to the bond energy between particles. means the distance between two particles.
The minimum of the Lennard-Jones potential is at the distance of   = 2 1 6 .For the convenience of MATLAB simulation, the equation ( 1) is changed to AB form I n the code, which is: (2) Here  = 4 12 and  = 4 6 .For simulating the interaction between two different types of atoms (for example, particle a and particle b), the parameter   and   used in the Lennard-Jones potential can be calculated from [7]: = √  ×   (4) Where   and   are used to describe intermolecular interactions between particles a.   and   are used to describe interaction between particles b.
For simulating the motion, the force form of the Lennard-Jones potential function in 2D is needed, which is: Where x and y mean the particles' distance along x axis and y axis respectively.  and   are the forces on the particle in the x and y directions respectively.
The accelerations calculated along x axis and y axis are: =    (8) Here, m is mass of the particle.The mass of a single particle can be calculated from the molecular mass of it: Where, M is molar mass of the particle corresponding to the value shown on the periodic table.  is Avogadro constant with the value of   = 6.022 × 10 23  −1 .
For simulating the motion process for every two steps, the modified forward Euler method is used.A half step calculation is used for reducing the error, shown in the equations (11) and (12) below: Where  is the length of each time step. is the starting time of the first step.  and   are the velocities of particles along x axis and y axis respectively.
In this program, the temperature of the system is determined by the average velocity of the particles, which is: Here,   means the average kinetic energy of all particles.v is the total velocity of the single particle.( = √  2 +   2 ) is the dimensions of the system.Since this project was simulated in 2D system, the value of is expected to be 2. k is Boltzmann's constant with the value of  = 1.38 × 10 −23 /.
For cooling or heating particles in the system, the velocity of the previous step is multiplied by a coefficient (the modified form of the equation ( 12)): ( + )) × 0.9999 (15) Where the coefficients 1.0001 and 0.9999 are used for heating and cooling respectively.

Method
This program is simulated in MATLAB with the using of the theories explained above.Here are two types of atoms simulated in this program.One type is the inert gas argon (Ar), which is used to simulate the air pressure.Another is this research object iron (Fe).According to the phase diagram of iron, iron will only undergo a phase transition from solid to liquid if there is a certain pressure [8].Therefore, argon particles were added to the simulation of the system to simulate the pressure so that the liquid state could be observed.
The Lennard-Jones potential function is used to simulate the interaction among particles.The parameters in equation ( 1 All particles are simulated in a closed box system with a length of 100  .For the initial setup, here are 50 Ar atoms, which are evenly distributed in the system, and 100 Iron particles separated by a length of   in order to get an initial solid state.The particles in the system are assigned random velocities based on a certain initial temperature.Due to the mutual conversion of kinetic energy and potential energy during the initial running of the code, the system needs to run for a few time steps (200 steps) before it can be in the most stable state.Due to the interconversion of kinetic energy and potential energy, even in the absence of heating or cooling, the temperature of the particles will be different from the initial temperature set after the system has been running for some time.For controlling the temperature of the particles, formulas (14) and (15) are used for heating and cooling, respectively.
Two groups of experiments were carried out in this project.One was to measure the melting point of iron by heating from low temperature to high temperature (from 200 K to 1050 K), and the another was to measure the melting point of iron by cooling from high temperature to low temperature (from 900K to 530 K).In both sets of simulations, 40,000 steps were performed respectively.The time step length in the simulation is determined by the initial particle maximum, following the equation:  =  5000•  .This setting of the time step length can prevent large error caused by the motion simulating process and a lot of computing caused by too short time step length.
In simulating the particles' trajectory, the modified forward Euler method is used (by using equations ( 10) -(12) in the code simulation).For the calculation of acceleration, the cut-off range is chosen to be 10  for avoiding excessive calculation.The selection of this cut-off range means that only the intermolecular force with a distance of less than 10  is included in the acceleration.Because the intermolecular potential used in this simulation is the Lennard-Jones potential, the intermolecular force is very small when the distance between molecules is greater than 10  .Therefore, the forces with the greater than 10  distance between molecules can be ignored, and the resulting error can be so small as to be negligible.The acceleration of a particle is determined by the combination of forces within the cutoff range.In each time step, formulas (10) -( 12) are calculated, and new accelerations are calculated based on the particle position and the Lennard-Jones potential.In addition, in each time step, formulas (14) or (15) are also used according to heating or cooling requirements.
After the simulation of the motion trajectories in each time step, the temperature and potential energy of iron particles can also be determined according to the distance between particles and velocity of particles by equations ( 2) and (13).The pair coefficient function can also be determined by checking how many particles are in different distance ranges for each particle.
After the completion of the two groups of simulation experiments, the motion trajectory, temperature, potential energy and pair coefficient function data of iron particles at different temperatures were obtained.The data begins to be analysed and processed, and more details will be shown in part 4.

The relationship between the particles' temperature and potential energy
In the process of simulation, the potential energy and kinetic energy of the whole system keep converting to each other.This leads to very large fluctuations in measured temperature and potential energy.In order to make it easier to show specific trend of change from the plot of data, the data used in the plot has been averaged.The value of the potential energy is the average of the potential energy every 200time step.The kinetic energy is the average value per thousand-time step.By averaging the data, a smooth plot of the relationship between particle temperature and potential energy is obtained (as shown in the Figures 1-2).According to the Figure 1 and 2, it is obvious that the slope of the curve changes significantly after a particular temperature.The change in slope is believed to be due to a phase transition in the metal.This particular temperature point is considered as the melting point.
In order to obtain the melting point, two best fitting straight lines can be obtained respectively for the changing trend before and after the melting point.The intersection of these two lines is the melting point in the two-dimensional simulation.The two best fitting lines in the cooling process are respectively: (19) By combining the equations ( 18) and (19) above, the melting point is gained to be 725.08K for the heating process.
The reason why this program gets two different melting points is because of overheating and overcooling.Before the system reaches equilibrium, the next time step of heating or cooling is carried out.So that the actual dramatic change occurs either before or after the measured melting point.Therefore, the expected melting point of this system is 698.82K -725.08 K.
Here are two reasons why the melting point obtained by the 2D simulation is much lower than the theoretical value of 1811 K [10].One reason is that the pressure in the theoretical simulations does not reach an atmosphere pressure needed to restrain the metal.Thus, the measured melting point will be smaller according to the phase diagram.Another reason is the limitations of two-dimensional simulations.Since the two-dimensional simulation only takes into account intermolecular forces within the same plane (the x-y plane), forces from the vertical direction (the z-axis) are not taken into account.Due to the absence of vertical forces, the constraint effect of intermolecular forces on the solid state of particles decreases exponentially.Thus, the simulated resulting melting point is expected to be much lower than the experimental value.With the increase of the temperature, the pair correlation value gradually becomes lower and lower.One of the reasons is the flattening of peaks.In order to integrate to obtain the same number, the highest value of the more widely distributed peak is smaller.Another reason for it is the transition to the gas phase.Due to the limitations of the two-dimensional simulation, the forces between the particles to maintain the solid and liquid state are not strong (due to the lack of vertical forces).Thus, as the temperature rises, the particles tend to disperse and some become gas phase.

Conclusion
Through two-dimensional simulation, this project finds that the melting point of iron is 698.82K -725.08K under certain pressure and the pair coefficient equations under different temperature conditions are successfully obtained.The found melting point is much smaller than the experimental melting point (1811 K).The pair coefficient equation found is generally consistent with the expectation from the solid state to the liquid state.The reason why the found melting point does not fit the experimental value is the limitation of the two-dimensional simulation.The vertical interactions are not considered during the simulation, so some physical properties such as melting point should be far away from the actual value.Another reason is that the simulated pressure does not reach the standard atmospheric pressure conditions, so the pressure does not bind the iron atoms enough to keep them in a solid state at higher temperatures.
Due to the limitation in computing power, the 2D simulation is used in the simulation.In the future study, the 3D simulation can be used for simulating more accurate physical properties of the metal atoms.The study simulates the properties of one kind metal atom (iron Fe).The similarly simulation method can then be used to investigate the chemical and physical properties of different metal atoms.

Figure 1 .
Figure 1.The relationship between the particles' temperature and potential energy for the cooling process.

Figure 2 .
Figure 2. The relationship between the particles' temperature and potential energy for the heating process.According to the Figure1and 2, it is obvious that the slope of the curve changes significantly after a particular temperature.The change in slope is believed to be due to a phase transition in the metal.This particular temperature point is considered as the melting point.In order to obtain the melting point, two best fitting straight lines can be obtained respectively for the changing trend before and after the melting point.The intersection of these two lines is the melting point in the two-dimensional simulation.The two best fitting lines in the cooling process are respectively:  = 1.8600 × 10 −21  − 9.0656 × 10 −18 (16)  = 6.1575 × 10 −21  − 1.2069 × 10 −17 (17) By combining the equations (16) and (17), the melting point calculated for the cooling process is 698.82K.The two best fit lines in the heating process are respectively:  = 1.8957 × 10 −21  − 9.3335 × 10 −18 (18)  = 1.3439 × 10 −20  − 11.7703 × 10 −17(19) By combining the equations (18) and (19) above, the melting point is gained to be 725.08K for the heating process.The reason why this program gets two different melting points is because of overheating and overcooling.Before the system reaches equilibrium, the next time step of heating or cooling is carried

𝑈 = 1 .
Figure 2. The relationship between the particles' temperature and potential energy for the heating process.According to the Figure1and 2, it is obvious that the slope of the curve changes significantly after a particular temperature.The change in slope is believed to be due to a phase transition in the metal.This particular temperature point is considered as the melting point.In order to obtain the melting point, two best fitting straight lines can be obtained respectively for the changing trend before and after the melting point.The intersection of these two lines is the melting point in the two-dimensional simulation.The two best fitting lines in the cooling process are respectively:  = 1.8600 × 10 −21  − 9.0656 × 10 −18 (16)  = 6.1575 × 10 −21  − 1.2069 × 10 −17 (17) By combining the equations (16) and (17), the melting point calculated for the cooling process is 698.82K.The two best fit lines in the heating process are respectively:  = 1.8957 × 10 −21  − 9.3335 × 10 −18 (18)  = 1.3439 × 10 −20  − 11.7703 × 10 −17(19) By combining the equations (18) and (19) above, the melting point is gained to be 725.08K for the heating process.The reason why this program gets two different melting points is because of overheating and overcooling.Before the system reaches equilibrium, the next time step of heating or cooling is carried

4. 2 .Figure 3 . 2 1 6
Figure 3.The pair coefficient equation under different temperature.From left to right, these figures (a)-(g) corresponding to 300 K, 400 K, 500 K, 600 K, 700 K, 800 K, and 900 K, respectively.Here are some discussion and analysis about the obtained pair correlation equation.In the low temperature situation, the pair correlation equation reaches the first peak at r=1.1  , which does make sense.The reason is that the Lennard-Jones potential has the minimum at   = 2 1 6  ≈ 1.12.As the temperature increases, the fluctuations in the end part of the image gradually decrease.(continue) This shows the change from solid to liquid state.With the increase of the temperature, the pair correlation value gradually becomes lower and lower.One of the reasons is the flattening of peaks.In order to integrate to obtain the same number, the highest value of the more widely distributed peak is smaller.Another reason for it is the transition to the gas phase.Due to the limitations of the two-dimensional simulation, the forces between the particles to maintain the solid and liquid state are not strong (due to the lack of vertical forces).Thus, as the temperature rises, the particles tend to disperse and some become gas phase.