Solid-liquid Interfacial Energy for Fe–Cr Alloy under Temperature Gradient from Molecular Dynamics Simulation

Abstract

The solid-liquid interfacial energy of Fe–Cr alloy under temperature gradient is investigated by molecular dynamics (MD) simulations in conjunction with a capillary fluctuation method including the effect of temperature gradient. It is revealed from the MD simulation that fluctuation of the solid-liquid interface decreases with increasing temperature gradient. This results in a large value of the solid-liquid interfacial energy under large temperature gradient. On the other hand, there is a competing effect reducing the solid-liquid interfacial energy with increasing temperature gradient in the formulation of the capillary fluctuation method including the effect of temperature gradient. As a result, the solid-liquid interfacial energy doesn’t change significantly at small temperature gradient. Moreover, it is confirmed that the solid-liquid interfacial energy of Fe–Cr alloy decreases with increasing Cr composition at Fe-rich composition regardless of the temperature gradient.

1. Introduction

It is essential to understand interfacial properties of solidification microstructure under various temperature conditions in order to fabricate practical metallic materials with desired properties since solidification microstructure is strongly affected by temperature conditon.^1,2,3,4,5,6) Especially, recent progress in additive manufacturing process shed light on the importance of interfacial properties under large temperature gradient.^7,8,9,10) For example, the microstructure obtained via additive manufacturing process is strongly influenced by the direction of heat flux and the large temperature gradient at solid-liquid interface. However, it is not straightforward to obtain such interfacial properties under large temperature gradient experimentally. Therefore, it is expected that computational approaches such as molecular dynamics (MD) simulation derive these properties.

There are several techniques for the estimation of solid-liquid interfacial energy on the basis of the MD simulation.¹¹⁾ For example, a capillary fluctuation method¹²⁾ is a popular technique, in which interfacial stiffness is extracted from the capillary fluctuation spectrum of roughness of the solid-liquid interface. A classical nucleation theory (CNT)-based technique¹³⁾ is also useful, in which the Gibbs-Thomson effect is considered for the threshold temperature dividing shrinkage and expansion of small solid particles in the undercooled melt. Actually, there are many reports on interfacial energy and mobility of pure metals^{12,13,14,15,16,17,18)} based on MD simulations in conjunction with these techniques. On the other hand, it is not straightforward to estimate the solid-liquid interfacial energy for alloys compared to pure metals since the solute partition at the solid-liquid interface should be taken into account properly for a given interatomic potential.⁸⁾ Several techniques have been proposed for tackling this problem. For example, a semi-grand canonical Monte-Carlo (SGCMC) simulation,^19,20,21,22) in which equilibrium composition is obtained for the given chemical potential difference, is used to derive solidus and liquidus compositions for binary alloys. Equilibrium solidus and liquidus compositions in a solid-liquid biphasic system can be also derived by the Metropolis Monte-Carlo (MMC) simulation in conjunction with the atoms swap technique.²³⁾ Moreover, we have proposed a simple technique, in which solidus and liquidus compositions are directly derived from equilibrium compositions around the solid–liquid interface in a long-time MD simulation.²³⁾ By utilizing above techniques, there are several recent reports on the interfacial energy of alloys.^{24,25,26,27,28)}

One of the remaining problems is how to treat the effect of temperature gradient on the solid-liquid interfacial energy. Recently, Brown and coworkers²⁹⁾ proposed a derivation process of the solid-liquid interfacial energy under the large temperature gradient. They modified the formulation of the capitally fluctuation method and applied this method for rapid solidification of Al system. Here, the effect of temperature dependent on the solid-liquid interfacial energy is discussed for iron-based alloy using this technique. In this study, Fe–Cr binary alloy is examined as one of the most practical iron-base alloy. Particularly, Fe–Cr alloys of Fe-rich composition (Fe-5at%Cr and Fe-19.1at%Cr) are examined since the effect of solute partition is very small at these compositions. For comparison, the calculation for pure Fe is also performed.

2. Simulation Methodology

The solid-liquid interface energy of Fe–Cr alloy under temperature gradient is investigated by MD simulations in conjunction with a capillary fluctuation method including the effect of temperature gradient.²⁹⁾ The simulation methodology basically follows our previous study²⁷⁾ except for the installation of temperature gradient in the MD simulation. MD simulations are performed by the Large-scale Atomic/Molecular Massive Parallel Simulator (LAMMPS).³⁰⁾ The embedded atom method (EAM) potentials fitted by Bonny et al.³¹⁾ is employed as the interatomic potential of Fe–Cr binary alloy. Parameter files of the EAM potential for simulations by LAMMPS are employed from the interatomic potential repository (IPR) at National Institute of Standards and Technology (NIST).³²⁾ The velocity-Verlet method is used to integrate the classical equation of motion with a time step of 1.0 fs. The Nose-Hoover thermostat and barostat^33,34) are employed to control temperature and pressure in MD simulations. In addition, the Langevin thermostat³⁵⁾ is employed for introducing temperature gradient in the system. The open visualization tool (OVITO)³⁶⁾ is used for the visualization of simulation result and the post-analyses. The polyhedral template matching (PTM) technique³⁷⁾ is employed for classifying atoms in crystalline or liquid, which is suitable for the classification of atoms at high temperature. A root-mean-square deviation (RMSD) cutoff parameter of the PTM is set to 0.155.²⁷⁾

Quasi two-dimensional (2D) solid-liquid biphasic systems²⁷⁾ are employed to obtain atomistic configuration of the solid-liquid interface for pure Fe, Fe-5at%Cr and Fe-19.1at%Cr, which are prepared as follows. Firstly, a liquid structure is prepared by annealing a bcc crystal (194400 atoms in a cell of 514.8 × 8.58 × 514.8 ų) at 5000 K for 10 ps with the NVT (the number of atoms, volume and temperature)-constant ensemble using the Nose-Hoover thermostat. Then, a solid bcc crystal of the same size is connected with the liquid facing (100) plane at the interface. The short side of the system is set to be [010] direction (notated as (100)[010] orientation hereafter). The energy minimization is then performed for the combined structure to avoid unexpected proximity of atoms at the solid-liquid interface. In the same manner, biphasic systems with (110)[001] and (110)[110] interfaces are prepared. Details of configuration of biphasic systems are listed in Table 1. All prepared biphasic systems are relaxed with the NPT (the number of atoms, pressure and temperature)-constant ensemble for 25 ps at 0 Pa using the Nose-Hoover thermostat. Then, biphasic systems are relaxed with the NPH (number of atoms, pressure and enthalpy)-constant ensemble using the Langevin thermostat for subsequent 9975 ps at 0 Pa to introduce the temperature gradient (G). The Langevin thermostats are set at the center of solid and liquid regions with 2 nm except for the case of G = 119.0 K/nm. For G = 119.0 K/nm, thermostats are set at the closer position from the solid-liquid interface to introduce large temperature gradient. Different temperatures at thermostats in solid (T_c) and liquid (T_H) regions introduce a temperature gradient in the system. Target temperatures of thermostats and corresponding temperature gradient are listed in Table 2.

Table 1. Size of simulation cells for three interface orientations. b and W are thickness and length of the interface, respectively.

Orientation	Number of Atoms	W [Å]	b [Å]
(100)[010]	388800	514.800	8.580
(110)[001]	387096	513.671	8.580
(110)[110]	548640	514.800	12.134

Table 2. Thermostat temperatures (T_C, T_H) and corresponding temperature gradient (G).

TC [K]	TH [K]	G [K/nm]
1550	1950	7.7
1300	2200	17.3
900	2500	30.8
500	3000	48.1
100	4000	75.0
500	3000	119.0

3. Results and Discussion

Figure 1 shows the snapshot of atomic configuration after 10000 ps calculation for the pure Fe system under G = 30.8 K/nm. Although the solid-liquid interface moves at the initial stage of calculation, movement of the solid-liquid interface stops by the end of the calculation, which represents the temperature distribution in the system converges at the end of calculation. Then, subsequent calculation is continued for 400 ps and atomic configuration of the system is extracted at every 1 ps (i.e., 401 snapshots). After the interval of 1000 ps, atomic configuration is again extracted at every 1 ps during 400 ps. Same procedure of the sampling is repeated for three times. Interfacial stiffness and solid-liquid interfacial energy are estimated three times separately for each condition from atomic configuration sampled as above. A series of above procedures are performed for pure Fe, Fe-5at%Cr and Fe-19.1at%Cr with various temperature gradients listed in Table 2.

Fig. 1.

Snapshot of atomic configuration after 10000 ps calculation for pure Fe system under G = 30.8 K/nm. The short side of the system is set to be [010] direction. Arrows and dotted lines represent the position of Langevin thermostats. and Blue and white atoms represents atoms with body-centered cubic (BCC) and liquid configurations, which are defined by polyhedral template matching technique. (Online version in color.)

Then, the position of the solid-liquid interface is extracted as follows. The MD simulation cell is divided into two-dimensional grids of 128 × 256. PTM configurations of all atoms (BCC or liquid) are assigned to the nearest each grid point and the ratio of BCC atoms to the number of assigned atoms is calculated for all grids. The solid-liquid interface is defined at the position between the grids with the minimum percentage of the solid and the maximum percentage of the liquid along the direction perpendicular to the solid-liquid interface. Technical detail of the technique is shown in our previous report.²⁷⁾ Figure 2 shows the solid-liquid interfaces for pure Fe with (100)[010] orientation under temperature gradient of 0, 30.8 and 75.0 K/nm, respectively. The solid-liquid interface defined as above approach is drawn as a red curve. Note that the red curve agrees well with the solid-liquid interface predicted from the visual observation of atomistic configuration. Coordinates of red curves are used for the capillary fluctuation method below. It is expected from Fig. 2 that the solid-liquid interface becomes smooth with increasing temperature gradient. This effect on the solid-liquid interfacial energy will be discussed later.

Fig. 2.

Snapshots of solid-liquid interface in the solid-liquid biphasic system of pure Fe with (100)[010] orientation under temperature gradients (G), 0, 30.8 and 75.0 K/nm. Blue and white atoms represents atoms with BCC and liquid configurations. The solid-liquid interface is drawn in the red curve. (Online version in color.)

The capillary fluctuation method with a modified term is employed to estimate the solid-liquid interfacial energy under temperature gradient. According to Brown and coworkers,²⁹⁾ the interface stiffness can be estimated from the amplitude of the fluctuation of the solid-liquid interface under temperature gradient on the basis of the following equation,   

$$⟨{\left|A(k)\right|}^2⟩=\frac{k_B{T}_{eq}}{bWS(k^2+a^{-2})},$$(1)

  

$$a=\sqrt{\frac{{\gamma }_{G0}{T}_{eq}}{LG}},$$(2)

where S is the solid-liquid interfacial stiffness, A(k) is the Fourier amplitude of the interface height profile, k is the wave number, k_B is Boltzmann’s constant and T_eq is the equilibrium temperature at which the solid-liquid interface does not move (i.e., melting point for pure metals and equilibrium temperature for alloys). b and W are the thickness and length of the interface, respectively. The brackets represent averages over the number of sampling times. a is additional contributions to the fluctuation spectrum that are independent of stiffness for a stationary solid-liquid interface under temperature gradient,²⁹⁾ which is derived from Karma’s theory on fluctuations in solidification.³⁸⁾ G is the temperature gradient, γ_G0 is the solid-liquid interfacial energy at no temperature gradient (i.e., G = 0), and L is the latent heat at the equilibrium temperature. The latent heat is separately calculated as follows. A bcc crystal consisting of 20 × 20 × 20 unit cells (16000 atoms) is annealed at its equilibrium temperature provided in the literature²³⁾ with the NPT-constant ensemble. A liquid structure obtained by pre-anneal of the bcc crystal at 5000 K for 10 ps is annealed at same temperature. The latent heat is derived from the difference of enthalpy for solid and liquid systems. These calculations are performed for pure Fe, Fe-5at%Cr and Fe-19.1at%Cr, respectively. Estimated values of the latent heat and corresponding equilibrium temperature are listed in Table 3.

Table 3. Latent heat (ΔH) and corresponding equilibrium temperature (T_eq).

Composition	ΔH [eV/atom]	Teq [K]
Pure Fe	0.161	1769
Fe-5at%Cr	0.156	1770
Fe-19.1at%Cr	0.149	1790

Figure 3 shows the relationship between $k_B{T}_{eq}/bW⟨{\left|A(k)\right|}^2⟩$ and k² + a⁻² for solid-liquid interfaces with (100)[010] orientation extracted from MD simulations under various temperature gradients for pure Fe system. The interfacial stiffness can be estimated from the slope of the fitted line according to Eq. (1). From the fitted lines for various temperature gradients, it is confirmed that the interfacial stiffness basically increases with increasing temperature gradient. Interfacial stiffness of the solid-liquid interface in (110)[001] and (110)[110] orientations is estimated in the same manner. Once the interfacial stiffness of solid-liquid interface for three different orientations are obtained, solid-liquid interfacial energy can be derived using following relation^27,39)   

$$S={\gamma }_0\left(1-\frac{18}{5}{\epsilon }_1-\frac{80}{7}{\epsilon }_2\right),\mathrm{\hspace{0.17em} }(\mathrm{f}\mathrm{o}\mathrm{r}\text{ (}100)[010]\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n})$$(3)

  

$$S={\gamma }_0\left(1+\frac{39}{10}{\epsilon }_1+\frac{155}{14}{\epsilon }_2\right),\mathrm{\hspace{0.17em} }(\mathrm{f}\mathrm{o}\mathrm{r}\text{ (}110)[001]\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n})$$(4)

  

$$S={\gamma }_0\left(1-\frac{21}{10}{\epsilon }_1+\frac{365}{14}{\epsilon }_2\right),\mathrm{\hspace{0.17em} }(\mathrm{f}\mathrm{o}\mathrm{r}\text{ (}110)[0\overline{1}0]\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n})$$(5)

where γ₀ represents the average solid-liquid interfacial energy and ε₁ and ε₂ are the anisotropy parameters. Practically, γ₀, ε₁ and ε₂ are specifically obtained by assigning values of interfacial stiffness for three different orientations into Eqs. (3), (4) and (5). These equations are derived from a low-order expansion consistent with cubic symmetry.³⁹⁾

Fig. 3.

Relationship between $k_B{T}_{eq}/bW⟨{\left|A(k)\right|}^2⟩$ and k² + a⁻² obtained for solid-liquid interfaces with (100)[010] orientation extracted from MD simulations under various temperature gradients for pure Fe system. The slopes of the solid lines represent the interfacial stiffness. (Online version in color.)

Figure 4 shows average solid-liquid interfacial energy as a function of temperature gradient for pure Fe, Fe-5at%Cr and Fe-19.1at%Cr, respectively. Error bar represents standard deviation of three calculations. At a glance, the solid-liquid interfacial energy decreases with increasing Cr composition in the range of our simulation (i.e., Fe-rich composition). This agrees with our previous simulation without consideration of temperature gradient.²⁷⁾ The solid-liquid interfacial energy generally increases with increasing temperature gradient. It is caused by the fact that the fluctuation of the solid-liquid interface decreases with increasing temperature gradient as shown in Fig. 2. This results in the increase of the slope of fitted lines in Fig. 3, that is, the increase of interfacial stiffness and solid-liquid interfacial energy as a result. On the other hand, the solid-liquid interfacial energy doesn’t change significantly at small temperature gradient (G < 30.8 K/nm). This is caused by another effect in the formulation of the capillary fluctuation method including temperature gradient. That is, G exists in the denominator of right side of Eq. (2). Therefore, the interfacial stiffness is inversely proportional to G at the constant Fourier amplitude. This factor decreases the solid-liquid interfacial energy with increasing temperature gradient. Therefore, it is considered that the solid-liquid interfacial energy doesn’t change significantly at small temperature gradient by competition of two competing factors. In order to confirm above discussion, plots in Fig. 4 are fitted by the following function   

$$\gamma =\frac{1}{pG+q}+rG+s$$(6)

where p, q, r, and s are fitting coefficients. Equation (6) includes two competing factors as a function of G. The first term of right side represents the effect from Eq. (2) and the second term represents the result from MD simulation that the solid-liquid interface becomes smooth with increasing temperature gradient. Since fitting curves in Fig. 4 basically follows the plots, it is considered that two competing as a function of G causes a plateau in the solid-liquid interfacial energy at small temperature gradient.

Fig. 4.

Solid-liquid interfacial energy, γ₀ as a function of temperature gradient for pure Fe, Fe-5at%Cr and Fe-19.1at%Cr. Error bar represents standard deviation of three calculations. Solid curves represent fitted curves using Eq. (6). (Online version in color.)

4. Conclusions

In this study, the solid-liquid interfacial energy of Fe–Cr binary alloy under temperature gradient is investigated on the basis of the MD simulation in conjunction with the capillary fluctuation method. In general, the solid-liquid interfacial energy becomes large with increasing temperature gradient since the solid-liquid interface becomes smooth with increasing temperature gradient at all compositions. However, the solid-liquid interfacial energy doesn’t change significantly at small temperature gradient since there are two competing effects on the interfacial energy as a function of temperature gradient. Moreover, it is confirmed that the solid-liquid interfacial energy of Fe–Cr alloy decreases with increasing Cr composition at Fe-rich composition. In summary, it is significant in this study to reveal the effect of temperature gradient on the solid-liquid interfacial energy for practical alloys since it is not straightforward to measure it directly from experimental approach. The solid-liquid interfacial energy for various practical iron-based alloys will be examined in the next step.

Acknowledgement

This work was supported by Grant-in-Aid for Scientific Research (B) [No. 16H04490] from Japan Society for the Promotion of Science (JSPS), Japan and the 26th ISIJ Research Promotion Grants from the Iron and Steel Institute of Japan. Also, this work was partially supported by 19th Committee (Steelmaking), JSPS, Japan.

References

• 1)   W.  Kurz and  D. J.  Fisher: Fundamentals of Solidification, Trans Tech Publications, Stafa-Zurich, (1998), 13.
• 2)   J. A.  Dantzig and  M.  Rappaz: Solidification, EPFL Press, Lausanne, (2009), 245.
• 3)   Y.  Shibuta,  M.  Ohno and  T.  Takaki: JOM, 67 (2015), 1793. https://doi.org/10.1007/s11837-015-1452-2
• 4)   T.  Ohashi and  W. A.  Fischer: Tetsu-to-Hagané, 61 (1975), 3077 (in Japanese). https://doi.org/10.2355/tetsutohagane1955.61.15_3077
• 5)   O.  Haida and  T.  Emi: Tetsu-to-Hagané, 63 (1977), 1564 (in Japanese). https://doi.org/10.2355/tetsutohagane1955.63.9_1564
• 6)   H.  Mizukami,  K.  Hayashi,  M.  Numata and  A.  Yamanaka: Tetsu-to-Hagané, 97 (2011), 457 (in Japanese). https://doi.org/10.2355/tetsutohagane.97.457
• 7)   I.  Gibson,  D. W.  Rosen and  B.  Stucker: Additive Manufacturing Technologies, Springer, New York, (2010), 33.
• 8)   S. L.  Sing,  L. P.  Lam,  D. Q.  Zhang,  Z. H.  Liu and  C. K.  Chua: Mater. Charact., 107 (2015), 220. https://doi.org/10.1016/j.matchar.2015.07.007
• 9)   P. C.  Collins,  D. A.  Brice,  P.  Samimi,  I.  Ghamarian and  H. L.  Fraser: Annu. Rev. Mater. Res., 46 (2016), 63. https://doi.org/10.1146/annurev-matsci-070115-031816
• 10)   J. H.  Martin,  B. D.  Yahata,  J. M.  Hundley,  J. A.  Mayer,  T. A.  Schaedler and  T. M.  Pollock: Nature, 549 (2017), 365. https://doi.org/10.1038/nature23894
• 11)   Y.  Shibuta: Mater. Trans., 60 (2019), 180. https://doi.org/10.2320/matertrans.ME201712
• 12)   J. J.  Hoyt,  M.  Asta and  A.  Karma: Phys. Rev. Lett., 86 (2001), 5530. https://doi.org/10.1103/PhysRevLett.86.5530
• 13)   Y.  Watanabe,  Y.  Shibuta and  T.  Suzuki: ISIJ Int., 50 (2010), 1158. https://doi.org/10.2355/isijinternational.50.1158
• 14)   D. Y.  Sun,  M.  Asta and  J. J.  Hoyt: Phys. Rev. B, 69 (2004), 174103. https://doi.org/10.1103/PhysRevB.69.174103
• 15)   J.  Monk,  Y.  Yang,  M. I.  Mendelev,  M.  Asta,  J. J.  Hoyt and  D. Y.  Sun: Modell. Simul. Mater. Sci. Eng., 18 (2010), 015004. https://doi.org/10.1088/0965-0393/18/1/015004
• 16)   R.  Hashimoto,  Y.  Shibuta and  T.  Suzuki: ISIJ Int., 51 (2011), 1664. https://doi.org/10.2355/isijinternational.51.1664
• 17)   S. K.  Deb Nath,  Y.  Shibuta,  M.  Ohno,  T.  Takaki and  T.  Mohri: ISIJ Int., 57 (2017), 1774. https://doi.org/10.2355/isijinternational.ISIJINT-2017-221
• 18)   N. T.  Brown,  E.  Martinez and  J.  Qu: Comput. Mater. Sci., 168 (2019), 65. https://doi.org/10.1016/j.commatsci.2019.05.059
• 19)   J. J.  Hoyt,  J. W.  Garvin,  E. B.  Webb III and  M.  Asta: Model. Simul. Mater. Sci. Eng., 11 (2003), 287. https://doi.org/10.1088/0965-0393/11/3/302
• 20)   C. A.  Becker,  M.  Asta,  J. J.  Hoyt and  S. M.  Foiles: J. Chem. Phys., 124 (2006), 164708. https://doi.org/10.1063/1.2185628
• 21)   G.  Bonny,  R. C.  Pasianot,  E. E.  Zhurkin and  M.  Hou: Comput. Mater. Sci., 50 (2011), 2216. http://doi.org/10.1016/j.commatsci.2011.02.032
• 22)   K.  Ueno and  Y.  Shibuta: IOP Conf. Ser. Mater. Sci. Eng., 529 (2019), 012037. https://doi.org/10.1088/1757-899X/529/1/012037
• 23)   K.  Ueno and  Y.  Shibuta: Materialia, 4 (2018), 553. https://doi.org/10.1016/j.mtla.2018.11.011
• 24)   P.  Saidi,  R.  Freitas,  T.  Frolov,  M.  Asta and  J. J.  Hoyt: Comput. Mater. Sci., 134 (2017), 184. http://doi.org/10.1016/j.commatsci.2017.03.044
• 25)   C.  Qi,  B.  Xu,  L. T.  Kong and  J. F.  Li: J. Alloy. Compd., 708 (2017), 1073. https://doi.org/10.1016/j.jallcom.2017.03.077
• 26)   J. J.  Hoyt,  S.  Raman,  N.  Ma and  M.  Asta: Comput. Mater. Sci., 154 (2018), 303. https://doi.org/10.1016/j.commatsci.2018.07.050
• 27)   K.  Ueno and  Y.  Shibuta: Comput. Mater. Sci., 167 (2019), 1. https://doi.org/10.1016/j.commatsci.2019.05.023
• 28)   S.  Kavousi,  B. R.  Novak,  M. I.  Baskes,  M. A.  Zaeem and  D.  Moldovan: Model. Simul. Mater. Sci. Eng., 28 (2020), 015006. https://doi.org/10.1088/1361-651X/ab580c
• 29)   N. T.  Brown,  E.  Martinez and  J.  Qu: Acta Mater., 129 (2017), 83. http://doi.org/10.1016/j.actamat.2017.02.033
• 30)   S.  Plimpton: J. Comput. Phys., 117 (1995), 1. https://doi.org/10.1006/jcph.1995.1039
• 31)   G.  Bonny,  R. C.  Pasianot,  D.  Terentyev and  L.  Malerba: Philos. Mag., 91 (2011), 1724. https://doi.org/10.1080/14786435.2010.545780
• 32)   C. A.  Becker,  F.  Tavazza,  Z. T.  Trautt and  R. A.  Buarque de Macedo: Curr. Opin. Solid State Mater. Sci., 17 (2013), 277. https://doi.org/10.1016/j.cossms.2013.10.001
• 33)   S.  Nosé: J. Chem. Phys., 81 (1984), 511. https://doi.org/10.1063/1.447334
• 34)   W. G.  Hoover: Phys. Rev. A, 31 (1985), 1695. https://doi.org/10.1103/PhysRevA.31.1695
• 35)   T.  Schneider and  E.  Stoll: Phys. Rev. B, 17 (1978), 1302. https://doi.org/10.1103/PhysRevB.17.1302
• 36)   A.  Stukowski: Model. Simul. Mater. Sci. Eng., 18 (2010), 015012. https://doi.org/10.1088/0965-0393/18/1/015012
• 37)   P. M.  Larsen,  S.  Schmidt and  J.  Schiøtz: Model. Simul. Mater. Sci. Eng., 24 (2016), 055007. https://doi.org/10.1088/0965-0393/24/5/055007
• 38)   A.  Karma: Phys. Rev. E, 48 (1993), 3441. https://doi.org/10.1103/PhysRevE.48.3441
• 39)   W. R.  Fehlner and  S. H.  Vosko: Can. J. Phys., 54 (1976), 2159. https://doi.org/10.1139/p76-256