Young’s Modulus Calculation of Some Metals Using Molecular Dynamics Method Based on the Morse Potential

It has been investigated computationally Young's modulus of some metals: nickel, copper, silver, gold, and aluminum. The offset method can graphically determine Young's modulus property by determining the elastic region based on the straight line intersection formed at a 0.2% strain against the stress-strain curve. In this simulation work, Young’s modulus calculation was performed by using the LAMMPS molecular dynamics software. The interatomic potential used to represent the interactions among atoms of materials in this simulation is the Morse potential. The metals under-investigated in this work are nickel, copper, silver, gold, and aluminum, and we got the results are 209.2 GPa, 110.8 GPa, 83.8 GPa, 79.2 GPa, and 70.3 GPa, respectively. The Young's modulus of the materials was also computed as temperature variations from 300K to the melting point to determine the effect of temperature on Young's modulus, and it is tensile strength. From our work we can found that the higher the temperature, the lower Young's modulus value. In addition, it can be seen that nickel metal has good temperature resistance. This is evidenced by the change in the nickel-metal phase near its melting point.


Introduction
The development of new/novel materials utilizing computational and simulation methods is one of the creative ways to optimize and streamline material research before direct experimental synthesis. Material research by reviewing the microscopic structure of materials can be used to predict materials macroscopically. By modeling and simulating materials in micro sizes, the macro-physical properties of materials can be estimated using Young's modulus of copper, silver, gold, aluminum, and nickel metals.
The elasticity property is the tendency of solid material to return to its original shape before permanently deformed [1]. Deformation in a solid material occurs when a force is applied to it. A measure of the degree of elasticity of a material indicates the material's resistance to elastic deformation due to external forces known as Young's modulus [2].
Investigation of Young's modulus of metallic materials has ever also been previously carried out by Camprubi [3] for nano-sized copper wire material. Ferguson [4] also studied the Young modulus of pure aluminum metal-based on EAM potential. Even many studies developed the EAM potential for metals. In our work, the calculation of Young's modulus is carried out using the molecular dynamics method based on metals' Morse potential. The fitting process is carried out to get proper morse potential parameters based on available experimental data (at 300K) before further evaluation to determine Young's modulus's temperature dependence. The use of Morse potential will simplify our simulation as a preliminary study of metals. We use the LAMMPS molecular dynamics code for simulating the phenomena (lammps.sandia.gov).

Young's Modulus
The Young's modulus of the material is closely related to material characteristics based on the degree of stiffness when external forces influence it. Materials with large Young's modulus are relatively non-elongated so that tremendous stress is required to produce deformation [5].
Young's modulus of material can be determined from a macroscopic perspective by knowing the ratio of stress to strain, known as Hooke's law.
where: E = young's modulus (GPa) σ = stress (GPa) ɛ = strain F = tensile force on metal (eV/Å) A = surface area of metal (Å 2 ) = length of metal before being subjected to a load (Å) = length of metal after being subjected to a load (Å) = change in length of metal before and after being subjected to a load (Å) Calculation of Young's modulus of pure metal is determined using the offset method by calculating the stress's slope to the elastic region's strain curve, namely the strain 0.002 [9]. The limit of the elasticity of the material is known as the yield stress (Yield Experimentally, the study of the effect of temperature on Young's modulus can be carried out using a tuning fork made of pure aluminum metal as the material [7]. Based on this research, it can be seen that the increase in material temperature affects Young's modulus value. The higher the temperature of the material, the lower Young's modulus value. The load applied to the material continuously will cause the stress and strain relationship to be non-linear. This area is called the plastic area, where the material cannot return to its original state when the applied load is removed. After passing the material's yield point, the stress will continue to increase for a specific limit called maximum stress. In this condition, the material can accept. The maximum stress is usually referred to as the tensile strength or UTS (Ultimate Tensile Strength) [8]. After passing through the tensile strength, the material's ability to accept loads will decrease until the material breaks (fracture).

Metals
Each metal has different characteristics so that several tests can be carried out, for example, tensile testing, impact tests, and hardness tests to determine the properties of the metal, especially its mechanical properties [9]. The following are some explanations of the metal materials used.

a) Copper
Copper is a metal that has an FCC crystal structure with a lattice constant of 3.6062 Å [10] and a melting point of 1358 K [11], which is often applied as nanowires. Copper has Young's modulus of 110 GPa [5] and is often used in uranium alloy compositions as an aircraft component.

b) Silver
Silver is a metal with a relative mass of 107.87 grams/mol, composed of an FCC crystal structure with a lattice constant of 4.0729 Å [10]. Based on its mechanical properties, silver has a lower Young's modulus value than copper and nickel, namely 82.5 GPa at 293 K [4].

c) Gold
Gold is an ideal metallic material for wires and electrodes at the nanoscale [12]. At the size of nanoparticles, gold and silver metals have been widely used as sensors, catalysis, biochemistry, optics, and electronics [13]. Gold is composed of an FCC crystal structure with a lattice constant of 4.0588 Å [14] and has a melting point of 1338 K [11].

d) Aluminum
Aluminum is a light metal with an atomic mass of 26.981 grams/mol and has resistance to corrosion [15]. The metal comprises an FCC crystal structure with a lattice constant of 4.0305 Å [10], Young & Freedman [5] stated in their book that aluminum metal has Young's modulus value of 70 GPa. Aluminum will undergo a phase change to become liquid at a temperature of 933.5 K [11]. Nickel is a metal that has an FCC crystal structure with a lattice constant of 3.5214 Å [10] and is often used in industry as a building block for steel. One of the advantages of nickel is that it has corrosion resistance properties. However, in a pure state, it is soft so that when combined with several other metals such as chrome and iron, it will be the right choice to form hardened corrosion-resistant steels. The metal has a relative mass of 58.71 grams/mol [16]. Nickel is very easy to combine with other metal elements, so its use is essential as a constituent of metal alloys for corrosion and heat resistance considering its high melting point is 1728 K [16].

Molecular Dynamics Methods
Molecular dynamics simulations predict the physical quantities of interactions at the atomic or molecular level based on the designed material model and based on the input simulation data given [17]. The research method used in this research is the classical molecular dynamics method, which can solve or find solutions to Newton's equations of motion using a potential function by the results of the atomic trajectory solution, which in this study is based on Morse's potential. Based on this, there is a relationship between the macroscopic properties of matter and the interactions at the atomic or molecular level.

Interactions Among Atoms
At the microscopic scale, each atom that is close to each other will interact with each other due to the atoms' forces. So it requires an interaction force between atoms, which is equivalent to the interaction potential used. The selection of interaction potentials greatly determines the results of the simulation performed. Rapaport [18] stated two main principles in the interaction between atoms: attractive forces and repulsive forces between atoms. First, if the atomic pairs are separated at a close distance, the resultant interaction forces will repel each other. Second, if the atomic pairs are separated at a great distance, the resultant interaction forces will attract each other. However, at a certain distance, the interaction styles that occur will cancel each other out so that the resultant interaction force will be zero. There is a dependence between the attractive interactions and the repulsive interactions between atoms on the potential energy, according to Figure 1   The interaction potential between atoms consists of a position function representing the atoms' potential energy when located in a particular configuration. This can be described using the relative position of an atom to another atom, which can be explained by the equation: The force that makes the atoms in an atomic system interact unbound when the atoms are close together or neighboring can survive in composing the crystal is the Van der Waals's force [19].

Morse Potential
Apart from the Lennard Jones potential and the EAM (Embedded atomic Method), another potential is often used in molecular dynamics simulations, namely the Morse potential. This potential has been widely used to explain the various properties of crystals. The characteristics of these crystals include cohesive energy, the lattice constant, compressibility, and elasticity constant. The usual Morse potential function is written in the following equation: where, U(r ij ) = interaction energy between atoms i and j separated at a distance (eV) D = potential well depth and is the energy parameter (eV) α = potential width controlling parameter (Å-1) r = distance between atoms (Å) r 0 =equilibrium distance between atoms i and j at minimum potential (Å) The Morse potential parameters of some pure cubic metals are composed of the FCC crystal structure, according to

Equation of Motion
In classical molecular dynamics, each atom and molecule is treated as a point mass that satisfies Newton's motion equations. The simulation is based on classical mechanics for a system of N atoms interacting through a potential function. This equation of motion is known as Newton's second law, namely: is the net force in newtons on atom i, determined from the negative gradient of potential energy based on atomic coordinates.
N is the number of atoms in the simulation [3]. Information about the correlation between the microscopic and macroscopic states of the material can be explained using an ensemble. The system under review is simulated in an isobaric-isothermal ensemble. In this ensemble, the number of particles, pressure, and temperature as macroscopic properties of the system can be maintained at a constant value while the volume of the system changes. Volume changes occur as a result of the load on the system.

Results and Discussion
The selection of the potential function and its parameters is essential for obtaining simulation result data suitable for the experiment. To get results consistent with experiments on macrosized materials, the system is conditioned to a pressure of 1 atm. The results obtained in this study are:

Correction of Morse Potential Parameters (D, α and r 0 )
Correction of the Morse potential parameter in the pure metal system was carried out to determine Young's modulus value at a temperature of 300K. The Young modulus value data was obtained, which was by the experimental data. Furthermore, these three parameters are used to determine Young's modulus value with temperature variations leading to its melting point. These parameters are obtained by running the simulation based on the input given. Based on the correction of the Morse potential parameter that has been done. It is obtained that the parameter values D, α, and r 0 are appropriate for the case under study. Parameter D shows the potential well depth expressed in energy units (eV). The parameter α is a parameter controlling the potential width expressed in (Å -1 ). Meanwhile, r 0 shows the equilibrium distance between atoms i and j when the minimum potential is expressed in units (Å). Parameters D, α, and r 0 on the calculated Young's modulus, along with the percentage of the discrepancy, are by

Determination of Young's Modulus Value of Pure Metal at 300 K
The calculation of Young's modulus values for the five materials is carried out based on the deformation system obtained. The deformation system simulates the material's mechanical properties by using a tensile test on the x-coordinate vector so that the system changes shape and size. The deformation system consists of stress (σ) at coordinates (X, Y, Z) and strain (ε).
The following is the stress-strain curve generated from the five materials deformation system at a temperature of 300K. According to Figures 2,3,4,5 and 6.  Figure 2. Pure copper stress-strain curve at a temperature of 300K Based on Figure 2, it is shown that the limit of the elasticity area occurs at a stress of 4.9 GPa with a strain of 0.05. Also, it can be seen that the tensile strength of the copper metal system is 11.4 GPa. which occurs at a strain of 0.13. Figure 2 shows that the system experiences a sharp stress drop when it passes through the tensile strength region of up to 4.5 GPa, which occurs in the strain area between 0.14 to 0.16. This shows that in this area, the system has almost lost resistance to a given load. At a strain of 0.17 to 0.3, it can be seen that the system has a small load resistance, which is indicated by a stress value range of 2.5 to 3 GPa.  Figure 3, it is shown that the limit of the elasticity area occurs at a stress of 3.1 GPa with a strain of 0.04. In addition, it can be seen that the tensile strength of the copper metal system is 6.9 GPa, which occurs at a strain of 0.1. Figure 3 shows that the system has experienced a sharp stress drop from the tensile strength area to 0.7 GPa, which occurs in the strain area between 0.11 to 0.12. However, at a strain of 0.13 to 0.3, it can be seen that the system goes into a state where the resistance to loading is zero. This indicates that the system has fractured.  Figure 4. it is shown that the limit of the elasticity area occurs at a stress of 3.4 GPa with a strain of 0.04. In addition, it can be seen that the tensile strength of the copper metal system is 7.1 GPa, which occurs at a strain of 0.1. Figure 4.8 shows that the system has experienced a sharp stress drop from the tensile strength area to 0.8 GPa, which occurs in the strain area between 0.12 to 0.14. However, at a strain of 0.15 to 0.3, it can be seen that the system goes into a state where the resistance to loading is zero. As with silver material, it indicates that the system has fractured. Figure 5. The stress-strain curve of pure aluminum at a temperature of 300 K Based on Figure 5, it is shown that the limit of the elasticity area occurs at a stress of 2.6 GPa with a strain of 0.04. In addition, it can be seen that the tensile strength of the copper metal system is 5.3 GPa, which occurs at a strain of 0.09. In Figure 4.9. it can be seen that the system has experienced a sharp stress drop from the tensile strength region to 0.3 GPa, which occurs in the strain area between 0.1 to 0.11. However, at a strain of 0.11 to 0.3, it can be seen that the system goes into a state where the resistance to loading is zero. This shows that after a sharp drop in stress, the system fractures.  Figure 6. The stress-strain curve of pure aluminum at a temperature of 300K Based on Figure 6, it is shown that the boundary of the elasticity area occurs at a stress of 8.3 GPa with a strain of 0.04. Also, it can be seen that the tensile strength of the copper metal system is 18.9 GPa, which occurs in a strain of 0.12. Figure 6 shows that the system experiences a sharp stress drop from its tensile strength area to 5.3 GPa, which occurs in the strain area between 0.12 to 0.14. However, at a strain of 0.15 to 0.3. it can be seen that the system has a small load resistance, with a stress value range of 4 to 5 GPa. Based on Figures  2,3,4,5 and 6, it can be concluded that the increase in the tensile strength of the material is proportional to the increase in the Yield stress of each material. Also, the copper and nickel metal systems do not fracture as in silver, gold, and aluminum at the end of the simulation.
Referring to the Morse potential parameter in Table 3, the Young modulus value for each pure metal material at a temperature of 300K is obtained along with the value of the discrepancy, according to Table 4 Based on the simulation results, the system's stress has a GPa scale, considering the material used is nano-sized metal, which has a regular arrangement of atoms. Whereas for materials such as macro-sized metals and composites, the stress produced has an MPa scale, which is undoubtedly far from the simulation results. This is because the order of the atoms is not as good as the system in the simulation.
Based on Table 4, it is known that the Young gold modulus has a much better level of accuracy than other metals. The smaller discrepancy value of gold indicates this compared to other metals, namely 0.2%. The lowest level of accuracy occurs in silver metal, where the discrepancy value is 1.6%. However, each of the Young modulus values above can be said to be following experimental data considering that the discrepancy value obtained is less than 2%. The accuracy of Young's modulus values for these five materials is highly dependent on the suitability of the Morse potential parameters used.

Young's Modulus at Various Temperature Variations
In addition to the temperature of 300K, the determination of the Young modulus value of pure metal using molecular dynamics methods is also carried out based on temperature variations. The temperature used starts from 300K to close to the melting temperature. Temperature variation is carried out in one direction, from the smallest to the largest value. This was done to determine the resistance of the material to loading during extreme temperatures. Knowing Young's modulus value shows how much stress the material can accept up to the elasticity limit of the material. For that, a stress-strain curve is needed for each temperature variation. Figure 7 shows the stress-strain curve for pure copper metal with temperature variations.  Figure 7, it can be seen that the increase in temperature applied to copper metal causes a decrease in the value of its tensile strength. The more significant the decrease in tensile strength, the easier the material is deformed. The graph shows that copper metal undergoes a phase change to become liquid at a temperature of 1273K before reaching its melting point at 1358K. This is evidenced by the absence of a stress-strain curve pattern as at the previous temperature.
At a temperature of 1273K, the calculation of Young's modulus value cannot be carried out considering that the stress-strain curve is linear along the x-axis with a stress value of 0 GPa along with the strain. At an initial temperature, such as at a temperature of 300 K, it can be seen that there is a sharp drop in stress when the load continuously exceeds its tensile strength. This shows that the metal is almost losing its ability to defend itself from deformation due to the load received until, in the end, the metal experiences continuous elongation with relatively constant stress. However, when the metal temperature has approached its melting temperature, such as at 1223K, the stress gradually decreases. This indicates that at this temperature, copper metal has a low resistance to deformation. In addition to copper metal, the resistance of the material to a given load with temperature variations in silver metal is shown in Figure 8. The increase in temperature applied to silver metal causes a decrease in the value of its tensile strength, as in Figure 8. The graph shows that silver metal undergoes a phase change to become liquid at a temperature of 873K, long before reaching its melting point, which is at 1235K. This is evidenced by the absence of a stress-strain curve pattern as at the previous temperature. At a temperature of 873K, the calculation of Young's modulus value cannot be carried out considering that the stress-strain curve is linear along the x-axis with a stress value of 0 GPa along with the strain. The increase in temperature applied to the gold metal also causes a decrease in the value of its tensile strength, as shown in Figure 9. The graph shows that the gold metal undergoes a phase change to become liquid at a temperature of 873K, long before it reaches its melting point at 1338K. This is evidenced by the absence of a stress-strain curve pattern like at the previous temperature. At a temperature of 873K, the calculation of Young's modulus value cannot be carried out considering that the stress-strain curve is linear along the x-axis with a stress value of 0 GPa along with the strain. At a temperature of 823K, the voltage drops slowly. This indicates that at this temperature, gold metal has low resistance to deformation. Figure 10. Stress-strain curves on pure aluminum with temperature variations Based on Figure 10, it can be seen that the increase in temperature applied to aluminum metal also causes a decrease in the value of its tensile strength. The graph shows that aluminum metal undergoes a phase change to become liquid at a temperature of 723K before reaching its melting point (933.5K). At a temperature of 723 K, the calculation of Young's modulus value cannot be carried out considering that the stress-strain curve is linear along the x-axis with a stress value of 0 GPa along with the strain. Figure 11. The stress-strain curve in pure nickel with temperature variations Based on Figure 11, the metal's resistance to deformation at low temperatures, which is between 300K and 573K, is still quite good. This is evidenced by the metal's ability to maintain the strain that occurs at the point of its tensile strength, which is still quite good. On the other hand, at temperatures of 1423K to 1673K, the metal is seen to experience a relatively more significant decrease in tensile strength than at the previous temperature. Figure 11 shows that nickel-metal undergoes a phase change at a temperature of 1723K. Young's modulus can be determined when the material is still in a solid phase. Solid material has a tremendous inter-atomic binding energy compared to the liquid and gas phases. This is related to the orderly arrangement of atoms so that the bonds between atoms that occur are also stronger. In the liquid phase, the arrangement between atoms is more irregular. That is why it is difficult to determine the Young modulus when it approaches its melting point. The relationship between the Young modulus values of pure metals and their temperature variations is shown in Figure 12.  Figure 12, it can be seen that each pure metal has decreased with the increase in temperature. Microscopically, the liquid phase of the material has a greater distance between atoms than the solid phase. This causes the binding force between atoms in the liquid phase to be smaller than in the solid phase.

Visualization of the Effect of Temperature
The temperature variations are given to each pure metal material affect the structure of the material. This is indicated by the percentage change in crystal structure at any given temperature variation. Based on the simulation results that have been carried out. When given temperature variations under the effect of a given loading, the behavior of the material can be seen in the visualization display using the OVITO software. Based on this visualization, the system state of each material can be seen. The state of the system starts from a temperature of 300K, as in Figure 13, to a state where the material undergoes a phase change to become liquid, as in Figure 14.  Figure 13. Visualization of nickel (a) before loading with (b) after loading at a temperature of 300 K Figure 13 shows that at a temperature of 300 K, nickel metal is composed of a uniform crystal structure with a percentage of the FCC crystal structure of 100%. In contrast to the case when the crystal structure of nickel metal is no longer uniform after loading. This was evidenced by the CNA display, which showed the percentage of HCP's crystal structure was 14.3%, the FCC was 55.8%, and 30% for the crystal structure that was undefined. Green shows the FCC's crystal structure, red shows the crystal structure of the HCP, and gray shows the crystal structure that is undefined. At a temperature of 300 K, it can be seen that the atoms making up the nickel-metal material are still in the solid phase. This is indicated by the absence of a cavity or free space in the simulation box, where the space indicates that the metal system has undergone a phase change, although not entirely. In addition to the temperature of 300 K, the nickel-metal state at 1723 K can be observed based on Figure 14.
(a) (b) Figure 14. Visualization of nickel (a) before loading with (b) after loading at a temperature of 1723 K Figure 14 shows that at 1723 K, nickel metal is composed of a small portion of the FCC crystal structure, namely 1% and 99% of the crystal structure is undefined. In contrast to the case when after loading, the nickel-metal's crystal structure becomes uniform, namely 100% composed of an undefined crystal structure. At a temperature of 1723 K, it can be seen that the atoms making up the nickel-metal material have undergone a phase change from solid to liquid. This is indicated by a cavity or free space on the edge of the simulation box. Microscopically, metals in the solid phase are composed of atoms with binding energy greater than those in the liquid phase. This shows that the atoms' binding energy making up the nickel-metal at a temperature of 1723 K has a relatively small value, causing the distance between the atoms and the volume of the system to be much larger than in the solid phase. The profile of the effect of temperature variations on copper, silver, gold, and aluminum materials is the same as for nickel metal. The other four materials' visualization is shown in the visualization of the effect of temperature on the pure metal material system.

Conclusions
The research's success can be reviewed based on Young's modulus value of the five materials according to the experimental results, namely with a discrepancy value of less than 2%. Also, the accuracy of selecting the interatomic potential function based on the appropriate parameters is an essential factor in obtaining precise and accurate results, which in this study uses the Morse potential, which is suitable for metallic materials in the solid phase. The simulation of the calculation of Young's modulus value is carried out using the molecular dynamics method. This method can be developed in more complex cases, one of which is the stiffness of alloys and composites.