Improved energy minimization of iron–carbon systems: on the influence of positioning interstitial atoms

Nina Gunkelmann; Maximilian Merkert

doi:10.1088/1361-651X/ab6bb6

1. Introduction

The energy minimization in molecular dynamics simulations (MD) is a highly complex optimization problem with a large number of optima comprising of differing qualities. Depending on the potential energy surface, numerous local energy minima exist which can be separated by high energy barriers. In general, one is interested in the global minimum to find the mechanically stable ground state of the systems. Due to the small timescales of MD, simulations of such systems can suffer from inadequate sampling because they are trapped in local minima.

There exist many different optimization schemes used in MD simulations including steepest descent, conjugate gradient (CG) and Newton–Raphson methods [1–3]. Bitzek et al [4] introduced a simple optimization algorithm which is significantly faster than standard implementations of the CG method. However, by using these optimization methods only local minima can be found.

One common method to improve the sampling of atomistic simulations is simulated annealing (SA) [5]. In SA simulations, a system is given an initially high temperature, allowing it to overcome high energy barriers, and then gradually cooled to the desired temperature. The SA technique finds a global optimum of any finite problem with a probability converging to 1 as the time taken to cool the system tends to infinity [6]. However, the time scale in MD is limited to several 100 ns.

Also other Monte Carlo methods are used for energy optimization—either for optimizing the configuration of C-atoms [7] or for the whole alloy [8], obviating the need for a following local descent phase. In both cases, from a randomly constructed starting configuration, a Markov process with a Boltzmann factor may be used to obtain a thermodynamic state which obeys the ground-state configuration at sufficiently small temperature [9]. Grand-canonical ensemble Monte Carlo schemes have proven very useful for dealing with density fluctuations. With the number of atoms fixed, they, however, share similar strengths and weaknesses with SA. In particular, there is no guarantee for Monte Carlo algorithms to reach a global minimum (not even a local minimum for many versions) in a finite amount of time.

In iron–carbon martensite and ferrite the C-atoms are not randomly distributed. It was found by Moessbauer spectroscopy that neighboring carbon atoms in octahedral sites repel [10]. This is coincident with the fact that carbon atoms obey a high mobility in ferrite. In reality, a random distribution of atoms is only expected at very high temperatures where the entropy overwhelms any ordering of the atoms.

However, many researchers randomly distribute carbon atoms in octahedral sites of the bcc lattice followed by a conventional optimization step using a CG method [11–13]. On the one hand, this procedure runs the risk that the configuration ends up in a secondary local minimum. On the other hand, the statistical distribution of interstitial atoms has to be taken into account resulting in increased computational cost. Especially for large-scale MD simulations it is often not possible to include statistical averaging of the results in a reasonable computation time.

In this work, we compare three different optimization schemes for the energy minimization of bcc or bct Fe–C martensite with high carbon concentrations using periodic boundary conditions. Martensitic steels are of interest in many applications because of their high strength and toughness. The mechanical properties of the metastable martensite systems usually formed by rapid cooling of the fcc structure are evaluated after minimization using the following.

CG with randomly distributed interstitial atoms in the initial state.
CG with interstitial atoms initially in globally optimal interstitial sites with respect to a long-range repulsive potential for the carbon–carbon interaction in Fe–C [14], computed by exact methods of discrete optimization.
SA with randomly distributed interstitial atoms in the initial state.

Although the global minimum is only achieved for T = 0 K, low-energy local minima will also be preferred for finite temperature. After SA to low temperature for sufficiently long durations, the system will achieve a local energy minimum which is close to the global minimum. Global optimization is also important to obtain detailed information about the energy landscape of molecular configurations and was frequently used for the energy minimization of complex molecular systems [15–17].

We present the simulation and optimization methods in section 2. The results of our simulations using different initial Fe–C configurations are compared in section 3, where we discuss energy differences in section 3.1 while section 3.2 addresses the elastic properties of the system. Finally, in section 3.3, we evaluate dynamic properties of the system by calculating the yield pressure.

2. Methods

Metals are described using the embedded-atom method potential. The total energy of the alloy can be subdivided into contributions from the different types of particle pairs, namely iron–iron, iron–carbon and carbon–carbon interactions. The energy contribution of the latter is given by

$\begin{eqnarray}&&E=\displaystyle \sum _{i=1}^{N}\left({F}_{{\rm{C}}}\left({\rho }_{i}\right)+\frac{1}{2}\displaystyle \sum _{j\ne i}{\phi }_{\mathrm{CC}}(\left|{\vec{r}}_{i}-{\vec{r}}_{j}\right|)\right),\end{eqnarray} \tag{ 1 }$

where the embedding energy F_C describes the bonding of atom i. It depends on the electron density at a given position:

$\begin{eqnarray}&&{\rho }_{i}=\displaystyle \sum _{j\ne i}{\rho }_{{\rm{C}}{\rho }_{{\rm{C}}}}\left(\left|{\overrightarrow{r}}_{i}-{\overrightarrow{r}}_{j}\right|\right).\end{eqnarray} \tag{ 2 }$

The function ϕ_CC(r) describes a pair-wise interaction. Contributions of iron–iron and iron–carbon interactions are described by similar terms, but different functions F_Fe, ρ_Fe, ϕ_FeC and ϕ_FeFe.

The problem of finding energy-optimal positions for all atoms is a very difficult optimization problem that we cannot expect to solve to global optimality for nontrivial sizes. Therefore, one usually employs a local search method to find a local optimum, such as the CG method [11–13, 18–20]. However, carbon atoms are very unlikely to move to other interstitial sites during a run of any local descent method. Instead, the structure of the problem suggests that there will be at least one local optimum for each assignment of carbon atoms to interstitial sites. Consequently, finding the best interstitial sites as a starting configuration for a CG algorithm is of crucial importance if we want to get close to a global energy optimum.

In order to measure the quality of an interstitial configuration, we assume that all lattice positions are occupied by Fe atoms. In this case, the contribution of iron–iron interactions is fixed. Furthermore, if carbon atoms occupy octahedral sites, the energy contribution of iron–carbon interactions is also constant due to periodic boundary conditions. Therefore, the only variable term to be optimized is the contribution by carbon–carbon interactions given in (1). Note that optimizing this value will in general not lead to a globally energy-optimal configuration due to the distortion of the Fe lattice during the CG method. Each C interstitial leads to a local lattice deformation that will affect the surrounding Fe–C bonds but also the energetics of the Fe lattice during the subsequent CG optimization considering all interactions between iron and carbon atoms. However, it gives a very good indication of the quality of a local optimum associated to the configuration. The reason is that the distortion after applying a CG method is small. For globally optimal initial C configurations with respect to ϕ_CC(r), the bcc lattice is less distorted. For this reason, the difference in interaction energy compared to all Fe atoms remaining in lattice positions (before CG minimization) remains small resulting in a small total energy state.

To find optimal interstitial positions with respect to (1) we use the long-range purely repulsive potential between carbon atoms in Fe–C by Belashchenko et al [14]. The authors used this potential to study high-temperature liquid Fe–C. The cutoff of this potential is at 8.6 Å. In Fe–C melt alloys the mixing enthalpy is negative [21]. As a consequence, there is no tendency for C-atoms to be clustered. Therefore, we expect that C-atoms do not appear at short distances after solidification but their distances are maximal. Note that any repulsive C–C potential with an appropriate long-range interaction could be used to optimize the interstitial positions. For the calculation of the elastic constants in solid Fe–C we use the potential by Becquart et al [18] which has been shown to fairly reproduce elastic properties in Fe–C [19, 20]. The cutoff of this potential amounts to 4.8 Å.

Carbon occupies octahedral sites in the bcc Fe lattice in the edge centers and in the face centers of the lattice ${\vec{r}}_{i}=l(a,b,c+1/2)$ , ${\vec{r}}_{i}=l(d+1/2,e+1/2,f)$ and equivalent sites ${\vec{r}}_{i}=l(c+1/2,a,b)$ , ${\vec{r}}_{i}=l(f,d+1/2,e+1/2)$ , ${\vec{r}}_{i}=l(a,c+1/2,b)$ , ${\vec{r}}_{i}\,=l(d+1/2,f,e+1/2)$ where a, b, f ∈ [0,1, ..., N_max] and $c,d,e\in [0,1,\ldots ,{N}_{\max }-1]$ with the lattice constant l (see figure 1 showing octahedral sites in bcc iron). The number of unit cells is N_max and the total number of C-atoms is fixed. We use iron–carbon crystals containing 128 and 250 Fe atoms and vary the number of C-atoms from 2 to 8.

**Figure 1.** Schematic view of octahedral sites in bcc Fe. Iron atoms are shown in red and octahedral sites are marked in black.
Download figure:
Standard image High-resolution image

We apply periodic boundary conditions affecting the distance between C-atoms occupying octahedral sites

$\begin{eqnarray}&&\left|{\vec{r}}_{i}-{\vec{r}}_{j}\right|=\left|\left[\begin{array}{c}\min \left(\left|{r}_{{\text{}}{ix}}-{r}_{{\text{}}{jx}}\right|,\left(-\left|{r}_{{\text{}}{ix}}-{r}_{{\text{}}{jx}}\right|+{r}_{\max }\right)\right)\\ \min \left(\left|{r}_{{\text{}}{iy}}-{r}_{{\text{}}{jy}}\right|,\left(-\left|{r}_{{\text{}}{iy}}-{r}_{{\text{}}{jy}}\right|+{r}_{\max }\right)\right)\\ \min \left(\left|{r}_{{\text{}}{iz}}-{r}_{{\text{}}{jz}}\right|,\left(-\left|{r}_{{\text{}}{iz}}-{r}_{{\text{}}{jz}}\right|+{r}_{\max }\right)\right)\end{array}\right]\right|.\end{eqnarray} \tag{ 3 }$

In this work, we focus on martensite Fe–C which is of special interest because of its high hardness [22]. It is formed by rapid cooling of the austenite fcc structure at high strain rates and is supersaturated with carbon. Note that martensite in high carbon steels (0.3–1.5 wt% steel) is metastable at room temperature [23].

In the iron–carbon system carbon atoms are usually situated at the octahedral sites distributed along one axis to assure tetragonality of the system. The tetragonal distortion of the Fe–C lattice can be explained by the Bain transformation from fcc to bcc: before the transition the carbon atoms are located in the octahedral sites of the fcc phase: in the middle edges and in the center of the cell. After transformation, the C-atoms are distributed in octahedral sites of the bcc phase, preferably along the c-axis. Since the carbon atoms are larger than the octahedral sites, the crystal will be tetragonally distorted (see figure 2).

**Figure 2.** Schematic view of the Bain transformation for the fcc-bct phase transition. Blue atoms belong to fcc and red circles denote the newly formed bcc unit cell. After the lattice rearrangement the C-atoms do not distribute randomly but preferably along the c-axis.
Download figure:
Standard image High-resolution image

In the following, we apply different simulation set-ups to distribute a given number of carbon atoms in the lattice.

(i)
Random cubic scheme: all octahedral positions of the crystal are equally probable to be occupied.
(ii)
Random tetragonal scheme: C-atoms are randomly set on interstitial octahedral positions of the same symmetry. This scheme has been previously applied by Janßen et al [20]. The justification of this scheme is the tetragonality of the martensiteFe–C lattice explained by the Bain transformation (see figure 2 above).
(iii)
Optimal cubic scheme: carbon occupies octahedral interstitial sites of all symmetries such that the total energy, equation (2) is minimized. For a detailed description see section 2.1 below.
(iv)
Optimal tetragonal scheme: carbon occupies octahedral interstitial sites such that the total energy, equation (2) is minimized. C-atoms are only set on interstitial octahedral positions of the same symmetry to ensure the tetragonal distortion explained by the Bain transformation in figure 2.
(v)
Single-C-atom scheme: a single C-atom is set in an octahedral site in the middle of the crystal. By varying the size of the ferrite matrix different C contents can be realized. Using periodic boundary conditions this configuration leads to a regular distribution of C-atoms. The scheme has been previously applied in the literature [18–20].

The simulations were carried out by the open-source MD code LAMMPS [24]. After creating these set-ups the potential energy is minimized using the conjugated-gradient method. We used a stopping energy tolerance of 10⁻²⁵ in LAMMPS (with the meaning that the energy change between successive iterations divided by the energy magnitude is less than or equal to the tolerance) and a force convergence value of 10⁻²⁵ eV Å⁻¹.

In addition, we compare the results with SA simulations where interstitial sites of any symmetry are chosen. The temperature is linearly decreased from T = 700 K to T < 1 K during 50 ps via velocity scaling. During SA we use a barostat to minimize the pressure. Note that we do not simulate the austenite-martensite transformation during annealing of fcc Fe–C at higher temperature because the potential used is not able to capture the temperature-induced phase transformation. Instead, we start with supersaturated bct iron–carbon to study SA of martensite.

2.1. Statically optimal interstitial schemes

Assuming Fe and C-atoms occupy lattice positions and a subset of the octahedral sites, respectively, the total energy of the system is completely determined by the octahedral sites chosen for the C-atoms. Hence, we aim to find an energy-optimal configuration among a finite—though very large—number of possible configurations. This is a discrete optimization problem that we can can model as a binary quadratic program: for each octahedral site i = 1,...,n we use a vector x of binary variables to express whether site i is occupied (x_i = 1) or empty (x_i = 0). This leads to the formulation

$\begin{eqnarray}&&\begin{array}{l}\mathop{\mathrm{minimize}}\limits_{x}\quad {x}^{T}{Ex}\\ \mathrm{subject}\ \mathrm{to}\quad \displaystyle \sum _{i=1}^{n}\ {x}_{i}=k\\ \qquad \qquad \qquad x\in \{0,1{\}}^{n},\end{array}\end{eqnarray} \tag{ QP }$

where E_ij denotes the energy contribution between two carbon atoms placed at positions i and j given by (1). These values can be computed beforehand and are therefore constants for the optimization problem. The single side constraint ensures that the correct number k of C-atoms representing the desired concentration will be placed; note that using an '≥'-inequality instead of an equation would also be sufficient if the potential is long-range and strictly repulsive, as it is in our case.

This model is employed for (iii) and (iv) from the list of simulation setups above, where for (iv) consequently only a subset of all octahedral sites is considered. We use the general-purpose mixed-integer programming solver Gurobi [25] (version 7.5 with standard parameter settings) to solve the resulting problems to global optimality. Note that formulation (QP) is highly symmetric. We could exploit this by several symmetry breaking techniques—for instance by setting x₁ = 1 to mention one very simple option—though this is not necessary as Gurobi 7.5 has automatic symmetry detection and handling. Furthemore, (QP) could be transformed to a mixed-integer linear program using the standard McCormick linearization [26], which the solver is also capable of doing by itself. The solver can also be used to find a given number of best solutions.

Note that an optimal solution of (QP) is only globally energy-optimal for the alloy assuming that all atoms stay in place. This will of course not be true in general, but the distortion after applying a CG method is assumed to be small enough such that it gives a very good indication of the overall quality of the chosen interstitial sites, see section 3.1.

3. Material properties

3.1. Energies

Figure 3 shows the distribution of 4 carbon atoms in a cubic box with 84 possible octahedral interstitial positions using periodic boundary conditions. This corresponds to an iron crystal with 128 iron atoms and a carbon content of 3.0%. The atoms appear to be evenly distributed maximizing the distance between carbon atoms.

We display in figure 4 boxplots of the energy difference to the minimal value after minimization with the CG method of 50 random carbon positions ordered along the z axis together with the optimized solution for the Becquart potential and the purely repulsive long-range potential. Boxplots show quartiles and the band inside the box represents the median of the energy values. Lines extending from the boxes indicate the standard deviation. Outliers are plotted as individual points. The energies for the optimized solution after CG are shown as lines for both potentials. We observe that for the repulsive potential, even after CG, the optimized solution has lowest energy. For the Becquart potential the energy value for our solution is not optimal but close to the minimal observed value with a distance of 0.02 eV. Although the cutoff for the interaction between C-atoms is rather small for the Becquart potential (4.81 Å), the energy for an optimal configuration according to (QP) where the long-range potential has been used is close to the minimal value. A reason might be that for these configurations the bcc lattice is less distorted so that the interaction energy between Fe atoms remains small resulting in a small total energy state.

In figure 5 we display the energy difference to the minimal observed value after minimization with the CG method of 50 random carbon octahedral positions, together with the corresponding result from 50 SA simulations, and the optimized solution for the Becquart potential. Again, we observe that the optimized solution is very close to the smallest known value. As expected, the SA results are superior than the random CG solution. However, the average value is still higher than that of the optimized solution.

3.2. Elastic properties

We calculate the elastic constants by deforming the cell in one direction and measuring the change in the stress tensor using the Becquart potential. A CG algorithm is used to restore a locally optimal atom configuration after each deformation. The elastic constants are derived by determining the slope of the stress-strain curve. We applied both a compressive and a tensile strain in each direction and averaged the corresponding slopes. To avoid statistical fluctuations we use a very small temperature T < 1 K.

The elastic components are shown in figure 6 for the averaged results of 100 random carbon atom configurations ordered along the z axis, the optimized solution and the elastic constants for a single C-atom set in the middle of the iron crystal. Error bars in figure 6 denote the maximum deviation from the average values.

The elastic constant C₁₁ stiffens with increasing carbon content up to concentrations of c_C = 0.05 due to the contraction of the lattice in a direction. C₃₃ softens with increasing C content due to the expansion of the lattice in c direction. For small concentrations the optimized and random solutions agree well and the dependence on C content is approximately linear. For larger concentrations we observe that for the optimized solution, the stiffness is increased for both elastic constants because the expansion of the lattice is smaller in c direction if C-atoms are evenly distributed. The deviation is largest for high carbon content where also the fluctuations of the random solution is very large. Here, the values for C₁₁ cover a range of up to 50% of the mean value. At these large concentrations deviations from linearity occur which can be explained by the overlap of the elastic fields of the carbon atoms. The single-atom interstitial scheme deviates from both the optimized and the random scheme. Here, the distance between C-atoms is maximal due to the ordered lattice so that no C interference can occur [20]. This results in significantly increased stiffness of the elastic constants. The elastic constants C₁₂ and C₁₃ increase for small C content and decrease for higher concentrations. The optimized solution has an energy smaller than the random scheme. A reason might be that the increase in stiffness results in a decreased shear resistance and high embrittlement. The constants for random C interstitials can differ by over 5%, highlighting the problems of the random scheme. The constants C₄₄ and C₆₆ follow the trend of softening of the elastic constants with C content where only small deviations between random and optimized solution occur.

The random cubic scheme in figure 7 shows large differences to the other schemes where the optimized solution is plotted together with the result from SA and random configurations. Note that we performed SA only for the cubic scheme because C-atoms diffuse to all available octahedral sites during quenching. Because of the symmetry of the lattice, there exist only the elastic constants C₁₁, C₁₂ and C₄₄. Here, we observe only a very mild dependence of the elastic constants with C content. The elastic constant C₄₄ slightly decreases with C content while C₁₁ slightly increases for the optimized solution and remains constant for the SA and the random scheme. The error bars of the results for all schemes overlap. However, the error of the random scheme is significantly larger and the range of values is approximately 65%. Note that here we also have error bars for the optimized scheme because several minimal energy configurations with different elastic constants exist. We conclude that the random scheme should not be used because it leads to large errors.

To compare our results with the available experimental data [19, 27], we use the Voigt–Reuss–Hill approximation to convert single-crystal elastic constants into polycrystalline elastic moduli. The Voigt and Reuss bounds are upper and lower bounds on the effective elastic constants and the Hill approximation is the average of both values. The dependence of the bulk modulus B and the shear modulus G on the elastic tensor can be found in the paper by Watt and Pelsnick [28].

From these averages we get Young's modulus Y via:

$\begin{eqnarray}&&E=9{BG}/(3B+G).\end{eqnarray} \tag{ 4 }$

The Young's modulus for our data is plotted in figure 8. The experimental result by Ledbetter [19] of the softening of the elastic modulus with respect to C content is well reproduced by the simulation results. Here, the softening is weaker for the optimized solution which fits better to the data in [19]. The softening is stronger for the experiments by Speich et al [27] but the results show also increased fluctuations. This problem was avoided in the measurements by Ledbetter by measuring all elastic constants many times [19]. The results using the single interstitial scheme also reproduce the experimental data by Ledbetter in the range up to c = 0.02. However, the large variety in the experimental predictions makes it difficult to decide which solution is most realistic.

3.3. Yield pressure

To evaluate the changes in the dynamics of optimized and random configurations distributed along the z axis to maximize tetragonality, we apply compressive strain by the following procedure used by Wang et al [29] using the Becquart potential.

After each simulation time period of 1 ps, the system is homogeneously compressed along the z direction by 0.1%; then the system relaxes for the next 1 ps of simulation time under NVE conditions. The strain rate is 10⁹ s⁻¹. This method is similar to a temperature-accelerated MD simulation [30]. Note that times and strain rates are not comparable to experiments resulting in unrealistically high pressures. However, our results can reveal the qualitative behavior of the material in dependence of the C concentration.

To measure yield strains in Fe–C, we evaluate the position in the curve where the difference in pressures normalized by the yield pressure p_Fe for pure iron between two consecutive strain values is larger than 5%.

At this strain, plastic yielding occurs and the material transforms from the body-centered cubic (bcc) to the hexagonal close-packed (hcp) phase. Iron is known to show this structural phase transition at a pressure around 13 GPa [31]. The shock wave structure consists of three waves, corresponding to elastic, plastic and phase changed regions [32, 33]. The potential used in this study does not correctly predict the phase transition pressure in iron–carbon because the transformation barrier is too high. For pure iron, a number of potentials exist which provide reasonable phase transition behavior [34, 35]. The modified Ackland version [35] has been applied for polycrystalline [36, 37] and single crystalline iron [38]. For single crystals iron was found to yield via twinning and the phase transformation was shown to depend on the plasticity history [38]. The potential by Voter [34] typically gives a large fraction of fcc phase, instead of the experimentally expected hcp phase and does not show plasticity before the phase transformation [35]. For Fe–C the modified potential by Ackland was coupled to several Fe–C interaction potentials to feature the phase transition at the correct pressure [39]. However, it was found that the transformation pressure is significantly increased due to phonon drag.

Figure 9 shows snapshots of an iron crystal containing 2 million atoms at different strains. The samples were visualized using polyhedral template matching [40] within the software OVITO [41] using a RMSD cutoff of 0.2. The crystal transforms to the hcp phase containing stacking faults identified as fcc. We also observe twinning as noted by Amadou et al [38]. The twin planes in figures 9(b)–(d) are (111) planes in fcc and appear as stacking faults. Twinning relieves the stress during the transformation. Note that our Fe–C samples have smaller sizes but show qualitatively similar deformation behavior. The yield curve for this sample is shown in figure 10 for pure iron and Fe–C with concentration c_C = 0.045. For pure iron, we observe yielding at a strain of 0.230 and a pressure of p = 73 GPa corresponding to the onset of the phase transformation to the hcp structure. For Fe–C the yield point is decreased to 0.20 and the pressure is p = 58 GPa for the optimized C configuration. Note that we do not observe dislocation emission before the phase transition. The reason is that in single crystals we do not have dislocation sources. It would be interesting to study polycrystalline iron–carbon samples where grain boundaries act as sources which could be activated to nucleate dislocations under stress.

**Figure 10.** Pressure versus strain for pure Fe and Fe–C with concentration c_C = 0.045 during uniaxial compression.
Download figure:
Standard image High-resolution image

The pressure at the yield strain is equal to the yield pressure. This quantity describes the driving pressure for the onset of plastic activity and the phase transition in the Fe–C crystal. The iron–carbon crystals contain 128 Fe atoms and 3–8 C-atoms resulting in C concentrations from 0.02 to 0.059. We observe in figure 11 that the yield pressure is increased for the optimized solution in comparison to the results for random C positions. This is consistent with the increased Young's modulus in figure 8. The strength of the crystal is thus higher if the position of the C-atoms is optimized. Note that the yield strength is decreasing with increasing C concentration. The reason is that the elastic constants are softened with increasing carbon content which was also observed in experiments [19, 27]. In macroscopic real steels, the increase in strength with C concentration is due to a distortion of the lattice because of C-induced anisotropy. In addition, various Fe–C phases exist modifying the hardness of the crystals [42].

**Figure 11.** Relative yield pressure normalized by the yield pressure p_Fe for pure iron as a function of carbon content. The averaged results of 50 random C positionings ordered along the z axis is given in red and the optimized solution is colored in green.
Download figure:
Standard image High-resolution image

The yield strain in figure 12 decreases with increasing C content. It is smaller for the optimized solution. This is consistent with the increased hardness resulting in decreased ductility. Note that a measure for ductility is the Cauchy pressure C_p = C₁₂ − C₄₄ [43] which is also smaller for the optimized solution resulting in increased brittleness. Error bars in figures 11 and 12 denote the maximum error and are only shown for exemplary measurement points.

4. Conclusion

For finding energy-minimal states in Fe–C via CG methods the initial positioning for the interstitial atoms plays a crucial role. We tested fully randomized as well as single-C-atom schemes that are common in the literature; in addition, we proposed solving an auxiliary discrete optimization problem as described in section 2.1 for energy minimization. The resulting configurations have been compared with respect to system energy as well as elastic properties and yield pressure. Using randomly distributed carbon atoms is computationally inexpensive but cannot be recommended, as it may lead to energy-wise severely suboptimal solutions even after performing a reasonable amount of runs. Furthermore, it shows large errors in the elastic constants. The solution by SA is superior but still causes significant deviations. Using statically optimized solutions using a long-range repulsive potential leads to very-low-energy solutions also after subsequently applying a CG method and even if the potential used for the simulation is low-range. Furthermore, optimized solutions tend to an increased stiffness of the Fe–C system. Hence, using these configurations reveals increased yield strength of the material. We conclude that the choice of the positioning scheme has decisive influence on the mechanical properties iron-based alloys will feature in a simulation.

Acknowledgments

N G gratefully acknowledges supports from Simulation Science Center Clausthal/Göttingen.

M M gratefully acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within GRK 2297 MathCoRe.

Improved energy minimization of iron–carbon systems: on the influence of positioning interstitial atoms

Article metrics

Submit

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Peer review information

Abstract

1. Introduction

2. Methods

2.1. Statically optimal interstitial schemes