Molecular Dynamics of Solidification

Abstract

Over many years, mesoscale analysis such as the phase-field method has been the mainstream for numerical simulation of solidification. In contrast, our group has taken the initiative in applying molecular dynamics (MD) simulation to various problems in solidification. In this review, recent advances and contributions of MD simulations for solidification are presented. The primary contribution of MD simulation is the derivation of solid-liquid interfacial properties since it is not easy to measure these properties experimentally with high precision. In addition, recent significant progress in computational environments has dramatically expanded the possibilities of MD simulations for solidification. Now, MD simulations with a scale of billion atoms at the micrometer-scale have become a reality, enabling the exploration of analyses previously dominated by mesoscale methods, such as grain growth and dendrite growth. In particular, the dendrite growth at the micrometer-scale presented here represents the first achievement of directly simulating a typical four-fold symmetrical dendrite structure solely through atomic-scale simulations, to the best of the author’s knowledge. Moreover, new attempts at the fusion of data-driven methods and MD simulations are presented in this review, aiming to contribute to the rapid development of the field of solidification in the future.

1. Introduction

Many practical metallic materials are manufactured using the casting and solidification processes. In order to achieve high functionality and high quality of materials, it is required to precisely control the solidification microstructure in the casting and solidification processes.¹⁾ However, direct observation of nucleation and solidification is generally not easy, although some advanced studies have succeeded in direct observation of solidification process.^2,3) Therefore, numerical simulations have been widely used for the research of solidification.^1,4,5) Mesoscale methods such as phase-field method (PFM),^6,7,8,9,10) cellar automata (CA)¹¹⁾ and dendritic needle network (DNN)¹²⁾ have been mainly used for numerical simulation of solidification. Recently, bottom-up numerical approaches such as molecular dynamics (MD) simulations^4,5,13) have become available for the investigation of solidification due to the significant progress in computational environments.^14,15) The main contributions of MD simulations in solidification are

(i) Derivation of solid-liquid interfacial properties,

(ii) Direct observation of solid-liquid interfacial structure at atomistic scale,

(iii) Direct observation of nucleation process at atomistic scale.

For (i), MD simulations are anticipated to systematically derive interfacial properties, such as solid-liquid interfacial energy,^{16,17,18,19,20,21)} kinetic coefficient^{17,19,20,21,23)} and partition coefficient.^23,24,25) These properties are experimentally challenging to measure with high precision at atomistic scale. It is also expected to play a role in the multiscale modeling by providing these interface properties obtained by the MD simulation as parameters for PFM.²⁶⁾

On the other hand, (ii) and (iii) represent efforts to directly perform MD simulations for phenomena that were previously investigated using mesoscale methods. This shift is attributed to the significant progress in computational environments, allowing the utilization of large-scale MD simulations for such investigations. Regarding (ii), for instance, the solidification process of a system that includes a triple point consisting of grain boundary and solid-liquid interface,²⁷⁾ as well as the anisotropic growth process of solid nuclei,²⁸⁾ can be examined solely through atomic trajectories derived from MD simulations, without the need for introducing phenomenological parameters. For (iii), it is not straightforward to handle nucleation as it is in the PFM by definition. Therefore, it is common to place the initial nuclei with a predetermined probability.^29,30) In the MD simulation, on the other hand, nucleation occurs spontaneously under certain conditions.^31,32) In other words, when the spatiotemporal scale is sufficiently large, it becomes feasible to investigate the formation process of polycrystalline structure through simultaneous nucleation by MD simulations.^14,15)

Considering the background mentioned above, this review focuses on recent advances and contributions of MD simulations for solidification. A brief overview of the MD simulation is first provided. Subsequently, examples of deriving the solid-liquid interfacial energy are presented, which are particularly crucial in the realm of solidification. Following this, a comprehensive discussion of nucleation and solidification from atomistic viewpoint is conducted incorporating numerous research cases. Moreover, the issues related to the spatiotemporal scale of MD simulations is discussed, and attempts to overcome these challenges are presented in the chapter on Challenges and Perspectives.

2. Brief Summary of Molecular Dynamics Simulation

MD simulation is a deterministic method of atomistic simulation that tracks the motion of all atoms in materials according to the equations of motion. Once the initial position and velocity of each atom and the force acting on the atom are given, the position and velocity of all atoms at a given time are uniquely determined by the equations of motion. The fundamental essence of MD simulation is to numerically integrate the equations of motion for a mass system composed of many atoms. This method has already been established as a numerical simulation technique.^33,34) From a practical standpoint, it is crucial to consider the definition of forces acting on each atom and understand how to derive macroscopic physical quantities from the atomic configuration.

An MD simulation that calculates the forces acting on an atom from the electronic structure for each atomic configuration at each time step by the molecular orbital (MO) method or density functional theory (DFT) is called ab initio MD or first-principle MD simulation.^34,35) While ab initio MD can accurately depict interatomic forces, it is computationally demanding. Therefore, in many cases, the force between atoms is expressed through analytical functions called the interatomic potential. An MD simulation using the interatomic potential for the force acting on atoms is called the classical MD simulation. Ab initio MD simulation is suitable for the quantitative analysis of elementary processes involving electron transfer,^35,36) such as bonding/dissociation of molecules^37,38) and redox reactions.^39,40) On the other hand, classical MD simulation is well-suited for phenomena in which the morphology is determined by cooperative interactions among many atoms, such as phase change^{14,15,28,41,42)} and grain growth.^43,44) Therefore, it is crucial to choose the appropriate method based on the intended purpose of the study. It should be noted that even in the ab initio MD simulation, the motion of atoms follows classical motion according to the equations of motion.

Generally, interatomic potentials are classified into two types: pair (two-body) potentials and many-body potentials.¹³⁾ In the pair potential, the total energy U of the system is described only by the distance r_ij between the two atoms i and j as

  

$$\begin{array}{c}U=\frac{1}{2}{\displaystyle \sum }_i{\displaystyle \sum }_j{U}_{ij}(r_{ij})\end{array}$$(1)

Representative two-body potentials include the Lennard-Jones⁴⁵⁾ and the Morse⁴⁶⁾ potentials. While two-body potentials are easy to handle, they are known to have limitations, such as the inability to accurately represent mechanical constants like elastic constants and the behavior of atoms on surfaces or in defects.¹³⁾ Therefore, numerous interatomic potentials that take into account many-body effects have been proposed. For example, bond-order type potentials such as the Tersoff potential,⁴⁷⁾ which takes into account the coordination number and angular term, are the most popular for interatomic potentials for covalent bonds. Furthermore, the ReaxFF potential,⁴⁸⁾ which incorporates charge information and can accurately reproduce chemical reactions, is widely utilized as the interatomic potential in reactive MD simulations. On the other hand, the embedded atom method (EAM) potential⁴⁹⁾ is the most popular for metallic systems. The EAM potential is composed of a two-body term V_ij and an embedding function F, which is a functional of the electron density ρ_i around atom i as

  

$$\begin{array}{c}U=\frac{1}{2}{\displaystyle \sum }_i{\displaystyle \sum }_j{V}_{ij}(r_{ij})+{\displaystyle \sum }_{i}F({\rho }_i)\end{array}$$(2)

  

$$\begin{array}{c}{\rho }_i={\displaystyle \sum }_{j\ne i}{\phi }_{ij}\end{array}(r_{ij})$$(3)

where ρ_i is the electron density around atom i, which is described by the sum of the electron densities φ_ij created at the position of atom i by neighboring atoms j excluding itself. In addition, a highly accurate modified EAM (MEAM) potential,^50,51) which takes into account angle-dependent terms, has also been proposed. To date, numerous EAM potentials have been proposed for pure metals and alloys, with many of them publicly available in databases.^52,53) Moreover, there has been a growing trend in constructing machine learning potentials using neural networks to accurately reproduce results from first-principles calculations. The machine learning-based potential (commonly referred to as neural network potential, NNP^54,55)) accurately represents the multidimensional potential energy surface of the target system by learning reference data obtained from first-principles calculations. To date, numerous NNPs have been proposed, and a recently introduced NNP can handle all combinations of 45 different elements.⁵⁶⁾

Once the interatomic potential is defined, the force acting on the atom is determined as the negative gradient of the potential function at the atomic position, taking the form of a conserved force. By numerically integrating the equations of motion, one can determine the positions and velocities of atoms for a small time step. In recent years, a large amount of numerical data on interatomic potentials has been made available as databases.^52,53) By integrating these databases with the open-source software such as Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS),⁵⁷⁾ significant reductions in the effort and time required for MD simulations can be achieved. Moreover, post-processing and visualization tools such as the Open Visualization Tool (OVITO)⁵⁸⁾ make it easy to perform quantitative analysis as well as intuitive understanding of the simulation results. Among simulation results presented in this review, the results in Sections 3.2, 3.3, 4.4, 4.5, 4.6, 5.2 and 5.3 were obtained using the LAMMPS simulator and the potential data from NIST-Interatomic Potential Repository (IPR),⁵²⁾ while the outcomes of large-scale MD simulations in Section 4.2 and 4.3 were obtained using a code of our own design.

3. Derivation of Solid-liquid Interfacial Properties

3.1. Overview

Thermodynamic and kinetic properties of the solid-liquid interface are known to be important quantities that govern the nucleation rate and morphology of solidification microstructure.^59,60) Therefore, significant efforts have been made over the years towards the high-precision measurement of solid-liquid interface properties. Turnbull^61,62) conducted measurements of solid-liquid interfacial energy through undercooling experiments. This method, based on the homogeneous nucleation theory, involves measuring the solid-liquid interfacial energy from the frequency of crystal nucleus formation and has been reported for various pure metals. Glicksman et al.⁶³⁾ have measured the solid-liquid interfacial energy by the dihedral angle measurement method. In this method, the solid-liquid interfacial energy is measured from the relationship between the solid-liquid interfacial energy and the grain boundary energy at the root of a tilt-boundary groove. It was reported that the original dihedral angle measurement method could not accurately determine the anisotropy of the solid-liquid interfacial energy, but recent improvements in measurement methods have made it possible to measure the solid-liquid interfacial energy with anisotropy taken into account.^64,65) Substantial efforts have continued to be made to experimentally measure solid-liquid interfacial energy. However, due to the challenges in achieving high-precision measurements at high temperatures, there have been limitations in reporting experimental values for solid-liquid interfacial energy.⁶⁶⁾

Therefore, theoretical models and numerical simulations are widely used to estimate the solid-liquid interfacial energy. In the early stages, Skapski⁶⁷⁾ proposed a solid-liquid interfacial energy model based on the contribution of latent heat of fusion. Waseda and Miller⁶⁸⁾ estimated the solid-liquid interfacial energy for a variety of metals based on the Ewing’s model⁶⁹⁾ where the interfacial entropy was evaluated from the radial distribution function of the liquid. Spaepan⁷⁰⁾ proposed that the origin of solid-liquid interfacial energy lies in the decrease of entropy in liquid-phase atoms near the solid-liquid interface using a rigid sphere model. This decrease occurs as a result of the orderly arrangement influenced by adjacent solid-phase atoms. Gránásy et al.⁷¹⁾ focused on the entropy loss of the liquid in contact with the crystal to study the dependence of the solid-liquid interfacial energy on the crystal structure. Regarding numerical approaches, MD simulations have mainly contributed to the estimation of solid-liquid interfacial energy. In early studies, the cleaving technique was frequently used to estimate the solid-liquid interfacial energy.^72,73) Then, the capillary fluctuation method¹⁶⁾ has gained popularity, where the interfacial stiffness is derived from the capillary fluctuation spectrum of the roughness of the solid-liquid interface. The classical nucleation theory (CNT)-based technique^18,19) is also proposed, in which the Gibbs-Thomson effect is considered in determining the threshold temperature that separates the shrinkage and expansion of solid nuclei in the undercooled melt. The capillary fluctuation method has been widely extended to various applications, such as the solid-liquid interfacial energy under temperature gradients.^74,75) Recently, a method for deriving solid-liquid interfacial free energy using metadynamics⁷⁶⁾ has also been proposed.^77,78,79,80) The advantage of this method lies in its capability to obtain the interface energy at non-equilibrium temperatures in which the solid-liquid interface is moving. In the remainder of this chapter, a case study of deriving solid-liquid interfacial energy for iron-based alloys using the capillary fluctuation method are presented in Section 3.2. Then, the temperature dependence of solid-liquid interfacial energy are explored using the metadynamics-based method in Section 3.3.

3.2. Capillary Fluctuation Method

First, a brief explanation of the methodology of the capillary fluctuation method is given by following the original paper¹⁶⁾ and our previous study.²¹⁾ In an equilibrium solid-liquid biphasic system, the solid-liquid interface fluctuates to minimize the free energy. The interface stiffness can be estimated from the amplitude of the fluctuation of the solid-liquid interface based on the following equation:

  

$$⟨{\left|A(k)\right|}^2⟩\begin{array}{c}=\frac{k_B{T}_e}{bW(\gamma +\gamma \text{'}\text{'})k^2}\end{array}$$(4)

where γ is the solid-liquid interfacial energy, A(k) is the Fourier amplitude of the interface height profile, k is the wave number, k_B is Boltzmann’s constant. b and W are the thickness and length of the interface, respectively. T_e is the equilibrium temperature at which the solid-liquid interface remains stationary. The brackets represent averages over the number of sampling times. The term $\gamma +\gamma \text{'}\text{'}$ corresponds to the interfacial stiffness, $\tilde{\gamma }$. From the relationship given in Eq. (4), the interfacial stiffness can be estimated by examining the slope of the relationship between k² and 1/$⟨{\left|A(k)\right|}^2⟩$ in the fluctuation spectrum of a quasi-one-dimensional interface. The interfacial stiffness is expressed by the following equation,⁸¹⁾ separate from the aforementioned relationship, as

  

$$\begin{array}{c}\tilde{\gamma }={\left(\gamma +\frac{d^2\gamma }{d{\theta }^2}\right)}_{\theta =0}\end{array}$$(5)

where θ represents the angle between the instantaneous local normal to the interface and the average orientation for the reference flat interface. The interfacial stiffness for some specific orientations⁸¹⁾ is expressed as

  

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

  

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

  

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

where γ₀ represents the average solid-liquid interfacial energy, and ε₁ and ε₂ represent anisotropy parameters, respectively. It is known that the crystal preferentially grows in the <100> orientation when ε₁ is dominant, whereas it preferentially grows in the <110> orientation when ε₂ is dominant.⁸²⁾ γ₀, ε₁ and ε₂ are specifically obtained by assigning values of interfacial stiffness for three different orientations into Eqs. (6), (7), (8). Finally, the interface energy, including its anisotropy, can be expressed using a low-order expansion consistent with cubic symmetry⁸³⁾ as

  

$$\begin{array}{c}\gamma ={\gamma }_0\left[1+{\epsilon }_1\left({\displaystyle \sum }_{i=1}^3{n}^4-\frac{3}{5}\right)+{\epsilon }_2\left(3{\displaystyle \sum }_{i=1}^3{n}^4+66n_1^2{n}_2^2{n}_3^2-\frac{17}{7}\right)\right]\end{array}$$(9)

where n_i represent the orientation indices of the solid-liquid interface, (n₁, n₂, n₃).

Figure 1 illustrates a case of deriving the solid-liquid interfacial energy of the pure iron using the capillary fluctuation method.²¹⁾ A quasi-two-dimensional biphasic system (Fig. 1(a)) was equilibrated at the melting point under the isobaric-isothermal condition for 100 ps at 0 Pa. The EAM potentials fitted by Bonny et al.^84,85) were employed for pure iron and iron-base alloys (Fe–Cr and Fe–Ni) described later. The polyhedral template matching (PTM)⁸⁶⁾ is used to classify the local atomic configuration. The details of the calculation procedure can be found in Ref. 21). Figure 1(b) depicts an enlarged view of the solid-liquid interface in the biphasic system of pure iron with (100)[010] orientation after 100 ps simulation. The solid-liquid interface was illustrated by the red curve, which can be expressed as a function of the z-direction, h(z). In a similar fashion, the position of the solid-liquid interface was extracted from the atomistic configuration of the MD simulation between 60 and 100 ps at 0.1 ps intervals to gather statistics. Then, fast Fourier transform (FFT) was performed on h(z) to derive the relationship between the wave number and its amplitude. Figure 1(c) shows the relationship between 1/bw$⟨{\left|A(k)\right|}^2⟩$ and k² for the case of the solid-liquid interface of pure iron with (100)[010] orientation. Note that it may lead to undersampling in the range where k is small.^81,87,88) Moreover, the relationship in Eq. (4) does not be satisfied due to the limitation of crystal lattice spacing for large k. Therefore, the solid-liquid interface stiffness was determined within the range of wave numbers k where the slopes of log₁₀$⟨{\left|A(k)\right|}^2⟩$ and log₁₀k² are −1.²¹⁾ From the slope of the fitted line, the interfacial stiffness is estimated to be 132.7 mJ/m². In a similar manner, the interfacial stiffness of biphasic systems with (110)[001] and (110)[110] orientations for pure iron is estimated to be 156.0 and 112.4 mJ/m², respectively. Substituting these values of interfacial stiffness for the three orientations into Eqs. (6), (7) and (8), the average solid-liquid interfacial energy γ₀ was derived as 143.4 mJ/m², and anisotropy parameters ε₁ and ε₂ were derived as 0.037 and −0.0053, respectively. The values of anisotropy parameters are close to the reported values in previous studies.²⁶⁾ Finally, the solid-liquid interfacial energy of (100), (110) and (111) orientations are derived by substituting γ₀, ε₁ and ε₂ into Eq. (9).

Fig. 1. Derivation of the solid-liquid interfacial energy of iron by the capillary fluctuation method after Ueno and Shibuta,²¹⁾ reprinted with permission from Elsevier. (a) Quasi-two-dimensional solid-liquid biphasic system of iron with the (100) interface. The short side of the system is set to be [010] direction. (b) Enlarged view of the solid-liquid interface after 100 ps simulation. (c) The relationship between 1/bw$⟨{\left|A(k)\right|}^2⟩$ and k², obtained from the fast Fourier transform for the solid-liquid interface. (Online version in color.)

Figure 2 illustrates the results depicting the composition dependence of the interfacial energy for Fe–Cr and Fe–Ni alloys using the aforementioned method.^21,89) The average values of five calculations (Fe–Cr) or three calculations (Fe–Ni) are plotted for each composition. The solid-liquid interfacial energy of the Fe–Cr alloy decreases with an increasing Cr composition, reaching a minimum at approximately 19.1 at%Cr in Fe, and then starts to increase. Regarding the Fe–Ni alloy, a face-centered cubic (fcc) structure naturally transforms into a body-centered cubic (bcc) one during relaxation for 15 at%Ni or less. Therefore, a threshold composition dividing the fcc and bcc phases is defined between 15 and 30 at%Ni, which is higher than the value reported in experiments (around 5 at%Ni) due to the characteristics of the EAM potential.⁸⁹⁾ The solid-liquid interfacial energy for Ni-rich compositions with an fcc solid structure is larger than that for Fe-rich compositions with a bcc solid structure. In the fcc solid region, the solid-liquid interfacial energy generally increases with the rising Ni composition until around 50 at%Ni, after which it becomes almost constant with some fluctuation beyond 50 at%Ni. On the other hand, in the bcc solid region, the solid-liquid interfacial energy typically decreases with the increasing Ni composition. The dependence of the interfacial energy on orientation is much smaller than the composition dependence for both Fe–Cr and Fe–Ni alloys. However, there is a clear trend with γ₍₁₀₀₎>γ₍₁₁₀₎>γ₍₁₁₁₎ for both alloys at all compositions. This trend agrees with previous reports.^90,91) It should be noted that the absolute value of the solid-liquid interfacial energy strongly depends on the interatomic potential used. However, these results from MD simulations are very useful for understanding the qualitative trend of the composition dependence of the solid-liquid interfacial energy, as it is challenging to determine experimentally.

Fig. 2. Composition dependence of the solid-liquid interfacial energy including the anisotropy effect for (a) Fe–Cr and (b) Fe–Ni alloys after Ueno and Shibuta,^21,88) reprinted with permission from Elsevier. The horizontal axis indicates mole fraction. Insert figures represent enlarged views of representative plots. (Online version in color.)

3.3. Metadynamics-based Method

As described above, the capillary fluctuation method has contributed to the estimation of solid-liquid interfacial energy of many metals and alloys. On the other hand, these estimations are basically confined to the solid-liquid interfacial energy at the equilibrium temperature, where the solid-liquid interface remains stationary. To overcome this problem, a metadynamics-based method^77,78,79,80) is recently proposed, in which the solid-liquid interfacial energy is estimated directly from a free energy profile of the solid-liquid biphasic system. Metadynamics⁷⁶⁾ is an accelerated sampling method used to calculate the free energy profile along the reaction coordinate, also referred to as the collective variable (CV). Metadynamics is widely employed for deriving free energy surfaces of various reaction systems^92,93) and accelerating MD simulations,^94,95,96) in conjunction with hyperdynamics.⁹⁷⁾ Here, an enhanced metadynamics-based method that we proposed for deriving solid-liquid interfacial energy is presented following our previous paper.⁸⁰⁾ This method builds upon the original approach.^77,78,79)

The free energy difference between a solid-liquid biphasic system with n solid atoms (G_n) and a reference liquid system with the same number of atoms (G_liquid) can be expressed as:

  

$$\begin{array}{c}{G}_n-G_{\text{liquid}}=n{\mu }_{\text{sl}}+A\gamma \end{array}$$(10)

where μ_sl is the difference in chemical potential between solid and liquid phases, A is the area of the solid-liquid interface, and γ is the solid-liquid interfacial energy. When two systems, each consisting of the same number of atoms but with different areas of the solid-liquid interface (A and 2A), are prepared, the free energy differences for the two systems are expressed as:

  

$$\begin{array}{c}\{\begin{array}{c}{G}_n^1-G_{\text{liquid}}^1=n{\mu }_{\text{sl}}+A\gamma \\ G_n^2-G_{\text{liquid}}^2=n{\mu }_{\text{sl}}+2A\gamma \end{array}\end{array}$$(11)

where superscripts 1 and 2 represent systems with solid-liquid interface areas of A and 2A, respectively. Because both systems consist of the same number of atoms, the probability of the liquid state appearing must be the same in both systems during the metadynamics simulation. This satisfies the relationship $G_{\text{liquid}}^1=G_{\text{liquid}}^2$. Therefore, Eq. (11) can be summarized as:

  

$$\begin{array}{c}{G}_n^2-G_n^1=A\gamma \end{array}$$(12)

Therefore, the solid-liquid interfacial energy can be estimated if the difference in the free energy, $G_n^2-G_n^1$, can be obtained from the metadynamics simulation.

In metadynamics, a bias potential in the following form of the Gaussian function as a function of reaction coordinate q is added to the free energy surface:

  

$$\begin{array}{c}{V}_{\text{bias}}(q,t)={\sum }_{k\tau <t}W\text{(}k\tau \text{)exp}\left[-\frac{{(q(t)-q(k\tau ))}^2}{2{\sigma }^2}\right]\end{array}$$(13)

where σ and W are the width and height of the Gaussian function, k is an integer and τ is the interval at which bias is added (referred to as the waiting time), respectively. The height is adjusted according to the well-tempered metadynamics scheme⁹⁸⁾ as follows:

  

$$\begin{array}{c}W(k\tau )=W_0\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{V_{bias}(q(k\tau ),k\tau )}{k_B\mathrm{\Delta }T}\right)\end{array}$$(14)

  

$$\begin{array}{c}\mathrm{\Delta }T=(\eta -1)T\end{array}$$(15)

where W₀ is the initial height of the bias and η is a bias factor.

Here, a study of deriving solid-liquid interfacial energy of iron using the metadynamics-based method is presented in line with our previous report.⁸⁰⁾ The number of solid atoms n was employed as the reaction coordinate. Bcc crystals of iron, comprising 18 × 8 × 8 unit cells (2304 atoms) and 36 × 8 × 4 unit cells (2304 atoms), were annealed at 1718 K. The (100) plane is oriented perpendicular to the x-axis (the long axis direction). The pre-annealing was performed for 25 ps using the isothermal-isobaric ensemble to achieve an equilibrium lattice size at the target temperature before the main calculation. The area of the (100) plane in the latter system (2A) is twice that of the former one (A). Then, the prepared systems were kept isothermal for 24975 ps calculation in conjunction with metadynamics using the isothermal-isobaric ensemble, with only the x-direction left free to fluctuate (i.e., the areas of the solid-liquid interface remain constant throughout the calculation). Parameters for the metadynamics simulation were set to be W₀ = 0.01 eV, σ = 600 and η = 90. To perform metadynamics calculations on LAMMPS, PLUMED plugin⁹⁸⁾ was employed. The EAM potentials fitted by Bonny et al.⁸⁵⁾ were employed for the interaction between iron atoms. The details of the calculation procedure can be found in Ref. 80).

Figure 3(a) illustrates the free energy profiles as a function of the number of solid atoms at 1718 K for two iron systems with distinct solid-liquid interface areas.⁸⁰⁾ The free energy profile exhibits two wells around the minimum and maximum values of n, representing liquid and solid states. Additionally, there is a plateau region between the two states, indicating the coexistence of solid and liquid phases. Figure 3(b) presents representative snapshots of the calculation system at the corresponding states indicated by arrows in Fig. 3(a). The snapshots confirm that solid, liquid, and solid-liquid coexisting states are accurately reproduced in the biased simulation with metadynamics. The free energy of the liquid state is higher than that of the solid state, and the free energy decreases with an increase in the number of solid atoms in the plateau region, which accurately reflects that 1718 K is below the melting point of iron. The free energy in the plateau region for the system with a solid-liquid interface area of 2A is higher than that for the system with an area of A, as the presence of the solid-liquid interface increases the excess free energy. The solid-liquid interfacial energy can be determined from the difference in free energy at the plateau region between systems with areas A and 2A, as per Eq. (12). Utilizing the average values of the free energy difference between n = 1000 and 1200, the solid-liquid interfacial energy of iron at 1718 K is estimated to be 220.8 mJ/m².⁸⁰⁾

Fig. 3. Derivation of the solid-liquid interfacial energy of iron by the metadynamics-based method.⁸⁰⁾ (a) Free energy profile with respect to the number of solid atoms, n for iron system at 1718 K. Blue and orange lines represent the free energy of the system with a solid-liquid interface area of A and 2A, respectively. (b) Representative snapshots of the calculation systems at the corresponding states indicated by arrows in (a). (Online version in color.)

In the same way, the solid-liquid interracial energy at various temperatures between 1668 and 1888 K at intervals of 50 K was estimated.⁸⁰⁾ Calculations were performed five times for each temperature to gather statistics. Figure 4 depicts the temperature dependence of the solid-liquid interfacial energy of bcc-iron. The plots represent the average value of five calculations, with error bars indicating the standard deviation. The solid-liquid interfacial energy takes values approximately between 220 and 270 mJ/m², which falls within the range of reported values for the solid-liquid interfacial energy at the melting point.^{61,62,68,99,100,101)} Simulation results suggest that the solid-liquid interfacial energy slightly decreases with increasing temperature, a trend that is particularly evident at temperatures below the melting point. On the other hand, the temperature dependence of the solid-liquid interfacial energy near the melting point is still a subject of ongoing debate. While some literatures report an increase in the solid-liquid interfacial energy with increasing temperature (i.e., the positive temperature dependence^71,102,103)), other sources suggest a decrease in the solid-liquid interfacial energy with increasing temperature (i.e., the negative temperature dependence^104,105,106)). Therefore, further research is essential to understand the temperature dependence of the solid-liquid interfacial energy.

Fig. 4. Temperature dependence of the solid-liquid interfacial energy of bcc-iron.⁸⁰⁾ The error bar represents the standard deviations of five calculations. Several reported values of the solid-liquid interface of iron at the melting point from experiments (EXP)^{61,62,100,101)} and theoretical calculations (CALC)^68,99) are plotted for comparison. (Online version in color.)

4. Nucleation and Solidification from Atomistic Viewpoint

4.1. Overview

Since many metallic materials are produced through solidification that commence with nucleation, it is essential to comprehend and control nucleation in the materials processing.⁶⁰⁾ Although the PFM has been mainly used for the numerical simulation of solidification, it is not straightforward to handle nucleation as it is in the PFM by definition. Therefore, MD simulations have been widely employed to elucidate the nature of nucleation from an atomistic perspective.¹⁰⁷⁾ For example, estimation of nucleation rate and barrier,^108,109,110) non-classical behavior of nucleation,^111,112) transient clusters during nucleation,^113,114,115) and polymorphism,^116,117,118) have been thoroughly investigated through MD simulations. In the early days, nucleation simulations were often restricted to investigating one or a few nuclei due to limitations in computational performance. However, recent dramatic progress in computational capabilities has enabled the simulation of simultaneous nucleation process, allowing for the formation of polycrystalline structures.^{14,15,119,120,121,122,123,124)} For further details, refer to the comprehensive review of theoretical and computational studies on nucleation in Ref. 107).

4.2. Homogeneous Nucleation and Solidification

Here, a study of MD simulation of simultaneous homogeneous nucleation from an undercooled iron melt, leading to the formation of polycrystalline structures, is presented in line with our previous reports.^4,119) Nucleation and subsequent solidification were induced through isothermal holding of a simulation cell filled with an iron melt at various undercooling temperatures. A quasi-two-dimensional system of iron measuring 53.4 × 53.4 × 4.3 nm (186 × 186 × 15 unit cells, 1037880 atoms) was employed. The Finnis–Sinclair potential¹²⁵⁾ was used for the interatomic potential between iron atoms. The periodic boundary condition was applied to all boundaries. Simulations were conducted on the GPU supercomputer, TSUBAME2.5, at the Tokyo Institute of Technology, using the original code developed with CUDA 3.0.¹²⁶⁾ Further details of the calculation procedure can be found in Ref. 119). Figure 5 depicts snapshots of the atomic configuration during a consecutive simulation of nucleation, solidification, and grain growth at 0.58 T_m for 10000 ps, where T_m represents the melting point of iron according to the Finnis-Sinclair potential. Numerous nuclei nucleate simultaneously before 150 ps and grow into spherical grains within the melt. Meanwhile, additional nuclei continuously form from the remaining melt. Grains nucleated later fill the spaces between those nucleated earlier, and the entire iron melt solidifies by 300 ps. Following solidification, the smaller grains gradually shrink and disappear, while the larger ones continue to grow.

Fig. 5. Consecutive molecular dynamics simulation of nucleation, solidification and grain growth.⁴⁾ An iron melt in a quasi-two-dimensional cell measuring 53.4 × 53.4 × 4.3 nm³ was isothermally held at 0.58 T_m for 10000 ps, where T_m represents the melting point of iron according to the Finnis-Sinclair potential. Purple and white atoms correspond to iron atoms with and without the bcc configuration, respectively. (Online version in color.)

Similar simulations were conducted at various temperatures ranging from 0.40 T_m to 0.67 T_m, and the number of nuclei in the calculation cell during the first 1000 ps was counted. The calculations were performed five times for each temperature to gather statistics. Figure 6 displays the nucleation rate as a function of temperature, obtained from five replicate calculations at each temperature. The curve connecting the plots exhibits a distinctive shape, featuring a peak at the critical temperature (i.e., the nose shape). The distinctive nose shape is indicative of the characteristics of thermally-activated nucleation.⁶⁰⁾ Specifically, the nucleation rate exhibits a rapid increase with decreasing temperature, corresponding to a reduction in the nucleation energy barrier. On the other hand, the mobility of atoms decreases with lowering temperature due to the thermally-activated nature of atomic mobility processes. The interplay between these competing effects manifests in the observed nose shape when plotting the nucleation rate against temperature. it should be stressed that thermally activated nucleation occurs naturally without any inducing factors in the MD simulation, whereas PFM generally requires predefined rules for nucleation.

Fig. 6. Nucleation rate as a function of temperature derived from the molecular dynamics simulation of homogeneous nucleation from undercooled melt iron.¹¹⁹⁾ The average number obtained from five replicate calculations for each temperature is plotted with error bars showing the standard deviation.

4.3. Heterogeneity in Homogeneous Nucleation

The morphology of microstructure obtained through nucleation is significantly dependent on the probability of nucleation. Therefore, it is essential to perform a large-scale simulation for a statistical discussion of the orientation relationship and particle size distribution of grains. Therefore, we have implemented our MD code on multiple graphics processing units (GPUs) and conducted over 500 GPU parallel computations for billion-atom MD simulations of nucleation.^14,15) Figure 7 shows snapshots of the simulation cell (236.8 × 236.8 × 236.8 nm³) with 1024000000 atoms during a consecutive simulation of nucleation and solidification from undercooled iron melt at 0.58 T_m. The calculation method is basically the same as in the previous section. Further details of the calculation procedure including the implementation on the GPU supercomputer can be found in Ref. 14). At 100 ps, numerous small grains appear in the simulation cell and subsequently grow larger over time. Additional nucleation also takes place in the remaining undercooled melt. By 200 ps, solidification is nearly complete, resulting in the filling of nearly all spaces with numerous grains exhibiting various orientations, forming a fine microstructure. Subsequently, the small grains observed in the 200 ps snapshot at the atomistic scale disappear by 2000 ps.

Fig. 7. Consecutive molecular dynamics simulation of nucleation and solidification performed on the GPU-rich supercomputer.¹⁴⁾ An iron melt consisting of 1024000000 atoms in a cell measuring 236.8 × 236.8 × 236.8 nm³ was isothermally held at 0.58 T_m for 2000 ps, where T_m represents the melting point of iron according to the Finnis-Sinclair potential. (Online version in color.)

The distribution of the grain orientation in the microstructure is investigated in detail. Figure 8(a) illustrates the distributions of disorientation angles between neighboring grains at representative time points. The dashed line represents the Mackenzie distribution function¹²⁷⁾ reflecting a random orientation. The distribution of disorientation angles between neighboring grains approximately conforms to the Mackenzie distribution, indicating that the grains are essentially randomly oriented in space. However, a distinct peak is observed at approximately 60°. Grain boundaries with Σ3 misorientations have a disorientation angle of 60°.¹²⁸⁾ Among these boundaries, there is a coherent twin boundary with an exceptionally low grain boundary energy, specifically the bcc(112) Σ3 plane for the <110> tilt axis with a tilt angle of 109.47°.¹²⁹⁾ Figure 8(b) illustrates snapshots of the atomic configuration of grains with twin boundaries extracted from the simulation cell. It is confirmed that grain boundaries with a tilt angle of 109.47° exist along with a <110> tilt axis between the two grains in contact, which compose a coherent twin boundary. Given that a coherent twin boundary in bcc iron possesses an extremely low grain boundary energy,¹²⁹⁾ it is anticipated that grains with subgrains containing a coherent twin boundary are readily formed during the solidification process. That is, heterogeneous nucleation occurs on the surface of the previously existing grain. The heterogeneity in homogeneous nucleation is a phenomenon that cannot be explained by the classical nucleation theory. Therefore, using MD simulation is advantageous for uncovering a detailed understanding of nucleation on an atomistic scale.

Fig. 8. Orientation relationship for all grains during the molecular dynamics simulation of nucleation in Fig. 7.¹⁴⁾ (a) Distributions of disorientation angle between neighboring grains. The dashed line is the Mackenzie distribution function, representing a random orientation. (b) Snapshots of atomic configuration of grains with twin boundaries extracted from the simulation cell. (Online version in color.)

4.4. Polymorphism

The nucleation of a metastable phase under a large undercooling temperature, known as Ostwald’s step rule,¹³⁰⁾ has been long recognized and has been widely observed experimentally over many years.^131,132,133) Following Ostwald’s step rule, the free energy hierarchy represents the relative magnitude of the driving force between possible phase transitions, leading to a stepwise phase transition through a metastable phase.^132,133) The stepwise phase transition has been extensively discussed from a thermodynamic viewpoint.^132,133) Moreover, polymorphism during nucleation has been studied in MD simulations.^116,117,134) in which bcc-like nucleus first appeared and transformed into fcc crystals from inner of the preordered cluster even the fcc phase should be stable in these systems.

Here, a study of MD simulation of polymorphism during homogeneous nucleation from undercooled nickel melt is presented in line with our previous report.¹¹⁸⁾ A nickel melt consisting of 4000000 atoms in a cell measuring 35.2 × 35.2 × 35.2 nm³ was isothermally held under the isobaric-isothermal condition at 1100 K for 3000 ps, The EAM potentials fitted by Purja Pun and Mishin¹³⁵⁾ were employed for the interaction between nickel atoms. Further details of the calculation procedure can be found in Ref. 118). Figure 9(a) depicts snapshots of stepwise phase transition from bcc to fcc phase within a nucleus during the nucleation from undercooled nickel melt. A bcc nucleus emerges from the undercooled melt and transforms into an fcc nucleus from within the cluster, aligning with previous findings from MD simulations.^116,117,134)

Fig. 9. Molecular dynamics study of polymorphism in nucleation from undercooled nickel melt.¹¹⁸⁾ (a) Snapshots of stepwise phase transition from fcc to bcc within a nucleus. Scale bars represent 5 Å. Blue, green and red atoms represent atoms with bcc, fcc and hcp configurations. (b) Difference in free energy as a function of nucleus radius for bcc nucleation from an undercooled melt (blue curve), fcc nucleation from an undercooled melt (green curve), and fcc nucleation from a bcc crystal (black curve). (Online version in color.)

Following the classical nucleation theory,⁶⁰⁾ the nucleation of a spherical nucleus with a radius, r in the undercooled melt results in a difference in free energy, ΔG as follows:

  

$$\begin{array}{c}\mathrm{\Delta }G=4\pi r^2\gamma -\frac{4}{3}\pi r^3\mathrm{\Delta }{G}_V\end{array}$$(16)

where ΔG_V represents the difference in Gibbs free energy between liquid and solid. The first term indicates energy loss due to the appearance of the solid-liquid interface, while the second term indicates energy gain due to the increase in volume of the solid phase. If the nucleus is too small, the energy gain from volume increase cannot surpass the energy loss due to surface energy, and nucleation is not promoted. When the nucleus size exceeds the critical radius, r*, at which ΔG is maximized in Eq. (16), nucleation is facilitated. In the case of the nickel system, ΔG_V between fcc and liquid phases becomes larger than that between bcc and liquid phases since fcc is the most stable phase for solid nickel. On the other hand, the solid-liquid interfacial energy of bcc-liquid interface is lower than that of fcc-liquid interface in general.^118,136) Therefore, it is possible that bcc nuclei may first appear in the undercooled nickel melt in the early stages of nucleation since the effect of the interfacial energy term is dominant. Figure 9(b) depicts ΔG as a function of nucleus radius for following three cases: (1) bcc nucleation from undercooled melt at 1100 K, (2) fcc nucleation from undercooled melt at 1100 K, and (3) fcc nucleation from bcc crystal, as estimated from Eq. (16). Parameters derived from a separate MD calculation were used for the derivation. Further details of the parameter derivation can be found in Ref. 118). Comparing the first and second cases, ΔG for bcc nucleation from undercooled melt is lower than that for fcc nucleation for nuclei with a radius smaller than approximately 4 Å. Therefore, it is reasonable to expect that a bcc nucleus first appears in the undercooled melt in the MD simulation. However, the difference in ΔG between bcc and fcc nucleation from the liquid is not so significant, which does not rule out the possibility of direct fcc nucleation from the undercooled melt. Moreover, ΔG for fcc nucleation from a bcc matrix becomes lower than that for bcc nucleation from undercooled melt at radii larger than approximately 9 Å, where these two lines intersect. In other words, fcc nucleation can take place from the interior of a pre-existing bcc nucleus with a radius of approximately 10 Å or larger. This implies that the stepwise phase transition can occur through polymorphism within the nucleus rather than mere competition between bcc and fcc nucleation from an undercooled melt. It should be noted that the origin of polymorphism during nucleation is still a subject of ongoing debate and the discussion presented here is one possible explanation based on MD simulations.¹¹⁸⁾

4.5. Heterogeneous Nucleation

MD simulations introduced so far have been focused on homogeneous nucleation. On the other hand, the actual solidification process is dominated by heterogeneous nucleation from impurities and walls. According to the classical nucleation theory,⁶⁰⁾ the change in free energy due to heterogeneous nucleation of a solid phase in a spherical-cap shape with a radius r on the substrate with a contact angle, θ is given as

  

$$\begin{array}{c}\mathrm{\Delta }{G}_{hetero}=\left(4\pi r^2\gamma -\frac{4}{3}\pi r^3\mathrm{\Delta }{G}_V\right)\times \frac{2-3cos\theta +{cos}^3\theta }{4}\end{array}$$(17)

The contact angle is determined by the balance of the solid-liquid interfacial energy, the surface energy of the substrate in the solid phase, and the surface energy of the liquid phase at the triple point. However, as seen in Chapter 3, it is challenging to experimentally determine interface properties with high precision. That is, a quantitative discussion on heterogeneous nucleation using experimental values is not straightforward. Therefore, MD simulation is again a powerful tool for investigating heterogeneous nucleation from an atomistic viewpoint. For example, the lattice misfit between the substrate and solidified metals^137,138,139) and the effect of grain refiners on solidification^31,32) are widely investigated by MD simulations. For further details, refer to an overview of atomistic mechanisms of heterogeneous nucleation in Ref. 139).

In practical production processes of metallic materials, grain refiners are frequently inoculated into undercooled melt to promote and control nucleation.^140,141) Although the narrow definition of nucleation is a stochastic process driven by thermal activation, it is considered that the onset of heterogeneous nucleation on inoculated fine particles is not stochastic but deterministic process dependent on undercooling temperature and particle size.¹⁴²⁾ This type of nucleation is called the free growth¹⁴¹⁾ or athermal nucleation.¹⁴²⁾ Here, a study of MD simulations of athermal heterogeneous nucleation via grain refiners is presented in line with our previous report.³¹⁾ A liquid aluminum was prepared by annealing an fcc crystal of aluminum consisting of 80 × 80 × 80 unit cells (2048000 atoms) at 2000 K for 10 ps with canonical ensemble. Then, a disk-shaped titanium grain refiner with a diameter of 8 nm was inserted in the center of aluminum melt, omitting liquid atoms located within 2.5 Å from titanium atoms. Energy minimization was carried out on the combined structure to prevent unintended proximity between liquid aluminum and solid titanium atoms at the interface. The prepared system was isothermally held at 830 or 810 K for 600 ps using the isothermal-isobaric ensemble, respectively. Further details of the calculation procedure can be found in Ref. 31).

Figure 10(a) depicts snapshots during the heterogeneous nucleation of undercooled aluminum melt on the surface of a disk-shaped titanium grain refiner at 830 K and 810 K, respectively. In the initial stage, a thin solid layer forms on the top and bottom surfaces of the grain refiner, and it develops into a spherical cap consisting of fcc-aluminum at both temperatures. At 830 K, the growth of the spherical cap stagnates around 200 ps, and the cap size remains constant thereafter. Conversely, the spherical cap continues to grow during the 600 ps calculation at 810 K, which is regarded as the free growth. Hence, there are two growth modes depending on temperature, with a threshold temperature between 810 K and 830 K. Similarly, changes in temperature and particle size were systematically applied to investigate the growth behavior of solid aluminum on the surface of cubic titanium grain refiners of various sizes. Figure 10(b) illustrates the phase diagram of growth modes, including stagnation, transition, or free growth, in relation to particle size and temperature. All examined sizes exhibit two growth modes, and the threshold temperature increases, i.e., the threshold undercooling temperature decreases with increasing particle size. The threshold temperature can be approximately estimated by determining the midpoint between the highest temperature of free growth and the lowest temperature of stagnation for each particle size. It should be stressed that the temperature range of free growth at each particle size can be defined deterministically based on the phase diagram obtained from the MD simulations.

Fig. 10. Molecular dynamics study of athermal heterogeneous nucleation.³¹⁾ (a) Snapshots of atomic configuration during heterogeneous nucleation of undercooled Al melt on the surface of a disk-shaped Ti grain refiner at 830 K and 810 K, respectively. Green and red atoms represent Al atoms with fcc structure and Ti atoms with hcp structure, respectively. All other atoms are not shown for the clarity of the image. (b) Phase diagram of growth modes, including stagnation (×), transition (Δ), or free growth (○), in relation to particle size and temperature. (Online version in color.)

The threshold undercooling temperature appears to be inversely proportional to the particle size from the results of MD simulations in Fig. 10(b). Here, the physical meaning of the threshold curve in the figure is discussed. Figure 11 depicts schematic images illustrating nucleation on the surface of the grain refiner under conditions of stagnation and free growth. In accordance with the theoretical model of athermal nucleation,¹⁴²⁾ a thin layer of solid covers the surface of a grain refiner in the initial stage of nucleation in the inoculated melt. Subsequently, the thin layer transforms into a spherical cap, growing outward, as represented by its height, h. The radius of the spherical cap decreases as the height increases, reaching its minimum when the solid forms a hemisphere, with the radius corresponding to the particle radius, r_N. Therefore, if the critical radius of nucleation, r* is larger than the particle radius, r_N, growth of the spherical cap on the particle stops when the radius of spherical decreases to equal r*. The additional increase in solid at this stage decreases the radius of spherical cap, leading to the shrinkage of the spherical cap. Consequently, the growth of the spherical cap stagnates when the radius of the spherical cap reaches at r*. On the other hand, the spherical cap can grow freely when r* is smaller than r_N, as the radius of the spherical cap is always larger than r*, which is referred to as the free growth. Following the classical nucleation theory, r* is defined as the radius at which the derivative of ΔG by r is zero in Eq. (16) as

  

$$\begin{array}{c}{r}^{*}=\frac{2\gamma }{\mathrm{\Delta }{G}_V}=\frac{2\gamma T_m}{\mathrm{\Delta }H}\frac{1}{\mathrm{\Delta }T}\end{array}$$(18)

Here, the following relationship is utilized: ΔG_V=∆TΔH/T_m,^19,60) where ΔH is the latent heat and ΔT is the undercooling temperature. Substituting r* for r_N in Eq. (16), the threshold undercooling temperature, ΔT_fg, dividing the stagnation and the free growth, is defined as:

  

$$\begin{array}{c}\mathrm{\Delta }{T}_{\text{fg}}=\frac{2\gamma T_m}{\mathrm{\Delta }H}\frac{1}{r_N}\end{array}$$(19)

The threshold undercooling temperature is inversely proportional to r_N, consistent with the results of the MD simulation observed in Fig. 10(b). Therefore, it is concluded that the athermal nucleation appears spontaneously in the MD simulation under the conditions of free growth.

Fig. 11. Schematic images depicting nucleation on the surface of the grain refiner under conditions of (a) stagnation (r* > r_N) and (b) free growth (r* < r_N), where r* is the critical radius of nucleation, and r_N is the particle radius of the grain refiner.³¹⁾ (Online version in color.)

5. Challenges and Perspectives

5.1. Remaining Issues

The issue in MD simulations of solidification lies in the mismatch of size and time scales compared to typical experimental observation systems. Of course, this is not limited to the field of solidification but rather a general challenge in MD simulations. In general, expanding the time scale of MD simulations is not straightforward, as the temporal sequence of atomic motion cannot be parallelized so easily. Therefore, various techniques for time acceleration in atomistic simulations have been developed, including umbrella sampling,¹⁴³⁾ multicanonical algorithm,¹⁴⁴⁾ hyperdynamics,⁹⁷⁾ temperature accelerated dynamics,¹⁴⁵⁾ metadynamics,⁷⁶⁾ collective variable-driven hyperdynamics (CVHD),⁹⁴⁾ parallel replica dynamics,¹⁴⁶⁾ diffusive molecular dynamics (DMD),¹⁴⁷⁾ and adaptive-boost MD.¹⁴⁸⁾ However, most of these techniques only accelerate infrequent events that do not occur in normal MD simulations, and the remaining parts of the system are more often not accelerated. To address this issue, time acceleration methods for MD simulations based on machine learning approaches have recently been developed.^149,150,151) However, from a practical perspective, establishing them is still considered challenging. The application of machine learning approaches to MD simulation will be discussed later.

On the other hand, recent advancements in computational environments have enabled the achievement of spatial scales comparable to mesoscale methods through massive parallel computing using domain decomposition technique.^5,14,15,152) However, the size that can be handled even in very large-scale MD simulations of more than 10 billion atoms is on the order of micrometers,¹⁵⁾ which is much smaller than the scale of typical solidification microstructure. One idea to achieve a larger spatial scale is to use a coarse-grained method, in which multiple atoms are treated as large particles.^153,154) This method is called the coarse-grained (CG) MD simulations. The CG-MD simulation has been widely applied for studies of proteins,¹⁵⁵⁾ polymers,^156,157) deoxyribonucleic acid (DNA),^158,159) and peptide folding,¹⁶⁰⁾ where a monomer or its aggregate is mainly employed as the unit of CG beads. On the other hand, applying CG method to a metallic system is not straightforward, as metallic systems are not bond-oriented, and the units of CG are not well-defined. Meanwhile, Dongare proposed a quasi-coarse-grained dynamics (QCGD) model¹⁶¹⁾ to apply CG simulations to metallic materials. In this approach, scaling relationships for atomic-scale interatomic potentials were employed to characterize the interactions between these representative atoms, and the equations of motion for representative atoms were solved. The QCGD simulations properly derive various mechanical and thermal properties of metallic systems.^161,162) Taking guidance from this model, we have performed a large-scale CG-MD simulation of the crystal growth on the supercomputer, in which a submicron-order crystal appears from the seed crystal in the undercooled melt.¹⁶³⁾ In the next section, a new study involving a micrometer-scale CG-MD simulation of dendritic growth from a seed crystal is presented, which is a new challenge aiming for the direct observation of dendritic growth solely from atomistic information.

5.2. Direct Observation of Dendritic Growth from Atomistic Simulation

First, an outline of the CG-MD method is given in accordance with our previous report.¹⁶³⁾ A CG model for the fcc crystal of metallic material consists of a new fcc crystal with a lattice constant different from the physical one.^{161,162,163,164)} By selecting the lattice constant of the CG model as n times the lattice constant of the atoms, each CG bead represents n³ atoms and the mass of each CG bead can be expressed as n³ times the atomic mass. When applying the CG model to the EAM potential (Eqs. (2) and (3)), each term is adjusted to ensure that the potential energy attains the same value as:

  

$$\begin{array}{c}{V}_{ij}^{\text{CG}}(r_{ij})=n^3{V}_{ij}^{\text{MD}}\left(\frac{r_{ij}}{n}\right)\end{array}$$(20)

  

$$\begin{array}{c}{F}^{\text{CG}}({\rho }_i)=n^3{F}^{\text{MD}}({\rho }_i)\end{array}$$(21)

  

$$\begin{array}{c}{\phi }_{ij}^{\text{CG}}(r_{ij})={\phi }_{ij}^{\text{MD}}\left(\frac{r_{ij}}{n}\right)\end{array}$$(22)

where the superscripts CG and MD indicate the terms specific to the EAM potential for the CG simulation and the original EAM potential for MD simulations, respectively. Regarding the kinetic energy, $\frac{1}{2}{k}_{B}T$ per degree of freedom is distributed according to the equipartition theorem of energy. That is, 3n³ degrees of freedom of n³ atoms reduce to 3 degrees of freedom for each CG bead when mapping atoms into CG beads. The loss of degrees of freedom resulting from coarse-graining can be considered by multiplying the temperature by n³. Therefore, an equivalent temperature, $T\text{'}=\frac{1}{n^3}T$ is defined as the temperature of a system with CG beads. Figure 12 illustrates the total energies of fcc nickel crystals, measuring 126.72 × 126.72 × 126.72 ų, composed of varying sizes of CG beads (n = 2, 3, and 4) during a 100 ps relaxation at 500 K. The simulations were conducted using the canonical ensemble and the EAM potential developed by Purja Pun and Mishin.¹³⁵⁾ Periodic boundary conditions were applied in all three directions. The average total energy was consistent across all systems, confirming that the beads in the CG model effectively coarse-grained the atoms in the MD simulation for the bulk crystalline system from an energetic perspective. The calculation performance of CG models (n = 2, 3, and 4) was approximately 8, 24, and 55 times faster, respectively compared to the MD simulation, when using a single core of the Intel(R) Xeon(R) CPU E5-2620 v4 (2.10 GHz).¹⁶³⁾ On the other hand, the simulation resolution becomes coarser as the value of n increases. Moreover, reducing the number of beads in the system results in larger thermal fluctuations due to a decreased degree of freedom. To balance calculation performance and simulation resolution, the optimal n must be chosen based on the purpose of the simulation.

Fig. 12. Comparison of the total energy of fcc nickel systems with dimensions of 126.72 × 126.72 × 126.72 ų, comprising different sizes of coarse-grained (CG) beads, during a 100 ps relaxation at 500 K using the canonical ensemble.¹⁶³⁾ (a) 186624 atoms for molecular dynamics (MD), (b) 23328 beads for CG (n = 2), (c) 6912 beads for CG (n = 3) and (d) 2916 beads for CG (n = 4) simulations, respectively. Insets represent snapshots of each calculation system at 100 ps. (Online version in color.)

Using this CG model, a large-scale CG simulation of crystal growth from a seed crystal in undercooled nickel melt in a quasi-two-dimensional system is performed. Here, CG beads with n = 3 were employed to efficiently handle a large system on the order of micrometers. The liquid structure was prepared by heating an fcc crystal consisting of consisting of 1000 × 1000 × 4 unit cells of CG (n=3) (16000000 CG beads) at 3000 K with the canonical ensemble for 10 ps. A 4 nm radius solid circular nucleus was placed at the center of the liquid as a seed crystal, and all liquid atoms within 3.0 Å of a solid atom were removed. The prepared system consists of 16000553 CG beads, which is equivalent to 432014931 nickel atoms. The composite structure was maintained at 1380 K and 0 Pa using the isothermal-isobaric ensemble for a duration of 16000 ps, employing a time step of 1 fs (i.e., 16000000 steps). The Nosé-Hoover thermostat^165,166) was employed. The simulation was performed by 10752-core parallel computation on the supercomputer, Oakbridge-CX at the University of Tokyo. Figure 13 shows snapshots of micrometer-scale CG-MD simulations of dendrite growth in the undercooled nickel melt of 1380 K. The corresponding movie (Movie S1) can be found in Supporting Information. In the early stages, around 2000 ps, the circular solid seed transforms into rhombic shapes, driven by preferential growth in <100> directions. This phenomenon aligns with observations from our previous MD and CG-MD simulations conducted with smaller systems.^106,163) Subsequently, at around 6000 ps, the solid-liquid interface in the <110> direction has negative curvature and the tip in the <100> direction begins to protrude. The protruding tip gradually enlarges and a typical 4-fold symmetrical dendrite structure appears after about 10000 ps.

Fig. 13. Snapshots of micrometer-scale coarse-grained (CG) molecular dynamics simulations of dendrite growth in the undercooled nickel melt of 1380 K. The calculation system consists of 16000553 CG beads, which is equivalent to 432014931 atoms. (Top) Green and white CG beads represent ones with fcc and liquid configurations identified by polyhedral template matching and (bottom) temperature distribution in the calculation system. (Online version in color.)

Focusing on the temperature distribution, the latent heat occurs at the solid-liquid interface and the temperature near the solid-liquid interface has recovered to the melting point. Therefore, a descending temperature gradient is generated from the solid phase towards the liquid phase (i.e., the negative temperature gradient) due to the undercooling on the liquid side. In the case of pure metals, a negative temperature gradient at the solid-liquid interface is required to enhance interfacial instability, leading to the formation of dendritic structures.^167,168) In order to investigate the effect of negative temperature gradient on the shape of the solidification interface, further calculations of solidification in a homogeneous temperature field was performed by adding the Langevin thermostat.¹⁶⁹⁾ Note that when using the Nosé-Hoover thermostat, which controls the average temperature of the entire system, local temperature gradients can occur in systems with heat generation or absorption. On the other hand, with the Langevin thermostat, where temperatures at specific locations are explicitly defined, it is possible to uniformly control the temperature of the entire system, even in cases involving heat generation or absorption. Figure 14 compares the solid structure obtained from the CG-MD simulations of solidification with and without the use of the Langevin thermostat. The simulation method is the same as previous one except for the addition of the Langevin thermostat. When using the Langevin thermostat, the temperature within the system is uniform, and thus, no negative temperature gradient occurs. As a result, the interface instability is not enhanced, and the dendritic structure do not emerge. Moreover, it is confirmed that when using the Langevin thermostat, the latent heat is forcibly extracted even from the solid nucleus, resulting in a faster crystal growth rate compared to the case without using the Langevin thermostat. Through comparative calculations using the Langevin thermostat, it is confirmed from atomistic viewpoint that a negative temperature gradient is necessary for dendrite growth in the crystalline growth of pure metal.

Fig. 14. Comparison of the solid structure obtained from the coarse-grained (CG) molecular dynamics simulations of solidification from the solid nuclei in the undercooled nickel melt of 1380 K (a) with the Langevin thermostat and (b) without the Langevin thermostat (same as in Fig. 13). (Top) Snapshots in which green and white CG beads represent ones with fcc and liquid configurations and (bottom) temperature distribution in the calculation system. (Online version in color.)

5.3. Application of Machine Learning Approaches to MD Simulation

Machine learning approaches have gained popularity and are now widely used across many fields. These approaches are also applied in various aspects of MD simulations,¹⁷⁰⁾ such as neural network potentials,^54,55,56) estimation of local atomic configuration,^{171,172,173,174)} and prediction of materials properties.^{175,176,177,178,179,180)} In particular, generative models are gaining attention in all fields as artificial intelligence tools for a new era. While there are not many examples of applying generative models to MD simulations yet, some pioneering studies have been reported.^149,150,151) For example, the variational autoencoder (VAE) encodes molecular structures into latent spaces and can subsequently decode them back.¹⁸¹⁾ Generative adversarial networks (GANs) has successfully predicted dynamic statistics, demonstrating high accuracy in tasks such as forecasting vibrational spectra for water molecules and simulating diffusion phenomena in polymeric polymers.^149,150) Moreover, the author’s group has proposed a novel approach using a combination of deep generative models and recurrent neural networks to predict multi-atomic cooperative phenomena at the atomic scale.¹⁵¹⁾ Here, a study of prediction of microstructure evolution of the polycrystalline iron at the atomic scale based on the deep generative model is presented in line with our previous report.¹⁵¹⁾ Due to space limitations, details of machine learning methodologies are not delved into in this paper. The methodology and implementation details of the machine learning procedures can be found in Ref. 151).

Figure 15 depicts the flowchart for predicting the microstructure evolution using atomic coordinates from the MD simulation as training data of the deep generative model. The calculation procedure is as follows. The preprocessed data of an MD simulation at equal time intervals are encoded into latent variables of 10 dimensions using the VAE model. Then, these latent variables are utilized in a long short-term memory (LSTM) model¹⁸²⁾ to predict future latent variables in the time series. Thereafter, the predicted latent variables are decoded using the VAE model to predict the microstructure at a future time not covered by the MD simulation. The unique aspect of this study lies in learning the temporal evolution trends of latent variables obtained by dimensionally compressing MD data using LSTM and in predicting the latent variables for future time.

Fig. 15. Flowchart for predicting the microstructure evolution using atomic coordinates from the MD simulation as training data of the deep generative model.¹⁵¹⁾ (Online version in color.)

This approach is applied to the data of MD simulation of grain growth in Fig. 5. The LSTM model was trained with 10-dimensional latent variables encoded from MD data at intervals of 100 ps from 1000 ps to 9000 ps. Subsequently, the LSTM model predicted the time evolution of latent variables up to 15000 ps based on the learned trends. These predicted latent variables in the future were then reconstructed using a decoder of the VAE model to predict microstructure evolution between 10000 ps and 15000 ps. It should be noted that there are no MD simulation results available beyond 10000 ps. Figure 16 displays snapshots of the microstructure obtained from the MD simulation up to 9000 ps and predicted microstructures up to 15000 ps. The predicted microstructures were reconstructed from predicted latent variables in the future time. The predicted microstructure remains relatively stable between 10000 and 15000 ps. Identifying the reason, whether it’s insufficient learning by the LSTM model or physical convergence due to the disappearance of small grains, is difficult at this stage.

Fig. 16. Snapshots of microstructure of polycrystalline bcc iron from the MD simulation up to 9000 ps (ground truth) and predicted microstructures up to 15000 ps. Predicted images are reconstructed from latent variables in the future time predicted by long short-term memory (LSTM) model.¹⁵¹⁾ (Online version in color.)

Next, an ongoing study of encoding and decoding atomic configuration from MD simulation of solidification is presented. Data from an MD simulation similar to the calculation in Fig. 13, but performed on a smaller system measuring 91.5 × 91.5 × 3.5 nm³,¹⁰⁶⁾ was used. Figure 17(a) shows snapshots of solidification from a seed crystal in the undercooled nickel melt at 1480 K, whose atomic configurations are used as training data. After converting the atomic configuration into array data with 100 × 100 elements as part of the data preprocessing, the array data is divided into 70% training data and 30% test data for training of the VAE model. Using the trained VAE model, all array data are compressed dimensionally into latent variables of 10 dimensions. Subsequently, t-distributed stochastic neighbor embedding (t-SNE)¹⁸³⁾ is applied to these 10-dimensional latent variables to visualize them in three dimensions, as shown in Fig. 17(b). Obtained latent variables are input to the decoder to restore the image of atomic configuration. Figure 17(c) shows images restored from 10-dimensional latent variables by the VAE decoder. From a macroscopic viewpoint, the restored images faithfully reproduce the original MD coordinates, suggesting that the encoding and decoding processes performed by the VAE are functioning as intended. However, it was difficult to identify the atomic positions since the probability of the existence of an atom is represented by the distribution of a real number between 0 and 1, as per the definition of VAE. In particular, the reconstruction of the liquid phase tends to be blurred at the atomic scale, and even with binarization, it is difficult to identify atomic position of the liquid phase. The extraction of quantitative atomic coordinates is still a challenge. Further enhancements in model refinement and parameter tuning are expected to enhance the accuracy of this approach.

Fig. 17. Encoding and decoding of atomic configuration from MD simulation of solidification. (a) Snapshots of solidification from a seed crystal in the undercooled nickel melt at 1480 K.¹⁰⁶⁾ (b) 10-dimensional latent variables encoded by the VAE model, which is visualized in the three-dimensional space by t-SNE. (c) Images restored from 10-dimensional latent variables by the VAE decoder. (Online version in color.)

5.4. Toward the Use of Experimental Data for Simulation

In recent years, data-driven approaches have gained momentum and are widely adopted in the field of materials science.¹⁸⁴⁾ In the computational materials science, numerous success stories have been demonstrated in the exploration of novel materials through high-throughput calculations.^{185,186,187,188,189,190)} Moreover, there have been attempts to incorporate experimental data into simulations. In recent years, data assimilation¹⁹¹⁾ has gained attention as one such method. Data assimilation is a mathematical approach that involves combining numerical simulations with observation data. Its purpose is to estimate the states and parameters of a system and to improve the accuracy of simulations. Data assimilation had achieved great success, primarily in geophysics,¹⁹²⁾ in its early stages. In recent years, data assimilation has been actively applied in the field of materials science. Examples include microstructure formation,^193,194) sintering,¹⁹⁵⁾ crystal growth of two-dimensional materials,¹⁹⁶⁾ thermal property predictions in solidification experiments,^197,198) evaluation of the finite-temperature magnetization,¹⁹⁹⁾ and parameter predictions of electrode reaction.²⁰⁰⁾ On the other hand, the spatiotemporal scale of MD simulations is significantly smaller than the typical experimental observation scale. As a result, there are not many research examples of data assimilation between MD simulations and experimental data. In this context, attempts such as data assimilation between X-ray diffraction data and MD simulations have been reported, particularly for the determination of accurate crystal structures.^201,202) The fusion of experimental data and simulations is undoubtedly becoming important in the future. Fusion methods utilizing data assimilation are currently under development and are expected to be widely expanded in the near future.

6. Conclusions

This review focused on the recent advances and contributions of MD simulations in the field of solidification. Main contributions of the MD simulation are derivation of solid-liquid interfacial properties and direct observation of nucleation and interfacial structures during solidification at atomistic scale. Regarding the derivation of the solid-liquid interfacial property, an overview of the capillary fluctuation method and the metadynamics-based method was given, and examples of studies on the derivation of solid-liquid interfacial energy for iron and iron-based alloy systems were presented. Subsequently, concerning MD simulations of nucleation and solidification, a comprehensive discussion was conducted, incorporating numerous research cases such as simultaneous nucleation and subsequent microstructure formation, heterogeneity in homogeneous nucleation, polymorphism, heterogeneous nucleation, and athermal nucleation. Not only were these success cases described, but also the issues related to the spatiotemporal scale of MD simulations were discussed. Attempts to overcome these challenges are presented in the chapter on Challenges and Perspectives. In the large-scale CG-MD simulation, dendrite growth was successfully achieved since there was a sufficient size for the emergence of a negative temperature gradient and a complex-shaped solid-liquid interface. To the best of the author’s knowledge, this study represents the first achievement of directly simulating dendrite growth solely through atomic-scale simulations without relying on any interfacial parameters. Moreover, a novel study of prediction of microstructure evolution at the atomic scale based on a combination of deep generative models and recurrent neural networks was presented. Importantly, decoding future-predicted latent variables by the LSTM model allows for obtaining microstructure predictions beyond the MD simulation timescale. This approach holds the potential to achieve a completely new method for accelerating MD simulations. The fusion of data-driven methods and MD simulations has just begun, and I anticipate that these new endeavors will contribute to the rapid development of the field of solidification in the future.

Supporting Information

Movie of micrometer-scale CG-MD simulation of dendritic growth. This material is available on the Journal website at https://doi.org/10.2355/isijinternational.ISIJINT-2024-010.

Acknowledgement

Part of this work was supported by Grant-in-Aid for Scientific Research (B) (No. 22H01754)) from Japan Society for the Promotion of Science, and the 30th ISIJ Research Promotion Grants from the Iron and Steel Institute of Japan. Part of the research was conducted using the Fujitsu PRIMERGY CX400M1/CX2550M5 (Oakbridge-CX) at the Information Technology Center, The University of Tokyo. Many of the studies presented here are the results of collaborative research with the following colleagues, and I appreciate their kind support and great contributions: Kensho Ueno, Satoru Fukuhara, Shin Okita, Shunsuke Orihara, Takuya Fujinaga, Loïc Chalamet, Kohei Sase, Munekazu Ohno, Tomohiro Takaki, Shinji Sakane, Eisuke Miyoshi, David Rodney and Tetsuo Mohri.

References

• 1)   W.  Kurz,  M.  Rappaz and  R.  Trivedi: Int. Mater. Rev., 66 (2021), 30. https://doi.org/10.1080/09506608.2020.1757894
• 2)   C. L. A.  Leung,  S.  Marussi,  R. C.  Atwood,  M.  Towrie,  P. J.  Withers and  P. D.  Lee: Nat. Commun., 9 (2018), 1355. https://doi.org/10.1038/s41467-018-03734-7
• 3)   H.  Yasuda,  K.  Morishita,  N.  Nakatsuka,  T.  Nishimura,  M.  Yoshiya,  A.  Sugiyama,  K.  Uesugi and  A.  Takeuchi: Nat. Commun., 10 (2019), 3183. https://doi.org/10.1038/s41467-019-11079-y
• 4)   Y.  Shibuta,  M.  Ohno and  T.  Takaki: JOM, 67 (2015), 1793. https://doi.org/10.1007/s11837-015-1452-2
• 5)   Y.  Shibuta,  M.  Ohno and  T.  Takaki: Adv. Theory Simul., 1 (2018), 1800065. https://doi.org/10.1002/adts.201800065
• 6)   T.  Takaki: ISIJ Int., 54 (2014), 437. https://doi.org/10.2355/isijinternational.54.437
• 7)   T.  Takaki,  S.  Sakane,  M.  Ohno,  Y.  Shibuta,  T.  Shimokawabe and  T.  Aoki: ISIJ Int., 56 (2016), 1427. https://doi.org/10.2355/isijinternational.ISIJINT-2016-156
• 8)   S.  Sakane,  T.  Takaki,  M.  Ohno,  Y.  Shibuta,  T.  Shimokawabe and  T.  Aoki: J. Cryst. Growth, 483 (2018), 147. https://doi.org/10.1016/j.jcrysgro.2017.11.029
• 9)   M.  Ohno: ISIJ Int., 60 (2020), 2745. https://doi.org/10.2355/isijinternational.ISIJINT-2020-174
• 10)   D.  Tourret,  H.  Liu and  J. LLorca: Prog. Mater. Sci., 123 (2022), 100819. https://doi.org/10.1016/j.pmatsci.2021.100810
• 11)   M.  Yamazaki,  Y.  Natsume,  H.  Harada and  K.  Ohsasa: ISIJ Int., 46 (2006), 903. https://doi.org/10.2355/isijinternational.46.903
• 12)   D.  Tourret,  A. J.  Clarke,  S. D.  Imhoff,  P. J.  Gibbs,  J. W.  Gibbs and  A.  Karma: JOM, 67 (2015), 1776. https://doi.org/10.1007/s11837-015-1444-2
• 13)   Y.  Shibuta: Mater. Trans., 60 (2019), 180. https://doi.org/10.2320/matertrans.ME201712
• 14)   Y.  Shibuta,  S.  Sakane,  E.  Miyoshi,  S.  Okita,  T.  Takaki and  M.  Ohno: Nat. Commun., 8 (2017), 10. https://doi.org/10.1038/s41467-017-00017-5
• 15)   Y.  Shibuta,  S.  Sakane,  E.  Miyoshi,  T.  Takaki and  M.  Ohno: Model. Simul. Mater. Sci. Eng., 27 (2019), 054002. https://doi.org/10.1088/1361-651X/ab1d28
• 16)   J. J.  Hoyt,  M.  Asta and  A.  Karma: Phys. Rev. Lett., 86 (2001), 5530. https://doi.org/10.1103/PhysRevLett.86.5530
• 17)   D. Y.  Sun,  M.  Asta and  J. J.  Hoyt: Phys. Rev. B, 69 (2004), 174103. https://doi.org/10.1103/PhysRevB.69.174103
• 18)   X.-M.  Bai and  M.  Li: J. Chem. Phys., 124 (2006), 124707. https://doi.org/10.1063/1.2184315
• 19)   Y.  Watanabe,  Y.  Shibuta and  T.  Suzuki: ISIJ Int., 50 (2010), 1158. https://doi.org/10.2355/isijinternational.50.1158
• 20)   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
• 21)   K.  Ueno and  Y.  Shibuta: Comput. Mater. Sci., 167 (2019), 1. https://doi.org/10.1016/j.commatsci.2019.05.023
• 22)   L.  Wang and  J. J.  Hoyt: Acta Mater., 219 (2021), 117259. https://doi.org/10.1016/j.actamat.2021.117259
• 23)   K.  Ueno and  Y.  Shibuta: Materialia, 4 (2018), 553. https://doi.org/10.1016/j.mtla.2018.11.011
• 24)   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
• 25)   S.  Kavousi,  B. R.  Novak,  J. J.  Hoyt and  D.  Moldovan: Comput. Mater. Sci., 184 (2020), 109854. https://doi.org/10.1016/j.commatsci.2020.109854
• 26)   J. J.  Hoyt,  M.  Asta and  A.  Karma: Mater. Sci. Eng. R, 41 (2003), 121. https://doi.org/10.1016/S0927-796X(03)00036-6
• 27)   Y.  Shibuta,  K.  Oguchi and  T.  Suzuki: ISIJ Int., 52 (2012), 2205. https://doi.org/10.2355/isijinternational.52.2205
• 28)   Y.  Shibuta,  K.  Oguchi and  M.  Ohno: Scr. Mater., 86 (2014), 20. https://doi.org/10.1016/j.scriptamat.2014.04.021
• 29)   T.  Takaki,  S.  Sakane,  M.  Ohno,  Y.  Shibuta and  T.  Aoki: Comp. Mater. Sci., 171 (2020), 109209. https://doi.org/10.1016/j.commatsci.2019.109209
• 30)   S.  Sakane,  T.  Takaki,  M.  Ohno,  Y.  Shibuta and  T.  Aoki: Comp. Mater. Sci., 178 (2020), 109639. https://doi.org/10.1016/j.commatsci.2020.109639
• 31)   T.  Fujinaga and  Y.  Shibuta: Comp. Mater. Sci., 164 (2019), 74. https://doi.org/10.1016/j.commatsci.2019.03.061
• 32)   T.  Fujinaga,  Y.  Watanabe and  Y.  Shibuta: Comp. Mater. Sci., 182 (2020), 109763. https://doi.org/10.1016/j.commatsci.2020.109763
• 33)   M. P.  Allen and  D. J.  Tildesley: Computer Simulation of Liquids, 2nd ed., Oxford University Press, Oxford, (2017), 71. https://doi.org/10.1093/oso/9780198803195.001.0001
• 34)   W.  Andreoni and  S.  Yip: Handbook of Materials Modeling (2nd ed.), Springer, Cham, (2020), 493. https://doi.org/10.1007/978-3-319-44677-6
• 35)   J.-L.  Fattebert and  F.  Gygi: Comput. Phys. Commun., 162 (2004), 24. https://doi.org/10.1016/j.cpc.2004.06.068
• 36)   F.  Shimojo,  A.  Nakano,  R. K.  Kalia and  P.  Vashishta: Phys. Rev. E, 77 (2008), 066103. https://doi.org/10.1103/PhysRevE.77.066103
• 37)   Y.  Shibuta,  R.  Arifin,  K.  Shimamura,  T.  Oguri,  F.  Shimojo and  S.  Yamaguchi: Chem. Phys. Lett., 565 (2013), 92. https://doi.org/10.1016/j.cplett.2013.02.038
• 38)   T.  Oguri,  K.  Shimamura,  Y.  Shibuta,  F.  Shimojo and  S.  Yamaguchi: J. Phys. Chem. C, 117 (2013), 9983. https://doi.org/10.1021/jp403006m
• 39)   R.  Sato,  S.  Ohkuma,  Y.  Shibuta,  F.  Shimojo and  S.  Yamaguchi: J. Phys. Chem. C, 119 (2015), 28925. https://doi.org/10.1021/acs.jpcc.5b09026
• 40)   S.  Sugiura,  Y.  Shibuta,  K.  Shimamura,  M.  Misawa,  F.  Shimojo and  S.  Yamaguchi: Solid State Ion., 285 (2016), 209. https://doi.org/10.1016/j.ssi.2015.06.008
• 41)   S.  Tateyama,  Y.  Shibuta and  T.  Suzuki: ISIJ Int., 50 (2010), 1211. https://doi.org/10.2355/isijinternational.50.1211
• 42)   S.  Tateyama,  Y.  Shibuta,  T.  Kumagai and  T.  Suzuki: ISIJ Int., 51 (2011), 1710. https://doi.org/10.2355/isijinternational.51.1710
• 43)   S.  Okita and  Y.  Shibuta: ISIJ Int., 56 (2016), 2199. https://doi.org/10.2355/isijinternational.ISIJINT-2016-408
• 44)   S.  Okita,  E.  Miyoshi,  S.  Sakane,  T.  Takaki,  M.  Ohno and  Y.  Shibuta: Acta Mater., 153 (2018), 108. https://doi.org/10.1016/j.actamat.2018.04.060
• 45)   J. E.  Lennard-Jones: Proc. Phys. Soc., 43 (1931), 461. https://doi.org/10.1088/0959-5309/43/5/301
• 46)   P. M.  Morse: Phys. Rev., 34 (1929), 57. https://doi.org/10.1103/PhysRev.34.57
• 47)   J.  Tersoff: Phys. Rev. B, 37 (1988), 6991. https://doi.org/10.1103/PhysRevB.37.6991
• 48)   A. C. T.  van Duin,  S.  Dasgupta,  F.  Lorant and  W. A.  Goddard., III: J. Phys. Chem. A, 105 (2001), 9396. https://doi.org/10.1021/jp004368u
• 49)   S.  Daw and  M. I.  Baskes: Phys. Rev. B, 29 (1984), 6443. https://doi.org/10.1103/PhysRevB.29.6443
• 50)   M. I.  Baskes: Phys. Rev. B, 46 (1992), 2727. https://doi.org/10.1103/PhysRevB.46.2727
• 51)   B.-J.  Lee,  J.-H.  Shim and  M. I.  Baskes: Phys. Rev. B, 68 (2003), 144112. https://doi.org/10.1103/PhysRevB.68.144112
• 52)   C. A.  Becker,  F.  Tavazza,  Z. T.  Trautt and  R. A.  Buarque de Macedoc: Curr. Opin. Solid. State. Mater. Sci., 17 (2013), 277. https://doi.org/10.1016/j.cossms.2013.10.001
• 53)   E. B.  Tadmor,  R. S.  Elliott,  J. P.  Sethna,  R. E.  Miller and  C. A.  Becker: JOM, 63 (2011), 17. https://doi.org/10.1007/s11837-011-0102-6
• 54)   J.  Behler and  M.  Parrinello: Phys. Rev. Lett., 98 (2007), 146401. https://doi.org/10.1103/PhysRevLett.98.146401
• 55)   E.  Kocer,  T. W.  Ko and  J.  Behler: Annu. Rev. Phys. Chem., 73 (2022), 163. https://doi.org/10.1146/annurev-physchem-082720-034254
• 56)   S.  Takamoto,  C.  Shinagawa,  D.  Motoki,  K.  Nakago,  W.  Li,  I.  Kurata,  T.  Watanabe,  Y.  Yayama,  H.  Iriguchi,  Y.  Asano,  T.  Onodera,  T.  Ishii,  T.  Kudo,  H.  Ono,  R.  Sawada,  R.  Ishitani,  M.  Ong,  T.  Yamaguchi,  T.  Kataoka,  A.  Hayashi,  N.  Charoenphakdee and  T.  Ibuka: Nat. Commum., 13 (2022), 2991. https://doi.org/10.1038/s41467-022-30687-9
• 57)   S.  Plimpton: J. Comput. Phys., 117 (1995), 1. https://doi.org/10.1006/jcph.1995.1039
• 58)   A.  Stukowski: Model. Simul. Mater. Sci. Eng., 18 (2010), 015012. https://doi.org/10.1088/0965-0393/18/1/015012
• 59)   W.  Kurz and  D. J.  Fisher: Fundamentals of Solidification 4th revised ed., Trans Tech Publication, Aedermannsdorf, (1998), 21. ISBN978-0878498048
• 60)   J. A.  Dantzig and  M.  Rappaz: Solidification, EPFL Press, Lausanne, (2009), 287. ISBN978-0849382383
• 61)   D.  Turnbull and  R. E.  Cech: J. Appl. Phys., 21 (1950), 804. https://doi.org/10.1063/1.1699763
• 62)   D.  Turnbull: J. Appl. Phys., 21 (1950), 1022. https://doi.org/10.1063/1.1699435
• 63)   M. E.  Glicksman and  C. L.  Vold: Acta Metall., 17 (1969), 1. https://doi.org/10.1016/0001-6160(69)90157-6
• 64)   L.  Wang,  X.  Lin,  M.  Wang and  W.  Huang: Scr. Mater., 69 (2013), 1. https://doi.org/10.1016/j.scriptamat.2013.03.028
• 65)   L.  Wang,  X.  Lin,  M.  Wang and  W.  Huang: J. Cryst. Growth, 406 (2014), 85. https://doi.org/10.1016/j.jcrysgro.2014.08.017
• 66)   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
• 67)   A. S.  Skapski: Acta Metall., 4 (1956), 576. https://doi.org/10.1016/0001-6160(56)90159-6
• 68)   Y.  Waseda and  W. A.  Miller: Trans. Jpn. Inst. Met., 19 (1978), 546. https://doi.org/10.2320/matertrans1960.19.546
• 69)   R. H.  Ewing: J. Cryst. Growth, 11 (1971), 221. https://doi.org/10.1016/0022-0248(71)90088-1
• 70)   F.  Spaepen: Acta Metall., 23 (1975), 729. https://doi.org/10.1016/0001-6160(75)90056-5
• 71)   L.  Gránásy,  M.  Tegze and  A.  Ludwig: Mater. Sci. Eng. A, 133 (1991), 577. https://doi.org/10.1016/0921-5093(91)90138-D
• 72)   J. Q.  Broughton and  G. H.  Gilmer: J. Chem. Phys., 84 (1986), 5759. https://doi.org/10.1063/1.449884
• 73)   R. L.  Davidchack and  B. B.  Laird: J. Chem. Phys., 118 (2003), 7651. https://doi.org/10.1063/1.1563248
• 74)   N. T.  Brown,  E.  Martinez and  J.  Qu: Acta Mater., 129 (2017), 83. https://doi.org/10.1016/j.actamat.2017.02.033
• 75)   K.  Ueno and  Y.  Shibuta: ISIJ Int., 60 (2020), 2301. https://doi.org/10.2355/isijinternational.ISIJINT-2019-769
• 76)   A.  Laio and  M.  Parrinello: Proc. Natl. Acad. Sci. USA, 99 (2002), 12562. https://doi.org/10.1073/pnas.202427399
• 77)   S.  Angioletti-Uberti,  M.  Ceriotti,  P. D.  Lee and  M. W.  Finnis: Phys. Rev. B, 81 (2010), 125416. https://doi.org/10.1103/PhysRevB.81.125416
• 78)   S.  Angioletti-Uberti: J. Phys. Condens. Matter, 23 (2011), 435008. https://doi.org/10.1088/0953-8984/23/43/435008
• 79)   B.  Cheng,  G. A.  Tribello and  M.  Ceriotti: Phys. Rev. B, 92 (2015), 180102. https://doi.org/10.1103/PhysRevB.92.180102
• 80)   K.  Ueno,  S.  Fukuhara and  Y.  Shibuta: Mater. Trans., 63 (2022), 209. https://doi.org/10.2320/matertrans.MT-M2021156
• 81)   R. L.  Davidchack,  J. R.  Morris and  B. B.  Laird: J. Chem. Phys., 125 (2006), 094710. https://doi.org/10.1063/1.2338303
• 82)   G. W.  Kim,  T.  Takaki,  Y.  Shibuta,  S.  Sakane,  K.  Matsuura and  M.  Ohno: Comp. Mater. Sci., 162 (2019), 76. https://doi.org/10.1016/j.commatsci.2019.02.027
• 83)   W. R.  Fehlner and  S. H.  Vosko: Can. J. Phys., 54 (1976), 2159. https://doi.org/10.1139/p76-256
• 84)   G.  Bonny,  R. C.  Pasianot,  D.  Terentyev and  L.  Malerba: Philos. Mag., 91 (2011), 1724. https://doi.org/10.1080/14786435.2010.545780
• 85)   G.  Bonny,  R. C.  Pasianot and  L.  Malerba: Model. Simul. Mater. Sci. Eng., 17 (2009), 025010. https://doi.org/10.1088/0965-0393/17/2/025010
• 86)   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
• 87)   J. R.  Morris: Phys. Rev. B, 66 (2002), 144104. https://doi.org/10.1103/PhysRevB.66.144104
• 88)   R.  Freitas,  T.  Frolov and  M.  Asta: Phys. Rev. E, 96 (2017), 043308. https://doi.org/10.1103/PhysRevE.96.043308
• 89)   K.  Ueno and  Y.  Shibuta: IOP Conf. Ser.: Mater. Sci. Eng., 861 (2020), 012064. https://doi.org/10.1088/1757-899X/861/1/012064
• 90)   R. L.  Davidchack and  B. B.  Laird: Phys. Rev. Lett., 94 (2005), 086102. https://doi.org/10.1103/PhysRevLett.94.086102
• 91)   Y.  Mu,  A.  Houk and  X.  Song: J. Phys. Chem. B, 109 (2005), 6500. https://doi.org/10.1021/jp046289e
• 92)   K. M.  Bal,  S.  Fukuhara,  Y.  Shibuta and  E. C.  Neyts: J. Chem. Phys., 153 (2020), 114118. https://doi.org/10.1063/5.0020240
• 93)   S.  Fukuhara and  Y.  Shibuta: Nanoscale Adv., 3 (2021), 6191. https://doi.org/10.1039/D1NA00377A
• 94)   K. M.  Bal and  E. C.  Neyts: J. Chem. Theory Comput., 11 (2015), 4545. https://doi.org/10.1021/acs.jctc.5b00597
• 95)   S.  Fukuhara,  K. M.  Bal,  E. C.  Neyts and  Y.  Shibuta: Comput. Mater. Sci., 177 (2020), 109581. https://doi.org/10.1016/j.commatsci.2020.109581
• 96)   H.  Ebina,  S.  Fukuhara and  Y.  Shibuta: Comput. Mater. Sci., 196 (2021), 110577. https://doi.org/10.1016/j.commatsci.2021.110577
• 97)   A. F.  Voter: Phys. Rev. Lett., 78 (1997), 3908. https://doi.org/10.1103/PhysRevLett.78.3908
• 98)   G. A.  Tribello,  M.  Bonomi,  D.  Branduardi,  C.  Camilloni and  G.  Bussi: Comput. Phys. Commun., 185 (2014), 604. https://doi.org/10.1016/j.cpc.2013.09.018
• 99)   Y.  Waseda and  Y.  Shiraishi: Bull. Res. Inst. Miner. Dressing Metall., Tohoku Univ., 33 (1977), 109 (in Japanese).
• 100)   K. F.  Kelton: Solid State Phys., 45 (1991), 75. https://doi.org/10.1016/S0081-1947(08)60144-7
• 101)   B.  Vinet,  L.  Magnusson,  H.  Fredriksson and  P. J.  Desre: J. Colloid Interface Sci., 255 (2002), 363. https://doi.org/10.1006/jcis.2002.8627
• 102)   Z.  Jian,  N.  Li,  M.  Zhu,  J.  Chen,  F.  Chang and  W.  Jie: Acta Mater., 60 (2012), 3590. https://doi.org/10.1016/j.actamat.2012.02.038
• 103)   Y.  Sun,  F.  Zhang,  H.  Song,  M. I.  Mendelev,  C. Z.  Wang and  K. M.  Ho: J. Chem. Phys., 149 (2018), 174501. https://doi.org/10.1063/1.5048781
• 104)   K.  Mondal,  A.  Kumar,  G.  Gupta and  B. S.  Murty: Acta Mater., 57 (2009), 3422. https://doi.org/10.1016/j.actamat.2009.03.051
• 105)   M.  Ohno,  Y.  Oka,  S.  Sakane,  Y.  Shibuta and  T.  Takaki: Phys. Rev. E, 101 (2020), 052121. https://doi.org/10.1103/PhysRevE.101.052121
• 106)   Y.  Nagatsuma,  M.  Ohno,  T.  Takaki and  Y.  Shibuta: Nanomaterials, 11 (2021), 2308. https://doi.org/10.3390/nano11092308
• 107)   G. C.  Sosso,  J.  Chen,  S. J.  Cox,  M.  Fitzner,  P.  Pedevilla,  A.  Zen and  A.  Michaelides: Chem. Rev., 116 (2016), No. 12, 7078. https://doi.org/10.1021/acs.chemrev.5b00744
• 108)   S.  Auer and Frenkel: Nature, 409 (2001), 1020. https://doi.org/10.1038/35059035
• 109)   R. S.  Aga,  J. R.  Morris,  J. J.  Hoyt and  M.  Mendelev: Phys. Rev. Lett., 96 (2006), 245701. https://doi.org/10.1103/PhysRevLett.96.245701
• 110)   J.  Bokeloh,  R. E.  Rozas,  J.  Horbach and  G.  Wilde: Phys. Rev. Lett., 107 (2011), 145701. https://doi.org/10.1103/PhysRevLett.107.145701
• 111)   F. J.  Cherne,  M. I.  Baskes,  R. B.  Schwarz,  S. G.  Srinivasan and  W.  Klein: Model. Simul. Mater. Sci. Eng., 12 (2004), 1063. https://doi.org/10.1088/0965-0393/12/6/001
• 112)   J.  Russo and  H.  Tanaka: MRS Bull., 41 (2016), 369. https://doi.org/10.1557/mrs.2016.84
• 113)   B.  Doliwa and  A.  Heuer: Phys. Rev. Lett., 80 (1998), 4915. https://doi.org/10.1103/PhysRevLett.80.4915
• 114)   W.  Kob: J. Phys.: Condens. Matter, 11 (1999), R85. https://doi.org/10.1088/0953-8984/11/10/003
• 115)   Z. Y.  Hou,  L. X.  Liu,  R. S.  Liu,  Z. A.  Tian and  J. G.  Wang: Chem. Phys. Lett., 491 (2010), 172. https://doi.org/10.1016/j.cplett.2010.04.003
• 116)   P. R.  ten Wolde,  M. J.  Ruiz-Montero and  D.  Frenkel: Phys. Rev. Lett., 75 (1995), 2714. https://doi.org/10.1103/PhysRevLett.75.2714
• 117)   C.  Desgranges and  J.  Delhommelle: J. Am. Chem. Soc., 128 (2006), 10368. https://doi.org/10.1021/ja063218f
• 118)   S.  Orihara,  Y.  Shibuta and  T.  Mohri: Mater. Trans., 61 (2020), 750. https://doi.org/10.2320/matertrans.MT-M2019353
• 119)   Y.  Shibuta,  K.  Oguchi,  T.  Takaki and  M.  Ohno: Sci. Rep., 5 (2015), 13534. https://doi.org/10.1038/srep13534
• 120)   Y.  Shibuta,  S.  Sakane,  T.  Takaki and  M.  Ohno: Acta Mater., 105 (2016), 328. https://doi.org/10.1016/j.actamat.2015.12.033
• 121)   X.  Sui,  Y.  Cheng,  N.  Zhou,  B.  Tang and  L.  Zhou: CrystEngComm, 20 (2018), 3569. https://doi.org/10.1039/C8CE00767E
• 122)   M.  Papanikolaou,  K.  Salonitis,  M.  Jolly and  M.  Frank: Metals, 9 (2019), 1217. https://doi.org/10.3390/met9111217
• 123)   L.  Separdar,  J. P.  Rino and  E. D.  Zanotto: Comput. Mater. Sci., 187 (2021), 110124. https://doi.org/10.1016/j.commatsci.2020.110124
• 124)   Y.  Gao,  B.  Wang,  J.  Huang,  T.  Gao,  W.  Yang,  Q.  Xie and  Q.  Chen: Cryst. Growth Des., 22 (2022), 2417. https://doi.org/10.1021/acs.cgd.1c01499
• 125)   M. W.  Finnis and  J. E.  Sinclair: Philos. Mag. A, 50 (1984), 45. https://doi.org/10.1080/01418618408244210
• 126)   K.  Oguchi,  Y.  Shibuta and  T.  Suzuki: J. Jpn. Inst. Met., 76 (2012), 462 (in Japanese). https://doi.org/10.2320/jinstmet.76.462
• 127)   J. K.  Mackenzie: Biometrika, 45 (1958), 229. https://doi.org/10.1093/biomet/45.1-2.229
• 128)   D. L.  Olmsted,  S. M.  Foiles and  E. A.  Holm: Acta. Mater., 57 (2009), 3694. https://doi.org/10.1016/j.actamat.2009.04.007
• 129)   Y.  Shibuta,  S.  Takamoto and  T.  Suzuki: ISIJ Int., 48 (2008), 1582. https://doi.org/10.2355/isijinternational.48.1582
• 130)   W.  Ostwald: Z. Phys. Chem., 22 (1897), 289. https://doi.org/10.1515/zpch-1897-2233
• 131)   G.  Tammann: Aggregatzustände, Verlag Von Leopold Voss, Leipzig, (1922), 234. ISBN978-1120566508
• 132)   K. N.  Ishihara,  M.  Maeda and  P. H.  Shingu: Acta Metall., 33 (1985), 2113. https://doi.org/10.1016/0001-6160(85)90172-5
• 133)   H.  Shingu and  K.  Ishihara: Bull. JIM, 25 (1986), 16 (in Japanese). https://doi.org/10.2320/materia1962.25.16
• 134)   H.  Wang,  H.  Gould and  W.  Klein: Phys. Rev. E, 76 (2007), 031604. https://doi.org/10.1103/PhysRevE.76.031604
• 135)   G. P.  Purja Pun and  Y.  Mishin: Philos. Mag., 89 (2009), 3245. https://doi.org/10.1080/14786430903258184
• 136)   J. J.  Hoyt,  M.  Asta,  T.  Haxhimali,  A.  Karma,  R. E.  Napolitano,  R.  Trivedi,  B. B.  Laird and  J. R.  Morris: MRS Bull., 29 (2004), 935. https://doi.org/10.1557/mrs2004.263
• 137)   H.  Men and  Z.  Fan: Comput. Mater. Sci., 85 (2014), 1. https://doi.org/10.1016/j.commatsci.2013.12.042
• 138)   Z.  Fan,  H.  Men,  Y.  Wang and  Z.  Que: Metals, 11 (2021), 478. https://doi.org/10.3390/met11030478
• 139)   Z.  Fan and  H.  Men: Metals, 12 (2022), 1547. https://doi.org/10.3390/met12091547
• 140)   M. A.  Easton,  M.  Qian,  A.  Prasad and  D. H.  StJohn: Curr. Opin. Solid State Mater. Sci., 20 (2016), 13. https://doi.org/10.1016/j.cossms.2015.10.001
• 141)   A. L.  Greer: J. Chem. Phys., 145 (2016), 211704. https://doi.org/10.1063/1.4968846
• 142)   T. E.  Quested and  A. L.  Greer: Acta Mater., 53 (2005), 2683. https://doi.org/10.1016/j.actamat.2005.02.028
• 143)   G. M.  Torrie and  J. P.  Valleau: J. Comput. Phys., 23 (1977), 187. https://doi.org/10.1016/0021-9991(77)90121-8
• 144)   B. A.  Berg and  T.  Neuhaus: Phys. Rev. Lett., 68 (1992), 9. https://doi.org/10.1103/PhysRevLett.68.9
• 145)   M. S.  Sørensen and  A. F.  Voter: J. Chem. Phys., 112 (2000), 9599. https://doi.org/10.1063/1.481576
• 146)   A. F.  Voter: Phys. Rev. B, 57 (1998), R13985. https://doi.org/10.1103/PhysRevB.57.R13985
• 147)   J.  Li,  S.  Sarkar,  W. T.  Cox,  T. J.  Lenosky,  E.  Bitzek and  Y.  Wang: Phys. Rev. B, 84 (2011), 054103. https://doi.org/10.1103/PhysRevB.84.054103
• 148)   A.  Ishii,  S.  Ogata,  H.  Kimizuka and  J.  Li: Phys. Rev B, 85 (2012), 064303. https://doi.org/10.1103/PhysRevB.85.064303
• 149)   R.  Kawada,  K.  Endo,  D.  Yuhara and  K.  Yasuoka: Soft Matter, 18 (2022), 8446. https://doi.org/10.1039/D2SM00852A
• 150)   R.  Kawada,  K.  Endo,  K.  Yasuoka,  H.  Kojima and  N.  Matubayasi: J. Chem. Inf. Model., 63 (2023), 76. https://doi.org/10.1021/acs.jcim.2c01316
• 151)   K.  Sase and  Y.  Shibuta: Acta Mater., 259 (2023), 119295. https://doi.org/10.1016/j.actamat.2023.119295
• 152)   E.  Miyoshi,  T.  Takaki,  Y.  Shibuta and  M.  Ohno: Comput. Mater. Sci., 152 (2018), 118. https://doi.org/10.1016/j.commatsci.2018.05.046
• 153)   S.  Izvekov and  G. A.  Voth: J. Phys. Chem. B, 109 (2005), 2469. https://doi.org/10.1021/jp044629q
• 154)   W. G.  Noid: J. Chem. Phys., 139 (2013), 090901. https://doi.org/10.1063/1.4818908
• 155)   M.  Pasi,  R.  Lavery and  N.  Ceres: J. Chem. Theory Comput., 9 (2013), 785. https://doi.org/10.1021/ct3007925
• 156)   G.  Rossi,  L.  Monticelli,  S. R.  Puisto,  I.  Vattulainen and  T.  Ala-Nissila: Soft Matter, 7 (2011), 698. https://doi.org/10.1039/C0SM00481B
• 157)   K.-H.  Shen,  M.  Fan and  L. M.  Hall: Macromolecules, 54 (2021), 2031. https://doi.org/10.1021/acs.macromol.0c02557
• 158)   N.  Korolev,  D.  Luo,  A. P.  Lyubartsev and  L.  Nordenskiöld: Polymers, 6 (2014), 1655. https://doi.org/10.3390/polym6061655
• 159)   J. J.  Uusitalo,  H. I.  Ingólfsson,  P.  Akhshi,  P. D.  Tieleman and  S. J.  Marrink: J. Chem. Theory Comput., 11 (2015), 3932. https://doi.org/10.1021/acs.jctc.5b00286
• 160)   R. D.  Hills, Jr.,  L.  Lu and  G. A.  Voth: Plos Comput. Biology, 6 (2010), e1000827. https://doi.org/10.1371/journal.pcbi.1000827
• 161)   A. M.  Dongare: Philo. Mag., 94 (2014), 3877. https://doi.org/10.1080/14786435.2014.961992
• 162)   G.  Agarwal,  R. R.  Valisetty,  R. R.  Namburu,  A. M.  Rajendran and  A. M.  Dongare: Sci. Rep., 7 (2017), 12376. https://doi.org/10.1038/s41598-017-12340-4
• 163)   L.  Chalamet,  D.  Rodney and  Y.  Shibuta: Comput. Mater. Sci., 228 (2023), 112306. https://doi.org/10.1016/j.commatsci.2023.112306
• 164)   P.  Delafrouz and  H. N.  Pishkenari: J. Theoretical Appl. Mech., 56 (2018), 601. https://doi.org/10.15632/jtam-pl.56.3.601
• 165)   S.  Nosé: J. Chem. Phys., 81 (1984), 511. https://doi.org/10.1063/1.447334
• 166)   W. G.  Hoover: Phys. Rev. A, 31 (1985), 1695. https://doi.org/10.1103/PhysRevA.31.1695
• 167)   W. W.  Mullins and  R. F.  Sekerka: J. Appl. Phys., 35 (1964), 444. https://doi.org/10.1063/1.1713333
• 168)   Y.  Saito: Materia Japan, 49 (2010), 431 (in Japanese). https://doi.org/10.2320/materia.49.431
• 169)   T.  Schneider and  E.  Stoll: Phys. Rev. B, 17 (1978), 1302. https://doi.org/10.1103/PhysRevB.17.1302
• 170)   F.  Musil,  A.  Grisafi,  A. P.  Bartók,  C.  Ortner,  G.  Csányi and  M.  Ceriotti: Chem. Rev., 121 (2021), 9759. https://doi.org/10.1021/acs.chemrev.1c00021
• 171)   T.  Fukuya and  Y.  Shibuta: Comput. Mater. Sci., 184 (2020), 109880. https://doi.org/10.1016/j.commatsci.2020.109880
• 172)   H. W.  Chung,  R.  Freitas,  G.  Cheon and  E. J.  Reed: Phys. Rev. Mater., 6 (2022), 043801. https://doi.org/10.1103/PhysRevMaterials.6.043801
• 173)   Y.  Deng,  Y.  Wang,  K.  Xu and  Y.  Wang: J. Chem. Theory Comput., 19 (2023), 8332. https://doi.org/10.1021/acs.jctc.3c00838
• 174)   A.  Allera,  A. M.  Goryaeva,  P.  Lafourcade,  J.-B.  Maillet and  M.-C.  Marinica: Comput. Mater. Sci., 231 (2024), 112535. https://doi.org/10.1016/j.commatsci.2023.112535
• 175)   M.  Tsunawaki,  S.  Fukuhara and  Y.  Shibuta: Mater. Trans., 62 (2021), 829. https://doi.org/10.2320/matertrans.MT-M2021032
• 176)   L.  Zhang,  K.  Qian,  B. W.  Schuller and  Y.  Shibuta: Metals, 11 (2021), 922. https://doi.org/10.3390/met11060922
• 177)   L.  Zhang,  K.  Qian,  J.  Huang,  M.  Liu and  Y.  Shibuta: J. Mater. Res. Tech., 13 (2021), 2043. https://doi.org/10.1016/j.jmrt.2021.06.021
• 178)   Q.  Kong and  Y.  Shibuta: Mater. Trans., 64 (2023), 1241. https://doi.org/10.2320/matertrans.MT-M2022204
• 179)   K.  Noda and  Y.  Shibuta: Comput. Mater. Sci., 229 (2023), 112448. https://doi.org/10.1016/j.commatsci.2023.112448
• 180)   J.  Yu,  F.  Yu,  Q.  Fu,  G.  Zhao,  C.  Gong,  M.  Wang and  Q.  Zhang: Nanomaterials, 13 (2023), 968. https://doi.org/10.3390/nano13060968
• 181)   A.  Karthikeyan and  U. D.  Priyakumar: J. Chem. Sci., 134 (2022), 2. https://doi.org/10.1007/s12039-021-01995-2
• 182)   S.  Hochreiter and  J.  Schmidhuber: Neural Comput., 9 (1997), 1735. https://doi.org/10.1162/neco.1997.9.8.1735
• 183)   L.  van der Maaten and  G.  Hinton: J. Mach. Learn. Res., 9 (2008), 2579.
• 184)   C.  Draxl and  M.  Scheffler: MRS Bull., 43 (2018), 676. https://doi.org/10.1557/mrs.2018.208
• 185)   S.  Curtarolo,  G. L. W.  Hart,  M. B.  Nardelli,  N.  Mingo,  S.  Sanvito and  O.  Levy: Nat. Mater., 12 (2013), 191. https://doi.org/10.1038/nmat3568
• 186)   A.  Jain,  S. P.  Ong,  G.  Hautier,  W.  Chen,  W. D.  Richards,  S.  Dacek,  S.  Cholia,  D.  Gunter,  D.  Skinner,  G.  Ceder and  K. A.  Persson: APL Mater., 1 (2013), 011002. https://doi.org/10.1063/1.4812323
• 187)   J. J.  de Pablo,  N. E.  Jackson,  M. A.  Webb,  L.-Q.  Chen,  J. E.  Moore,  D.  Morgan,  R.  Jacobs,  T.  Pollock,  D. G.  Schlom,  E. S.  Toberer,  J.  Analytis,  I.  Dabo,  D. M.  DeLongchamp,  G. A.  Fiete,  G. M.  Grason,  G.  Hautier,  Y.  Mo,  K.  Rajan,  E. J.  Reed,  E.  Rodriguez,  V.  Stevanovic,  J.  Suntivich,  K.  Thornton and  J.-C.  Zhao: npj Comput. Mater., 5 (2019), 41.
• 188)   G.  Brunin,  F.  Ricci,  V.-A.  Ha,  G.-M.  Rignanese and  G.  Hautier: npj Comput. Mater., 5 (2019), 63. https://doi.org/10.1038/s41524-019-0200-5
• 189)   A.  Takahashi,  Y.  Kumagai,  J.  Miyamoto,  Y.  Mochizuki and  F.  Oba: Phys. Rev. Mater., 4 (2020), 103801. https://doi.org/10.1103/PhysRevMaterials.4.103801
• 190)   S.-H.  Jang,  R.  Jalem and  Y.  Tateyama: J. Phys. Chem. A, 127 (2023), 5734. https://doi.org/10.1021/acs.jpca.3c00076
• 191)   W.  Lahoz,  B.  Khattatov and  R.  Menard: Data Assmilation - Making Sense of Observations, Springer, Berlin, (2010), 13. ISBN 978-3540747031
• 192)   M.  Ghil and  P.  Malanotte-Rizzoli: Adv. Geophys., 33 (1991), 141. https://doi.org/10.1016/S0065-2687(08)60442-2
• 193)   A.  Yamanaka,  Y.  Maeda and  K.  Sasaki: Mater. Des., 165 (2019), 107577. https://doi.org/10.1016/j.matdes.2018.107577
• 194)   A.  Yamamura,  S.  Sakane,  M.  Ohno,  H.  Yasuda and  T.  Takaki: Comput. Mater. Sci., 215 (2022), 111776. https://doi.org/10.1016/j.commatsci.2022.111776
• 195)   A.  Ishii,  A.  Yamanaka,  E.  Miyoshi,  Y.  Okada and  A.  Yamamoto: Model. Sim. Mater. Sci. Eng., 29 (2021), 065012. https://doi.org/10.1088/1361-651X/ac13cd
• 196)   X.  Qiang,  Y.  Iwamoto,  A.  Watanabe,  T.  Kameyama,  X.  He,  T.  Kaneko,  Y.  Shibuta and  T.  Kato: Sci. Rep., 11 (2021), 22285. https://doi.org/10.1038/s41598-021-01666-9
• 197)   Y.  Oka,  M.  Ohno and  K.  Matsuura: Tetsu-to-Hagané, 103 (2017), 755 (in Japanese). https://doi.org/10.2355/tetsutohagane.TETSU-2017-031
• 198)   T.  Oikawa,  Y.  Natsume and  M.  Ohno: ISIJ Int., 62 (2022), 1666. https://doi.org/10.2355/isijinternational.ISIJINT-2022-095
• 199)   Y.  Harashima,  K.  Tamai,  S.  Doi,  M.  Matsumoto,  H.  Akai,  N.  Kawashima,  M.  Ito,  N.  Sakuma,  A.  Kato,  T.  Shoji and  T.  Miyake: Phys. Rev. Mater., 5 (2021), 013806. https://doi.org/10.1103/PhysRevMaterials.5.013806
• 200)   K.  Sakaushi,  A.  Watanabe,  T.  Kumeda and  Y.  Shibuta: ACS Appl. Mater. Interfaces, 14 (2022), 22889. https://doi.org/10.1021/acsami.1c21038
• 201)   S.  Yoshikawa,  R.  Sato,  R.  Akashi,  S.  Todo and  S.  Tsuneyuki: J. Chem. Phys., 157 (2022), 224112. https://doi.org/10.1063/5.0125553
• 202)   Y.  Zhao,  R.  Sato and  S.  Tsuneyuki: J. Non. Cryst. Solids., 600 (2023), 122028. https://doi.org/10.1016/j.jnoncrysol.2022.122028