Dynamic Simulation of Structural Phase Transitions in Magnetic Iron

TL;DR

This study uses advanced spin-lattice dynamic simulations based on ab initio data to directly demonstrate the finite-temperature phase transitions between bcc and fcc structures in magnetic iron, highlighting the role of magnetic and lattice interactions.

Contribution

It introduces a novel dynamic simulation approach that couples magnetic fluctuations and atomic vibrations to accurately model phase transitions in iron.

Findings

Identified two phase transition points between bcc and fcc phases.

Quantified the free energy difference, maxing at 2 meV, consistent with experimental data.

Demonstrated the importance of magnetic-lattice coupling in phase stability.

Abstract

The occurrence of bcc-fcc ($\alpha$-$\gamma$) and fcc-bcc ($\gamma$-$\delta$) phase transitions in magnetic iron stems from the interplay between magnetic excitations and lattice vibrations. However, this fact has never been proven by a direct dynamic simulation, treating non-collinear magnetic fluctuations and dynamics of atoms, and their coupling at a finite temperature. Starting from a large set of data generated by ab initio simulations, we derive non-collinear magnetic many-body potentials for bcc and fcc iron describing fluctuations in the vicinity of near perfect lattice positions. We then use spin-lattice dynamics simulations to evaluate the difference between free energies of bcc and fcc phases, assessing their relative stability within a unified dynamic picture. We find two intersections between the bcc and fcc free energy curves, which correspond to $\alpha$-$\gamma$ bcc-fcc and $\gamma$-$\delta$ fcc-bcc phase transitions. The maximum fcc-bcc free energy difference over the temperature interval between the two phase transition points is 2 meV, in agreement with other experimental and theoretical estimates.

Tables

Table 1. Table 1: Parameters of the non-magnetic iron potential.

$\varphi$	$-4.483075702293698016e-04$
$r_0^t$	$2.0000000000000000e+00$
$r_1^t$	$2.2000000000000000e+00$
$r_2^t$	$2.6000000000000000e+00$
$r_3^t$	$3.2000000000000000e+00$
$r_4^t$	$3.8000000000000000e+00$
$r_5^t$	$4.6000000000000000e+00$
$r_6^t$	$5.3000000000000000e+00$
$t_0$	$2.5999782982854347e+00$
$t_1$	$2.9319480072508499e+00$
$t_2$	$-2.8388905185188360e+00$
$t_3$	$-1.0267419494754382e-01$
$t_4$	$1.5484736035888333e-02$
$t_5$	$-7.2805743511785065e-02$
$t_6$	$-3.6343523861565924e-03$
$r_0^V$	$2.3254531341916498e+00$
$r_1^V$	$2.3889005055990276e+00$
$r_2^V$	$2.5614990650026459e+00$
$r_3^V$	$2.5615004425308658e+00$
$r_4^V$	$2.8344513929093051e+00$
$r_5^V$	$2.8321879787700808e+00$
$r_6^V$	$2.6382534884695783e+00$
$r_7^V$	$3.4262631740080707e+00$
$r_8^V$	$3.8479639767860356e+00$
$r_9^V$	$3.8515517885908994e+00$
$r_{10}^V$	$4.3740210397021579e+00$
$r_{11}^V$	$4.4054035197078845e+00$
$r_{12}^V$	$4.5503412747697087e+00$
$r_{13}^V$	$4.7731075757035732e+00$
$r_{14}^V$	$5.3000000000000000e+00$
$V_0$	$2.2831054190426084e+01$
$V_1$	$-2.1062362139531867e+01$
$V_2$	$5.6190823955741749e+00$
$V_3$	$8.0795758060570382e+00$
$V_4$	$-8.5213153270399573e+01$
$V_5$	$9.0355710040623180e+01$
$V_6$	$-8.3613137262443793e+00$
$V_7$	$-3.4250845501053456e-01$
$V_8$	$5.2035042290453923e+01$
$V_9$	$-5.1583785613198948e+01$
$V_{10}$	$4.0569674844835752e+00$
$V_{11}$	$-5.0779874829818361e+00$
$V_{12}$	$1.7323802861372730e+00$
$V_{13}$	$-3.5971267571846299e-01$
$V_{14}$	$-1.1478647839739256e-01$

Table 2. Table 2: Parameters for exchange coupling and Landau coefficients

bcc
$J_0$	$1.7613094778950000e-01$
$r_{cut}$	$5.3000000000000000e+00$
$a_0$	$-2.3827723674043900e-01$
$a_1$	$1.2945703172205700e-02$
$a_2$	$-1.1518969922985000e-04$
$b_0$	$1.0600315078586900e-02$
$b_1$	$1.6104913287021000e-03$
$b_2$	$-4.3178188078544200e-05$
fcc
$J_0$	$1.1095507874951400e-01$
$r_{cut}$	$5.3000000000000000e+00$
$b$	$1.6502332463388100e+00$
$c$	$-4.1373722623161200e+00$
$a_0$	$3.1803486683085200e-01$
$a_1$	$6.0141682907976200e-02$
${\rho }_a$	$2.2852502987397700e+01$
$b_0$	$1.4290243674270400e-02$
$b_1$	$0.0000000000000000e+00$
${\rho }_b$	$3.2563330708156800e+01$

Equations

$$⟨n_k⟩+\frac{1}{2}=\frac{1}{\exp \left(\frac{\hslash {\omega }_k}{k_{B}T}\right)-1}+\frac{1}{2}\approx \frac{k_{B}T}{\hslash {\omega }_k}+\frac{1}{12}\left(\frac{\hslash {\omega }_k}{k_{B}T}\right)+...$$

$$\delta {\mathcal{H}}_{um}={\mathcal{H}}^1-{\mathcal{H}}^0.$$

$$⟨\mathcal{O}{⟩}_0=\frac{\int \mathcal{O}\exp (-\beta {\mathcal{H}}^0)d\Omega }{\int \exp (-\beta {\mathcal{H}}^0)d\Omega }$$

$$⟨\mathcal{O}{⟩}_0=\frac{⟨\mathcal{O}\exp (\beta \delta {\mathcal{H}}_{um}){⟩}_1}{⟨\exp (\beta \delta {\mathcal{H}}_{um}){⟩}_1}.$$

$$F^0$$

$$\delta F_{um}$$

$${\mathcal{H}}_{ti}(\lambda )=(1-\lambda ){\mathcal{H}}^0+\lambda {\mathcal{H}}^1,$$

$$\delta {\mathcal{H}}_{ti}$$

$$\delta F_{ti}=F^1-F^0={\int }_0^1⟨\delta {\mathcal{H}}_{ti}{⟩}_{\lambda }d\lambda .$$

$$\lambda (\tau )$$

$${\mathcal{H}}_{HO}=\frac{{\mathbf{p}}^2}{2m}+\frac{1}{2}m{\omega }^2{\mathbf{x}}^2+C,$$

$$F_{HO}=-3k_{B}T\ln \left(\frac{k_{B}T}{\hslash \omega }\right)+C,$$

$${\mathcal{H}}_{LO}=A_{LO}{\mathbf{S}}^2+B_{LO}{\mathbf{S}}^4$$

$$F_{LO}=-k_{B}T\ln \left(4\pi {\int }_0^{\infty }\exp (-\beta {\mathcal{H}}_{LO})S^{2}dS\right)$$

$${\mathcal{H}}_l=\sum _i\frac{{\mathbf{p}}_i^2}{2m}+U(\mathbf{R})$$

$$F_l=F_{HO}-\delta F_l,$$

$$\delta F_l=k_{B}T\ln ⟨\exp (\beta \delta {\mathcal{H}}_l){⟩}_{HO},$$

$$\delta H_l$$

$$\frac{d{\mathbf{R}}_i}{dt}$$

$$\frac{d{\mathbf{p}}_i}{dt}$$

$${\mathbf{F}}_i=-\frac{\partial {\mathcal{H}}_{HO}}{\partial {\mathbf{R}}_i}.$$

$${\mathcal{H}}_{sl}={\mathcal{H}}_l+{\mathcal{H}}_s$$

$${\mathcal{H}}_{HO,s}={\mathcal{H}}_{HO}+{\mathcal{H}}_s,$$

$$\delta {\mathcal{H}}_l$$

$$=$$

$$\frac{d{\mathbf{S}}_i}{dt}=\frac{1}{\hslash }\left[{\mathbf{S}}_i\times {\mathbf{H}}_i\right]+{\gamma }_s{\mathbf{H}}_i+{\mathbf{\xi }}_i$$

$${\mathbf{H}}_i=-\frac{\partial {\mathcal{H}}_{HO,s}}{\partial {\mathbf{S}}_i}.$$

$${\mathbf{F}}_i=-\frac{\partial {\mathcal{H}}_{HO,s}}{\partial {\mathbf{R}}_i}.$$

$$\delta F_l$$

$$=$$

$${\mathcal{H}}_{HO,LO}={\mathcal{H}}_{HO}+{\mathcal{H}}_{LO},$$

$$\delta {\mathcal{H}}_s$$

Full text

I Introduction

Pure iron undergoes bcc-fcc ($\alpha$-$\gamma$) and fcc-bcc ($\gamma$-$\delta$) phase transitions at $T_{\alpha-\gamma}$=1185K and $T_{\gamma-\delta}$=1667K, respectively. They occur in the temperature interval between the Curie temperature $T_{C}=1043$K and the melting temperature $T_{M}$=1811K. These transitions have exceptionally significant practical implications as they are responsible for the formation of martensite in steels, and hence represent the most basic phenomena underpinning steel manufacturing and modern metallurgy. It has long been speculated that $\alpha$-$\gamma$-$\delta$ phase transitions in iron stem from the interplay between magnetic excitations and lattice vibrations. Still, there is no explicit proof, derived from a direct simulation, confirming this assertion. The position is somewhat unsatisfactory as it leaves open the fundamental question of whether, with modern materials modelling concepts and algorithms, it is actually possible to discover a magnetism-driven structural phase transition by exploring the dynamics of a discrete atomistic model.

Hasegawa and PettiforHasegawaPRL1983 investigated the relative stability of bcc, fcc, and hcp phases of iron as a function of temperature and pressure. They concluded that the relative stability of phases was primarily determined by the magnetic free energy contribution. Since they used a single-site spin-fluctuation approximation, which is a mean-field approach where the short-range magnetic order (SRMO) is neglected, the predicted phase diagram was rather qualitative than quantitative.

Recent experimentsNeuhausPRB2014 , which explored phonon dispersion in iron at high temperatures, showed that the stabilization of the high-temperature bcc $\delta$ phase is due primarily to vibrational entropy, whereas the fcc $\gamma$ phase is stabilized by the fine balance between electronic and vibrational entropy contributions. This agrees with calculations performed using Monte Carlo Magnetic Cluster Expansion (MCE) LavrentievCMS2010 ; LavrentievPRB2010 . Although in the MCE the phonon contribution to the free energy is derived from experimental data, the MCE analysis shows that magnetic excitations lower the free energy difference between the bcc and fcc phases $\Delta F^{fcc-bcc}$, and stabilize the $\gamma$ phase. The bcc $\delta$ phase again becomes more stable at higher temperatures, because the lattice vibrations part of the free energy at high temperature is greater than the part associated with magnetic excitations.

Several recent studies of phase stability of iron are based on ab initio calculationsHerperPRB1999 ; SoulairolJPCM2010 ; NguyenPRB2009 ; AlnemratJPCM2014 or the tight-binding Stoner modelSoulairolJPCM2010 ; NguyenPRB2009 ; AutesJPCM2006 . Most of them explore phase stability at 0K, since this is an intrinsic limitation associated with density functional theory (DFT). Treating magnetism in the framework of a tight-binding model is also not trivial since the fully non-collinear version of the Hamiltonian, including both spin and orbital magnetism, has been derived only recentlyCouryPRB2016 . In principle, the Coury HamiltonianCouryPRB2016 should enable fully self-consistent non-collinear magnetic dynamic simulations of atoms and magnetic moments, treated at the electronic scale.

Several approaches have been developed to describe finite-temperature magnetic excitations using ab initio techniques, see Ref. KormannCOSSMS2016, for a review. The disordered local moments (DLM)OkatovPRB2009 ; RazumovkiyPRL2011 ; PolesyaPRB2016 ; MankovskyPRB2013 approximation assumes randomly distributed collinear up and down oriented magnetic moments, to imitate a fully magnetically disordered paramagnetic (PM) state of a material. It has been applied to modeling bcc-fccOkatovPRB2009 and bcc-hcpMankovskyPRB2013 transitions. To treat temperatures lower than $T_{C}$, where iron is in a partially ordered ferromagnetic (FM) state, partial or uncompensated DLM appromationsRazumovkiyPRL2011 ; PolesyaPRB2016 ; MankovskyPRB2013 were proposed, where the net magnetization is constrained to a fixed value, matching experimental observations. This approach reproduces the elastic anisotropy of iron and Fe-Cr alloysRazumovkiyPRL2011 . However, the notion of SRMO remains undefined as the DLM approximation neglects magnetic non-collinearity.

Recently, Leonov et al.LeonovPRL2011 ; LeonovPRB2012 ; LeonovSR2014 investigated the structural stability of iron using a combination of DFT and dynamical mean field theory (DMFT). The treatment involved an explicit consideration of temperature-dependent electron correlations. In the DMFT formalism, excitations associated with electron-electron interactions were treated using a single-site mean-field approximation, and neither collective magnetic excitations nor the SRMO were taken into account.

SRMO can be treated in the random phase approximation (RPA) combined with rescaling KormannPRB2008 . This approach was applied to evaluate magnetic, electronic and phonon contributions to the free energy, and also to assess the pressure dependence of $T_{C}$ KormannPRB2009 in iron. The treatment has also been extended to quantum effects by rescaling the available classical solutions KormannPRB2010 . However, the mean-field nature of the approach gives rise to the predicted value of $T_{C}$ to be higher than the observed value. An alternative approach to modelling SRMO is the spin-wave methodRubanPRB2012 .

SRMO can be simulated using ab initio spin dynamics (SD)AntropovPRL1995 ; AntropovPRB1996 ; IvanovJAP1999 ; FahnleCMS2005 combined with constrained non-collinear calculationsMaPRB2015 ; KurzPRB2004 ; SingerPRB2005 ; UjfalussyJAP1999 ; StocksPMB1998 ; DederichsPRL1984 . However, this is a computationally highly demanding approach, applicable only to relatively small systems. A more practical way of treating SRMO is by deriving parameters from ab intio data and performing atomic scale SDMaPRB2011 ; MaPRB2012 or Monte Carlo RubanPRB2007 ; RubanPRB2013 ; LavrentievCMS2010 ; LavrentievPRB2010 simulations.

In the treatment of structural phase transitions, magnetic and phonon excitations, as well as coupling between them, appear significant. Körmann et al. proposed a spin space averaging procedure KormannPRB2012 ; KormannPRL2014 to evaluate effective interatomic forces at a finite temperature by interpolating wave functions between the FM and PM states and performing statistical averaging over many magnetic microstates. The procedure was applied to phonon spectra and changes in the spectra due to magnon-phonon interactions.

Conventional many-body interatomic potentials for molecular dynamics (MD) treat only the atomic degrees of freedom. Many of the potentials OlssonPRB2005 ; OlssonPRB2006 ; MalerbaJNM2010 ; MarinicaPRL2012 ; LeeJPCM2012 are fitted to ab initio data at 0K. Although some were successfulLeeJPCM2012 in reproducing bcc, fcc, hcp phases at different temperatures and pressures, they required adjusting the energies of fcc and hcp structures to compensate for the absence of an explicit treatment of magnetic effects. Improvements in the functional forms of interatomic potentials were proposedDudarevJPCM2005 ; DerletPMS2007 ; AcklandJNM2006 , but they still did not include magnetic degrees of freedom explicitly, and did not treat magnetic thermal excitations.

There were partially successful attempts to incorporate magnetic states explicitly in MD. For example, Lian et al. LianPRB2015 performed ab initio MD simulations to obtain the phonon dispersion of $\gamma$ and $\delta$ phases over a range of temperatures. However, only the antiferromagnetic state (AFM) was considered as representative of the PM state. Alling et al.AllingPRB2016 performed DLM+MD to study the effect of atomic vibrations on magnetic properties. No equations of motion for magnetic moments were considered, and the dynamics of magnetic excitations was modeled as stochastic spin flips. This leads to conceptual difficulties in the treatment of thermalization of atoms and spins.

In this study, we use spin-lattice dynamics (SLD)MaACP2008 ; MaPRB2008 to study structural phase transitions in magnetic iron. SLD treats the dynamics of lattice and magnetic subsystems within a unified framework. Lattice and spin temperaturesMaPRE2010 can be well controlled through Langevin thermostatsChandrasekharRMP1943 ; KuboRPP1966 ; MaPRB2012 . SLD also treats anharmonicity and the coupling of lattice vibrations to magnetic excitations. SLD is an efficient and versatile simulation approach, and it has been recently applied to a variety of phenomena including the anomalous self-diffusion in ironWenJNM2013 ; WenJNM2014 , magnetic excitations in thin filmsMaPM2009 and a broad range of other magnetic phenomena BeaujouanPRB2012 ; ThibaudeauPASMA2012 ; PereraJAP2014 ; PereraPRB2016 .

An assessment of structural phase stability requires the evaluation of magnetic free energy, where the treatment of SRMO is critically importantKormannPRL2014 ; RubanPRB2012 to the quantitative prediction of phase transition temperatures. The SRMO is treated fully by SLD through the use of dynamic spin equations of motion. In this work, interatomic interaction parameters for SLD simulations are derived from ab initio calculations, and are given in the form of a non-collinear magnetic many-body potential. Using SLD simulations, we are able to evaluate contributions to the free energy from lattice and spin excitations as functions of temperature. By a direct simulation, we find $\alpha-\gamma$ and $\gamma-\delta$ phase transitions, which manifest themselves as changes of sign of $\Delta F^{fcc-bcc}$. The maximum free energy difference between bcc and fcc phases over the temperature interval between the two phase transition points is close to 2 meV.

II Free Energy Calculations

We used two complementary techniques to carry out free energy calculations, the umbrella sampling and thermodynamic integration. Both methods are well established but they have not yet been applied to the treatment of phase transitions in magnetic systems. Brief summaries of the techniques are given below. We also describe our approach to the evaluation of the free energy of a harmonic oscillator and a Landau oscillator in the classical limit, and outline our sampling procedure.

Since we only treat the classical limit, our results are valid at temperatures that are sufficiently high, close to or above approximately one-third of the characteristic temperature for a particular subsystem, the Debye temperature for lattice vibrations and the Curie or the Néel temperature for magnetic excitations. Low temperatures classical results are only given for completeness, and they should not be treated as predictions. A quantum treatment, not considered below, should be applied if one is interested in the accurate low temperature values.

In lattice case, the average number of phonons in mode $k$, $\langle n_{k}\rangle$, is given by the Planck distribution, where at sufficiently high temperatures

[TABLE]

Taking as an estimate $\hbar\omega_{k}\sim k_{B}T_{D}$, we see that classical treatment applies at $T\geq T_{D}/\sqrt{12}\approx T_{D}/3.5$. The Debye temperature of iron is close to 470K. Therefore, one can argue that the temperature range of validity of classical molecular dynamics is defined by the condition $T>135$K. A broadly similar argument can also be applied to magnetic excitations, which we treat classically at temperatures above approximately one third of $T_{C}$ or $T_{N}$.

II.1 Umbrella Sampling

Umbrella samplingTorrieJCP1977 is a biased sampling technique. It is a re-weighting technique for evaluating the free energy difference between a reference and a target state. It is particularly useful for sampling metastable states.

We start by considering two classical Hamiltonians $\mathcal{H}^{0}$ and $\mathcal{H}^{1}$, and their difference:

[TABLE]

The ensemble average of an observable $\mathcal{O}$ with respect to Hamiltonian $\mathcal{H}^{0}$ at a particular temperature $T$ is:

[TABLE]

where $\beta=(k_{B}T)^{-1}$ and $d\Omega$ is an element of volume in classical phase space, which in this instance has 9N dimensions, and includes position vectors of all the atoms, their kinematic momenta, and vectors of all the atomic magnetic moments. Substituting (1) into (2), we arrive at:

[TABLE]

This formula recasts the calculation of an ensemble average of a classic observable $\mathcal{O}$ over the equilibrium defined by Hamiltonian $\mathcal{H}^{0}$ into calculations of ensemble averages of $\mathcal{O}\exp(\beta\delta\mathcal{H}_{um})$ and $\exp(\beta\delta\mathcal{H}_{um})$ over thermodynamic equilibrium defined by another Hamiltonian $\mathcal{H}^{1}$.

This shows a way of evaluating the difference between free energies associated with two classical Hamiltonians $\mathcal{H}^{0}$ and $\mathcal{H}^{1}$. For example, the expression for the free energy corresponding to Hamiltonian $\mathcal{H}^{0}$ can be written as

[TABLE]

Hence, the difference between free energies of two equilibrium configurations defined by Hamiltonians $\mathcal{H}^{0}$ and $\mathcal{H}^{1}$ is

[TABLE]

If one of the free energies $F^{1}$ is known, the other free energy $F^{0}$ can be computed by sampling the phase space with thermodynamic weights defined by $\mathcal{H}^{1}$, and no independent averaging over thermodynamic equilibrium defined by $\mathcal{H}^{0}$ is required.

II.2 Thermodynamic Integration

Another technique for evaluating the difference between free energies is the adiabatic switching thermodynamic integration method KirkwoodJCM1935 ; FrenkelJCP1984 ; CiccottiBOOK1987 . For any two Hamiltonians $\mathcal{H}^{0}$ and $\mathcal{H}^{1}$, we can define a Hamiltonian that is a linear combination

[TABLE]

where $\lambda$ is a switching parameter varying from 0 to 1. The difference between Hamiltonians $\mathcal{H}^{1}$ and $\mathcal{H}^{0}$ equals the derivative of $\mathcal{H}_{ti}$ with respect to $\lambda$.

[TABLE]

The free energy difference between the initial ($\lambda=0$) and final ($\lambda=1$) states can be calculated as an integral over the switching parameter, namely

[TABLE]

Brackets $\langle...\rangle_{\lambda}$ correspond to taking an ensemble average with respect to $\mathcal{H}_{ti}(\lambda)$. We evaluate this average using a dynamic simulation, by imposing the chain rule $d\lambda=(\partial\lambda/\partial t)dt$ and adopting a time-dependent switching function KoningPRE465 :

[TABLE]

where $\tau=t/t_{tot}$, $t$ is the elapsed time and $t_{tot}$ is the total switching time. One can check that if $t=0$ then $\lambda=0$, and if $t=t_{tot}$ then $\lambda=1$.

II.3 Harmonic oscillator and Landau oscillator

We use two reference states for free energy calculations. A harmonic oscillator is used as a reference state for the free energy of the lattice. A Landau oscillator is used as a reference state for the treatment of magnetic excitations.

The Hamiltonian of a three-dimensional harmonic oscillator is

[TABLE]

where $\mathbf{p}$ is the kinematic moment, $\mathbf{x}$ is the displacement, $\omega$ is the frequency, $m$ is the mass and $C$ is a constant. In the classical limit, the free energy can be evaluated analytically as

[TABLE]

where the Planck constant is introduced for dimensional convenience. In what follows, we assume that $\omega$ equals the Debye frequency of iron, $\hbar\omega=k_{B}T_{D}$, where $T_{D}=470$K.

The Landau spin Hamiltonian has the form

[TABLE]

where $A_{LO}$ and $B_{LO}$ are constants, and $\mathbf{S}$ is a dimensionless spin vector. The free energy of magnetic excitations in the Landau approximation can be written as

[TABLE]

where $S$ is the magnitude of $\mathbf{S}$. In this work, we choose $A_{LO}=-1.184$eV and $B_{LO}=0.578$eV to match the spectrum of longitudinal magnetic excitationsMaPRB2012 . The value of $F_{LO}$ is then computed numerically at various temperatures.

Hamiltonians $\mathcal{H}_{HO}$ and $\mathcal{H}_{LO}$ defined above refer to an individual atom and an individual spin. In the calculations below we will use $\mathcal{H}_{HO}$ and $\mathcal{H}_{LO}$ to represent all the atoms and spins $N$, assuming that they are independent of each other.

II.4 Sampling procedure - MD

In an MD simulation, we calculate the free energy using the umbrella sampling. The full lattice Hamiltonian has the form

[TABLE]

where $U(\mathbf{R})$ is the interatomic potential, $\mathbf{R}=\{\mathbf{R}_{i}\}$ are the coordinates of all the atoms, and $\mathbf{p}=\{\mathbf{p}_{i}\}$ are the kinematic momenta.

The free energy of the lattice system can be computed by using equation (4) and sampling over the thermodynamic equilibrium of harmonic oscillators, namely

[TABLE]

where

[TABLE]

and $\delta\mathcal{H}_{l}$ is defined as

[TABLE]

Here $\mathbf{x}_{i}$ is the displacement of atom $i$ from its position in the lattice $\mathbf{R}_{i}^{0}$, i.e. $\mathbf{x}_{i}=\mathbf{R}_{i}-\mathbf{R}_{i}^{0}$. The value of constant $C$ is chosen to minimize the variation of $\delta\mathcal{H}_{l}$. A suitable choice of $C$ helps ensure the numerical stability of umbrella sampling by eliminating large numerical values in the argument of exponential function in Eq. (16).

Sampling is performed using dynamic Langevin thermostat simulatonsChandrasekharRMP1943 ; KuboRPP1966 that generate the correct equilibrium energy distribution, assuming ergodicity. Langevin equations of motion have the form

[TABLE]

where the regular component of the force acting on atom $i$ is

[TABLE]

The damping constant $\gamma_{l}$ and the fluctuating force $\mathbf{f}_{i}$ are related through the fluctuation-dissipation theoremChandrasekharRMP1943 ; KuboRPP1966 , namely $\langle f_{\alpha i}(t)f_{\beta j}(t^{\prime})\rangle=2k_{B}T\gamma_{l}\delta_{\alpha\beta}\delta_{ij}\delta(t-t^{\prime})$, where indexes $\alpha$ and $\beta$ refer to Cartesian coordinates $x,y,z$.

By following the above procedure, we sample over thermodynamic equilibrium defined by the Einstein model for a solid, where the lattice points are ordered as either bcc or fcc lattices. A major advantage of umbrella sampling is that it overcomes the difficulties associated with sampling the spectra of excitations of an unstable structure. For example, the recently developed interatomic potentials for ironMalerbaJNM2010 ; MarinicaPRL2012 predict a stable bcc phase. The fcc phase is unstable, but since sampling is performed over an equilibrium defined by suitably spatially ordered harmonic oscillators, the fact that the crystal structure is unstable has no effect on the stability of the numerical procedure.

II.5 Sampling procedure - SLD

In a spin-lattice dynamic (SLD) simulation, we adopted a two-step approach to free energy calculations. We use the umbrella sampling, which is followed by thermodynamic integration. We write the spin-lattice Hamiltonian as a sum of the lattice and spin parts,

[TABLE]

where the spin part $\mathcal{H}_{s}=\mathcal{H}_{s}(\mathbf{R},\mathbf{S})$ depends on atomic coordinates $\mathbf{R}$ and atomic spin vectors $\mathbf{S}=\{\mathbf{S}_{i}\}$. Since we use various functional forms to represent the Hamiltonians, in what follows we discuss the choice of the specific functional forms adopted in simulations.

First, we apply the umbrella sampling. We define an auxiliary Hamiltonian,

[TABLE]

which is a sum of the harmonic oscillators Hamiltonian and the spin Hamiltonian. The difference between this Hamiltonian and the exact Hamiltonian is

[TABLE]

Notably, this expression is exactly the same as that investigated in the connection with the pure MD analysis. Sampling can again be performed using Langevin thermostat simulations. Since we now also need to take into account magnetic fluctuations, the full set of equations now includes equations of motion for the spinsMaPRB2012 :

[TABLE]

where the damping constant $\gamma_{s}$ and the fluctuation spin force $\bm{\xi}_{i}$ are related by the fluctuation-dissipation theorem $\langle\xi_{\alpha i}(t)\xi_{\beta j}(t^{\prime})\rangle=2k_{B}T\gamma_{s}\delta_{\alpha\beta}\delta_{ij}\delta(t-t^{\prime})$. The effective exchange field acting on spin $i$ is

[TABLE]

Forces in (18) now depend on the orientation of atomic spins

[TABLE]

Similarly to the MD case, we evaluate the difference between the free energies of an equilibrium configuration defined by the spin-lattice Hamiltonian, and a configuration defined by the auxiliary Hamiltonian

[TABLE]

As the second step, we perform thermodynamic integration. We define a reference Hamiltonian

[TABLE]

which is a sum of the harmonic oscillators Hamiltonian and the Landau Hamiltonian. The difference between the reference and the auxiliary Hamiltonians is

[TABLE]

The Hamiltonian required for carrying out thermodynamic integration can be written as

[TABLE]

Langevin equations of motion remain unchanged, but the effective field and the force now depend on the integration parameter $\lambda$, namely

[TABLE]

The free energy difference between the equilibrium states defined by the auxiliary and reference Hamiltonians is

[TABLE]

Combining the results derived using the umbrella sampling and thermodynamic integration, we find the free energy of the equilibrium configuration defined by the spin-lattice Hamiltonian

[TABLE]

We note that this free energy $F_{sl}$ can be represented as a sum

[TABLE]

where

[TABLE]

The above expression has a clear meaning since $F_{l}$ represents a part of the free energy associated primarily with lattice excitations, whereas $F_{s}$ is a part of the free energy derived primarily from spin fluctuations. Since the spin and lattice degrees of freedom are coupled, and we sample through an auxiliary step, it would be inaccurate to interpret $F_{l}$ and $F_{s}$ as independent contributions from the lattice and spin subsystems. However, the two quantities still provide some qualitative insight into the relative magnitude of contributions by the two coupled subsystems to the total free energy.

III Simulations using literature parameterizations

Using the methods described above, we performed MD and SLD simulations, using parameters taken from literature. All the simulations were performed using our MD and SLD program SPILADYMaCPC2016 . We used cubic simulation cells containing 16000 atoms in the bcc case and 16384 atoms in the fcc case. We explored temperatures in the range from 1K to 1400K. For each temperature, we simulated samples at nine different volumes close to an assumed equilibrium volume. A third-order polynomial was then fitted to the calculated free energies. The equilibrium volume was determined from the minimum of the polynomial. All the quantities that we describe below are interpolations corresponding to the equilibrium volume. The guessed equilibrium volume itself was computed using the same method, starting from a larger interval of trial volumes.

In Fig. 1 the free energy computed using SLD simulations is plotted as a function of volume of bcc crystal structure for 300K and 1000K. The polynomial fit interpolates the data points fairly well, although fluctuations are larger at higher temperatures, affecting the accuracy of evaluation of the equilibrium lattice constant. The dotted curves shown in blue indicate the standard deviation of the fitted curve shown in red. It is evaluated using the covariance matrix of the coefficients of the polynomial. The free energy minimum remains accurate at the sub-meV level, as illustrated by the scale of the y-axis.

All the simulation cells were thermalized to equilibrium before sampling. In MD, we take 100,000 data points when performing the umbrella sampling. In SLD, we take 200,000 data points for umbrella sampling, and 0.2ns as the total switching time for the adiabatic switching thermodynamic integration. The Marinica iron potentialMalerbaJNM2010 ; MarinicaPRL2012 was used for both MD and SLD simulations. SLD simulations are based on the spin Hamiltonian of the form MaPRB2008 ; MaPRB2012 :

[TABLE]

where $J_{ij}$ is the exchange coupling function, and $A_{i}$ and $B_{i}$ are the Landau coefficients for atom $i$. The form of $H_{s}$ guarantees that the energy difference between bcc and fcc structures at 0K is the same as in the non-magnetic MD potential case.

We assume that $J_{ij}$ is a pairwise function that depends only on the distance between atoms $i$ and $j$. It has the form $J_{ij}(r)=J_{0}(1-r/r_{c})^{3}\Theta(r_{c}-r)$, where $J_{0}=0.92$eV and $r_{c}=3.75$Å. The value of $J_{0}$ is slightly larger than the one that we derived in Ref. MaPRB2008, to match the experimental value of the Curie temperature $T_{C}$. Values of parameters $A_{i}=-0.744824$ eV and $B_{i}=0.345295$ eV are taken from Ref. MaPRB2012, . We note that the ground state of this Hamiltonian is ferromagnetic regardless of whether the underlying crystal structure is bcc or fcc.

Fig. 2 shows equilibrium lattice constants as functions of temperature, predicted by MD and SLD simulations. Magnetic excitations enhance thermal expansion of both bcc and fcc structures. The curves derived from SLD simulations flatten in the vicinity of $T_{C}$ (Fig. 3). Fluctuation of the curves result primarily from polynomial fitting.

Since we use the same spin Hamiltonian for bcc and fcc cases, they both adopt ferromagnetic ground states at temperatures below $T_{C}$. However, experimental data for fcc iron indicate that it has a relatively low Néel temperature $T_{N}$ of 67KGonserJAP1963 ; JohansonPRB1970 . Although the precise nature of magnetic configuration at temperatures below $T_{N}$ is debatable, the net magnetization is zero. This differs from our simulations, and we will address the issue in the following sections.

In Fig. 4, we plotted the free energy difference between fcc and bcc phases as a function of temperature. We define $\Delta F^{fcc-bcc}=F^{fcc}-F^{bcc}$ in such a way that when this value is positive, bcc phase is stable and vice versa. In MD, $\Delta F^{fcc-bcc}$ is always positive and does not approach zero when temperature increases. Interestingly, if we include magnetic excitations, then $\Delta F^{fcc-bcc}$ decreases significantly. We also show $\Delta F^{fcc-bcc}_{l}=F^{fcc}_{l}-F^{bcc}_{l}$ on the same graph, and $\Delta F^{fcc-bcc}_{s}=F^{fcc}_{s}-F^{bcc}_{s}$ in the inset. The free energy difference predicted by SLD simulations largely originates from the spin subsystem, though it is not sufficient to stabilize the fcc phase at high temperature. We also note that the derivative of $\Delta F^{fcc-bcc}_{s}$ with respect to temperature is small near the Curie temperature $T_{C}$ of the bcc phase. Standard deviations of the calculated free energies are also shown, which are all in the sub-meV level.

Although we find that the magnitude of the magnetic part of the free energy is significant, we see that parameters taken from literature have no chance of success in predicting the bcc-fcc phase transition. The main deficiency of existing parameterizations is that the effective interatomic potentials were all fitted without considering magnetic excitations. Even if we add a spin part to the Hamiltonian in an ad-hoc manner, this still does not fully account for the free energy contribution from the spin system. A better approach to deriving parameters for spin-lattice dynamics simulations is necessary to model magnetic iron on the atomic scale.

IV Parametrization

In what follows, we present a new derivation of parameters for spin-lattice dynamics simulations of bcc and fcc iron. We start by fitting a non-magnetic iron potential, and augment it by the Heisenberg-Landau Hamiltonian.

IV.1 Non-magnetic iron potential

We fitted a non-magnetic iron potential using an interatomic potential fitting program potfitBrommerMSMSE2015 ; BrommerMSMSE2007 ; BrommerPM2006 . It fits a many-body potential to a user-defined functional form. Parameters of the potential are fitted using the force matching methodErcolessiEL1994 , using the total energy and forces taken from ab initio data. All of our ab initio calculations were performed using VASPKressePRB1993 ; KressePRB1994 ; KresseCMS1996 ; KressePRB1996 . We use the GGA-PBEPerdewPRL1996 ; PerdewPRL1997 pseudo-potential with 14 valence electrons. The plane wave energy cutoff is 450 eV.

We first generate ab initio data for the non-magnetic iron. The structures include perfect bcc and fcc lattices, and simulation cells with distortions such as rhombohedral and tetragonal shape, at various volumes. We also produced ab initio data for amorphous structures and structures containing defects. The functional form and parameters of the fitted potential are given in Appendix. Since we are only interested in the energy and free energy differences between bcc and fcc structures, we are not going to discuss other features of our non-magnetic iron potential here.

Fig. 5 shows ab initio energies of non-magnetic bcc and fcc phases at various lattice constants. The minimum energy of the fcc phase is 0.312eV lower than the minimum energy of the bcc phase, in agreement with data from Ref. HerperPRB1999, . Fcc structure is more stable when magnetism is not taken into account. The curves computed using non-magnetic potential appear similar, and the energy difference is 0.317eV. The difference between the absolute values of ab initio data and non-magnetic potential data is due to the different choice of the reference points. Since we take the cutoff distance of the potential as 5.3Å, the energies of all the ab initio data points are reduced by the energy of a perfect bcc structure with lattice constant $a=5.3\times 2/\sqrt{3}=6.1199$Å, where the nearest neighbour distance is 5.3Å. Interatomic forces remain unaffected by this procedure.

IV.2 Magnetic contributions

Fig. 6 shows VASP data for the energy of magnetic bcc and fcc phases at various lattice constants. In the bcc case, we only show the data for the FM collinear ground state. In the fcc case, there are a number of magnetic configurations that all have comparable energies. We show ab initio data for the FM, AFM, and double layer AFM (DLAFM) magnetic configurations. If we impose a constraint and consider only the collinear magnetic configurations, the DLAFM state has the lowest energy. This also agrees with ab initio results given in Ref. HerperPRB1999, .

The energy of the FM bcc phase is now 0.11eV lower than that of the DLAFM fcc phase. We now need to find a way of describing magnetic excitations using a magnetic Hamiltonian added to the non-magnetic Hamiltonian. This magnetic Hamiltonian should also describe interactions between magnetic moments that are not explicitly evident from the VASP data.

There are various ways to describe interactions between magnetic moments. One can use the spin spiral methodOkatovPRB2011 or spin-cluster expansion methodSingerPRL2011 . We choose the spin Hamiltonian in the Heisenberg-Landau formLavrentievCMS2010 ; LavrentievPRB2010 ; MaPRB2012 ,

[TABLE]

where $\mathbf{M}=-g\mu_{B}\mathbf{S}_{i}$ is the magnetic moment of atom $i$, $g=2.0023$ is the electron $g-$factor, $\mu_{B}$ is the Bohr magneton, $A(\rho_{i})$ and $B(\rho_{i})$ are the Landau coefficients that depend on the effective electron density $\rho_{i}$. This is the same $\rho_{i}$ that enters the non-magnetic interatomic potential.

The first step is to calculate values of $J_{ij}$ from the lowest energy state of bcc and fcc phases. The exchange coupling functions are calculated using ab initio electronic structure multiple-scattering formalism. We use the method and program developed by van Schilfgaarde et al. LiechtensteinJMMM1987 ; SchilfgaardeJAP1999 . It is based on the linear muffin-tin orbital approximation combined with Green’s function technique (LMTO-GF). We calculated values of parameters $J_{ij}$ involving various neighbours over a range of variation of the lattice constant. Fig. 7 shows the calculated values and the fitted curves. Ab initio data for the bcc case are smoother, whereas the data for the fcc cases are more scattered. A possible reason is that magnetic configuration of collinear FM state in the bcc phase is fundamentally simpler than the DLAFM state of the fcc phase. To match the data, we used different functional forms for the bcc and fcc cases. The functional forms and values of parameters are given in Appendix.

We now evaluate the Landau coefficients. We define a temporary Hamiltonian with no magnetic moment interactions, assuming that on-site magnetic moments can be treated as order parameters, i.e.

[TABLE]

Since we know the difference between energies of magnetic and non-magnetic configurations (Fig. 8), and also the magnitude of magnetic moments on each atoms (Fig. 9) as functions of lattice constant, we can identify the energy difference per atom and the magnitude of magnetic moment as

[TABLE]

or $M_{0}=0$ if the non-magnetic state is more stable. We obtain values of $A^{\prime}$ and $B^{\prime}$ at various lattice constants by solving the above equations.

We now need to relate the values of the Landau parameters to the values of parameters characterizing the Hamiltonian that describes interacting magnetic moments. We equate Eq. 41 and 42, and find

[TABLE]

Using the fitted function $J_{ij}$ and considering magnetic configurations of perfect crystals, we find the values of Landau coefficients $A$ and $B$. These coefficients are plotted as functions of the effective electron density $\rho_{i}$ for perfect crystals (Fig. 10 and 11). Again, we use different functional forms for fitting results for bcc and fcc cases, to match various features of the curves. The functional forms and numerical parameters are given in Appendix. The strong scatter of values for $B(\rho_{i})$ is due to the small value of $M_{0}$ corresponding to the small size of the simulation box.

Using the above procedure, we generated several sets of parameters, which were derived using different methods and have various functional forms. The parameters that have been selected are those that match experimental results well. Although all of the parameters produce qualitatively similar predictions, Landau coefficients may need to be adjusted through the choice of fitting intervals to achieve sub-meV accuracy of free energy calculations. In all cases, ab initio data provide the foundation for the fitting procedure.

V Structural Phase Transitions

We calculated the free energies of bcc and fcc phases using the above new sets of parameters. We performed both MD and SLD simulations. In MD, we used the non-magnetic potential, and performed umbrella sampling calculations with 100,000 data points. In SLD, we used the non-magnetic potential with Heisenberg-Landau Hamiltonian, and used the two step approach in the free energy calculations. We took 300,000 data points for umbrella sampling, and 0.2ns as the total switching time in the adiabatic switching thermodynamic integration. The magnetic configuration of bcc and fcc cases are initialized as FM and DLAFM states, respectively. We explored a large temperature range from $1\times 10^{-5}$K to 2000K. All that samples are thermalized to equilibrium before sampling.

The most significant results of this paper are illustrated in Fig. 12 and 13. Fig. 12 shows the calculated free energies of bcc and fcc phases at the equilibrium volume. In the MD case, the free energy of the fcc phase is always lower than the free energy of the bcc phase. In SLD, the bcc phase has lower free energy initially, but the curves corresponding to bcc and fcc phases approach each other at higher temperature. There are two intersections between the curves, which can be seen if one follows the difference $\Delta F^{fcc-bcc}$ plotted in Fig. 13. The curve crosses the zero line at around 1130K and 1600K. These temperatures are close to the experimentally observed values of $\alpha$-$\gamma$ and $\gamma$-$\delta$ phase transitions at $T_{\alpha-\gamma}$=1185K and $T_{\gamma-\delta}$=1667K, respectively. The free energy difference at the minimum is close to 2meV. This value agrees with the MCELavrentievCMS2010 ; LavrentievPRB2010 and RPAKormannCOSSMS2016 results, and shows that $\alpha-\gamma-\delta$ transitions are associated with fairly small free energy differences between the competing phases, of the order of 1 meV.

The reason for the difference between MD and SLD simulations may be understood by considering free energy contributions from magnetic excitations and lattice vibrations. Fig. 14 shows plots of $\Delta F^{fcc-bcc}_{l}$ and $\Delta F^{fcc-bcc}_{s}$ as functions of temperature. For comparison, we also plotted $\Delta F^{fcc-bcc}_{MD}$ calculated using the non-magnetic potential. The value of $\Delta F^{fcc-bcc}_{l}$ is similar to $\Delta F^{fcc-bcc}_{MD}$. Both of them increase relatively linear as functions of temperature. On the other hand, $\Delta F^{fcc-bcc}_{s}$ decreases initially but flattens out near the Curie temperature $T_{C}$ of the bcc phase. Hence the change of $\Delta F^{fcc-bcc}_{s}$ is mainly due to the difference in the degree of disorder associated with magnetic configurations involved. If the temperature is higher than $T_{C}$, the bcc phase is in the paramagnetic state, where the long range magnetic order vanishes. Magnetic configurations of fcc and bcc phases become similar, making the entropy difference smaller. This interpretation is consistent with recent experimental findingsNeuhausPRB2014 on phonon dispersion of iron. Experimental data suggest that the $\gamma$ phase forms as a result of interplay between electronic and vibrational contributions to entropy, whereas the $\delta$ phase is due primarily to the contribution of vibrational entropy. The two crossing points result from the interplay between free energy contributions derived from lattice and spin excitations. The standard deviation of the free energy remains in the sub-meV level.

The magnetization curve of iron is plotted in Fig. 15. For the bcc phase, the existing model predicts a relatively high Curie temperature at around 1300K. For the fcc case, the predicted magnetization is zero. We investigated the Néel temperature $T_{N}$ of the fcc phase by changing the value of energy per atom as a function of temperature (Fig. 16). We thermalized two simulation cells to 165K and 205K. Then, we gradually increased the temperature of the simulation cell equilibrated at 165K to 205K. This was performed by increasing the thermostat temperature linearly over the time interval of 4 ns. Then, we decrease the temperature of the simulation cell equilibrated at 205K to 165K. We filtered the output data by averaging over every 1000 data points. We see that the magnetic subsystem undergoes a first order transition at around 185K. A jump of about 0.01 eV in the magnetic energy can be observed. The change of the total energy is fully accounted for by the magnetic energy. The calculated $T_{N}$ is higher than the experimental value of 67K. However, such value is obtained from experiments on small particlesGonserJAP1963 ; JohansonPRB1970 where the role of local stresses is unclear. Besides, the temperatures of the $\alpha-\gamma$ and $\gamma-\delta$ phase transitions are significantly higher than $T_{N}$.

We also show equilibrium lattice constants of bcc and fcc phases predicted by MD and SLD simulations in Fig. 17. The addition of the spin Hamiltonian changes the value of the equilibrium lattice constant even at 0K. Fundamentally, this agrees with ab initio data shown in Fig. 5 and 6. The spin subsystem of the material affects mechanic properties through its contribution to the total free energy.

VI Conclusions

Starting from a large amount of ab initio data, we fitted non-magnetic many-body potentials and Heisenberg-Landau Hamiltonians for bcc and fcc iron. We performed free energy calculations using umbrella sampling and thermodynamics integration. The free energy has been sampled by molecular dynamics and spin-lattice dynamics simulations. Our method provides a reasonably consistent way of assessing the phase stability of magnetic iron within a unified dynamic picture. It treats both magnetic excitations and lattice vibrations and their coupling self-consistently. The bcc-fcc ($\alpha$-$\gamma$) and fcc-bcc ($\gamma$-$\delta$) phase transitions in magnetic iron are reproduced using newly fitted potential and parameters. The structural phase stability of magnetic iron is governed by non-collinear magnetic excitations and lattice vibrations, in agreement with other experimental and theoretical results. The maximum free energy difference between bcc and fcc phases is about 2 meV.

Appendix A Non-magnetic iron potential

The functional form of the interatomic potential broadly follows the conventional embedded atom method (EAM) representation

[TABLE]

where $\mathbf{R}_{i}$ is the position of atom $i$, $\rho_{i}$ is the effective electron density and $V_{ij}$ is a pairwise function that depends only on the distance between atoms $i$ and $j$. The many-body part of the potential takes the same form as that proposed by Mendelev et al.MendelevPM2003 and Ackland et al.AcklandJPCM2004

[TABLE]

where $\phi$ is a parameter. The effective electron density $\rho_{i}$ is defined in a slightly different way from the conventional EAM potential. We write

[TABLE]

where $t_{ij}=t_{ij}(R_{ij})$ is a pairwise hopping integral, which we take a function of the distance between the atoms $R_{ij}$. We note that the derivative of $\rho_{i}$ with respect to $R_{ij}$ is

[TABLE]

$t_{ij}$ and $V_{ij}$ are given by the third-order splines

[TABLE]

where $n$ are knots, $t_{n}$, $V_{n}$ are parameters with dimensionality eVÅ*-3*, and $r_{n}^{t}$ and $r_{n}^{V}$ are given in Å  units. Their values are given in Table 1.

Appendix B Exchange coupling and Landau coefficients

The exchange coupling function and Landau coefficients, expressed as functions of electron density, have different functional forms for bcc and fcc cases. In the bcc case, we use the following form

[TABLE]

In the fcc case, we take

[TABLE]

The units of $J_{ij}$, $A$ and $B$ are eV$\mu_{B}^{-2}$, eV$\mu_{B}^{-2}$ and eV$\mu_{B}^{-4}$, respectively. The cutoff distance $r_{cut}$ is in Å  units. All the parameters are listed in Table 2.

Acknowledgements.

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053 and from the RCUK Energy Programme [grant number EP/P012450/1]. To obtain further information on the data and models underlying this paper please contact PublicationsManager@ccfe.ac.uk*. The views and opinions expressed herein do not necessarily reflect those of the European Commission. We also acknowledge EUROFusion for the provision of Marconi supercomputer facility at CINECA in Italy. The authors are grateful to M.-C. Marinica for stimulating discussions.