Micrometer-scale molecular dynamics simulation of microstructure formation linked with multi-phase-field simulation in same space scale

The micrometer-scale polycrystalline microstructure is directly obtained from a 10 billion atom molecular dynamics (MD) simulation of the nucleation and growth of crystals from an undercooled melt, which is performed on a graphics processing unit-rich supercomputer. The grain size distribution in the as-grown microstructure obtained from the MD simulation largely deviates from that resulting from steady-state growth in ideal grain growth, whereas the distribution of the disorientation angle between grains in contact with each other basically agrees with a random distribution. The atomistic configuration of the polycrystalline microstructure is then converted into a phase-field profile (diffuse interface description) of a phase-field model (PFM) and the subsequent grain growth is examined by multi-phase-field (MPF) simulation. A significant achievement in this study is direct mapping of the atomistic configuration into the phase-field profile used in the MPF simulation since only representative parameters for larger-scale model (e.g. interatomic potentials for MD and interfacial parameters for PFM) are extracted from a smaller‐scale simulation in conventional multi-scale modeling. Our new achievement supported by high-performance supercomputing can be regarded as an evolution of multi-scale modeling, which we call inter-scale modeling to differentiate it from conventional multi-scale modeling.

performance computing have enabled the MD simulation of multiple nucleation from an undercooled melt and its subsequent solidification, which forms a polycrystalline microstructure at the submicrometer scale [23][24][25][26]. Moreover, we have proposed a novel technique bridging atomistic and continuum-based simulations by converting an atomic configuration from an MD simulation into the interfacial profile of a PFM [27]. Using this methodology, million-atom MD and MPF simulations of grain growth were carried out from the same initial structure and the results were directly compared at the same spatiotemporal scale [27]. However, it is desirable to employ a much larger MD configuration to discuss grain growth kinetics from a statistical viewpoint. To this end, a 10 billion atom MD simulation of the nucleation and growth of crystals from an undercooled melt is performed on a GPU-rich computer to obtain the polycrystalline microstructure at the micrometer scale in this study. Then, an MPF simulation of the subsequent grain growth is performed using the MDgenerated microstructure consisting of 10 billion atoms. In this study, the iron system is examined as a model case since iron is one of the most practical metals. Finally, a new criterion of multi-scale modeling is discussed on the basis of our simulation results.

Simulation methodology and computational detail
The methodology of the MD simulation basically follows our previous studies [25,26]. The Finnis−Sinclair (FS) potential [28] is employed for the interatomic potential between iron atoms, which is one of the most established potentials for body-centered-cubic (bcc) metals. A leapfrog method [29] is used to integrate a classical equation of motion with a time step of 5.0 fs. A Berendsen thermostat [30] is applied to control the temperature in each step. The Andersen method [31] is applied to independently control the pressure in each direction. Periodic boundary conditions are employed for all computing domain boundaries. Figure 1 shows the initial configuration of the MD simulation prepared as follows. A bcc crystal of iron consisting of 3600×3600×400 unit cells (10 368 000 000 atoms) is annealed at 3000 K for 10 ps with an isobaric-isothermal ensemble. Then, the system is annealed at 2000 K for 200 ps and at 1800 K for 200 ps in stages with the isobaric-isothermal ensemble. In the main calculation, the prepared initial configuration is annealed isothermally for 500 ps under zero pressure with isobaric-isothermal ensemble at 1600 K. Note that our previous simulation using the same interatomic potential [23] confirmed that there is a critical temperature around 1400 K where the nucleation rate becomes maximum. Deviation from the critical temperature exponentially decreases the nucleation rate. Therefore, nucleation seldom The absolute temperature of 1600 K corresponds to a normalized temperature of 0.67T m since the melting point of bcc iron given by the FS potential, T m , is 2400 K [32].
After the main calculation, geometric information of grains is extracted from atom coordinates and visualized following our previous studies [23,25,33]. Firstly, the obtained atomic configurations are analyzed by common neighbor analysis [33,34] to identify solid and liquid structures. In particular, atoms satisfying the following two conditions are identified as belonging to the bcc configuration: (i) 8 neighbor atoms within a cutoff length of 2.75 Å and (ii) 14 neighbor atoms within a cutoff length of 3.4 Å. Next, the disorientation angle relative to the coordination axis is estimated for all atoms with the bcc configuration as follows [25]. The principal axis for each bcc atom is defined as the direction to secondnearest-neighbor atoms, which is the orthogonal direction. The Euler angles between the coordination axis and the principal axis are calculated for all atoms with the bcc configuration. Then, neighboring atoms with a difference in the crystal orientation of less than 3°are classified as belonging to the same grain. The average disorientation angle for all atoms in a grain is regarded as the disorientation angle of the grain. All atoms are colored in accordance with the disorientation angle of the grain in the range from 0°to 62.80°. Details of the methodology are given elsewhere [25].
The MD simulation is performed on the GPU-rich supercomputer TSUBAME3.0 at Tokyo Institute of Technology. An original parallel GPU code [25] is developed with Compute Unified Device Architecture and Message Passing Interface (MPI) for internode communication in the MD simulation. The whole simulation domain is divided into subdomains and each subdomain is assigned to one GPU. The domain decomposition method [29,35] is also applied to accelerate the search for neighboring atoms inside each subdomain. The whole domain is divided into 256 (16×16×1) subdomains, and 256 GPUs (NVIDIA P100) are parallelized for the MD simulation. It takes 78 h for a 500 ps (100 000 MD steps) calculation with approximately 10 billion atoms using 256 GPUs. The postanalysis is also performed on TSUBAME3.0 with an original code developed with MPI. The whole domain is divided into 1792 (16×112×1) subdomains and 1792 CPU cores are parallelized. It takes about 1 h for the postanalysis of the coordinates of all atoms per step. All allocated resources are used exclusively during all computations. The load imbalance is roughly estimated from the difference in the number of atoms allocated in the subdomains, which differs by less than approximately 6.4% in the MD simulation in this study. Figure 2 shows snapshots of the simulation cell during a 500 ps calculation of the nucleation and growth of crystals from an undercooled iron melt. All grains larger than 1.0 nm are colored in accordance with the disorientation angle relative to the coordination axis. Note that solid assemblies in the simulation cell are called grains in this study since it is not conceptually straightforward to distinguish nuclei and grains. Small grains sparsely appear via nucleation in the undercooled melt at 200 ps and they grow larger with time. The solid fraction, which is defined as the ratio of the volume of the solid region to the total volume, is 38.2% at 200 ps. An enlarged view of the snapshot at 250 ps reveals that individual grains are not completely dispersed, with some grains in contact with each other. This was also observed in our previous study, in which heterogeneous nucleation occurred on the surface of previously existing grains [25]. The disorientation angle between the grains in contact is examined later. After that, further nucleation occurs sporadically in the remaining part of the undercooled melt, which is followed by the formation of small grains in the space between  figure 2 shows that the grain size of the microstructure is less than the thickness of the simulation cell. That is, three-dimensional grains are distributed in the simulation cell. The number of grains in the simulation cell at 500 ps is approximately 16 400. This is approximately 10 times as many as that in our previous MD simulation with one billion atoms at the same temperature [25]. That is, the density of grains in the microstructure is almost the same. Therefore, it is considered that the cell size no longer affects the nucleation behavior in the MD simulation when we employ a submicrometer-scale simulation cell. Figure 3(a) shows the grain size distribution normalized by the average grain radius at 200 and 500 ps. The number of grains with a radius of 2.0 nm or more is counted. Here, the grain radius (r) is estimated from the volume of each grain (V ) using the relationship V=4πr 3 /3. The grain size distribution for steady-state growth in ideal grain growth under the conditions of isotropic grain boundary energy and mobility [36,37] with parameters obtained from an ultralarge-scale phase-field simulation [17] is shown for comparison. Note that steady-state growth represents a stage of grain growth in which the grain size distribution normalized by the average grain radius is time-independent. The grain size distribution in the as-grown microstructure via nucleation largely deviates from that resulting from the steadystate growth in ideal grain growth. Especially, the frequency of small grains in the microstructure obtained from the MD simulation is much larger than that in the steady-state growth. Moreover, the frequency of large grains (r/〈r〉 1.5) is also larger than that in the steadystate growth. These trends indicate that the sporadic nucleation [25,38] is dominant rather than spontaneous nucleation. That is, previously existing grains grow larger at the final stage of solidification, which causes the long tail toward large radius in the grain size distribution. This is different from the distribution from the theoretical prediction such as Poisson−Voronoi tessellation [39,40]. Moreover, abundant small grains appear in the sporadic nucleation since there are no enough spaces to grow for subsequent grains. Figure 3(b) shows the distribution of the disorientation angle between grains in contact at 200 and 500 ps. The dashed line represents a random orientation called the Mackenzie distribution [41]. The distribution of the disorientation angle basically agrees with the Mackenzie distribution, which means that grains in the as-grown microstructure via nucleation are randomly oriented in principle. However, a shoulder appears at approximately 60°in the distribution of the disorientation angle. It is known that the existence of abundant coherent twin boundaries with a small grain boundary energy increases the frequency of distribution of disorientation angle at 60° [17,42]. Indeed, some grains in contact with each other constitute twin boundaries in the same manner as in our previous simulation [17]. Interestingly, such twin boundaries appear preferentially at the initial stage of nucleation around 200 ps and the shoulder at 60°d ecreases at 500 ps when solidification is completed. Moreover, a small negative deviation from the Mackenzie distribution appears at small disorientation angles of less than 20°. The grain boundary energy at a low-angle grain boundary decreases with decreasing disorientation angle [32,43,44], which is known as the Read−Shockley relation [45]. Therefore, it is considered that grain rotation also occurs in the process of microstructure formation. However, there are no significant differences between the overall distributions of the grain size and orientation at 200 and 500 ps except for the above points. This implies that the nucleation behavior does not differ significantly between the middle and the end of solidification.

MPF model
The generalized MPF model proposed by Steinbach and Pezzolla [46] is employed for the subsequent MPF simulation of grain growth following our previous study [27]. A polycrystalline system including N grains is represented through N phase-field variables f i (i=1, 2, K, N), which take a value of 1 in the ith grain and 0 in the other grains with a smooth change from 1 to 0 across the grain boundary face. The migration of grain boundaries is reproduced by calculating the time evolution of f i at each spatial point under the constraint of free energy minimization as follows: å å where n is the number of non-zero phase-field variables at the specified point. f M , ij W ij , and a ij represent the phase-field mobility, barrier height, and gradient coefficient of the boundary between ith and jth grains, respectively. These parameters are correlated with the thickness (δ), energy (σ), and mobility (M) of the interface via the following equations: In the practical calculation, we employ δ=6Δx, σ=1 J m −2 , and M=6.34×10 −9 m 4 J −1 s −1 [47]. Δx denotes the mesh size for the finite-difference computation and Δx=1.196 nm is employed. The interfacial energy and mobility employed are from the MD simulation of the grain growth for a pure Fe system in a similar temperature range using the FS potential [47]. Therefore, the effects of the temperature and the material characteristics are reflected via the interfacial parameters, although the effect of anisotropy on the interfacial energy and mobility is not taken into account. The time evolution of equation (1) is numerically solved using the first-order forward difference scheme and the second-order central difference scheme explicitly for time and space, respectively, with a time increment of Δt=0.01 ns. A single GPU computation is employed for the subsequent MPF calculation.

Conversion of atomistic configuration into phase-field profile
Here, the atomistic configuration of the polycrystalline microstructure in the MD simulation at 500 ps (figure 4(a)) is converted into the phase-field profile of the PFM to bridge the MD and MPF simulations. The methodology basically follows previous studies [27]. The MD simulation cell is divided into three-dimensional difference grids of 864×864×96. A periodic boundary condition is employed for all boundaries. After extracting the bcc atoms in each grid, the most dominant grain is assigned to each grid point to obtain a voxel structure. All grains are numbered and are colored gradationally in numerical order ( figure 4(b)). There is no interfacial thickness in the obtained voxel structure but there is a discontinuous change in the phase-field variables. Therefore, relaxation by the MPF simulation without a curvature effect is performed to create diffuse interfaces with an equilibrium profile in the normal direction to the interface while maintaining the grain configuration. The curvature effect is removed by subtracting the contribution of curvature, from the Laplacian term in equation (1). The derivation of the curvature effect in the MPF model is and (c) are colored gradationally in numerical order to examine the changes in morphology due to grain growth later (see figure 5).
shown in [27]. The obtained structure (figure 4(c)) is employed as the initial configuration in the subsequent MPF simulation of grain growth. Figure 5 shows snapshots of the microstructure evolution in the MPF simulation of grain growth. In general, small grains shrink and disappear quickly and large grains grow with time, which is a typical picture of grain growth. Three-dimensional grain growth occurs in the beginning, where the grain size is less than the thickness of the calculation system in principle. However, most of the grains have passed through top and bottom boundaries by approximately 0.3 μs. Then, the growth mode changes from three-dimensional to twodimensional growth. Finally, only two-dimensional columnar grains survive in the polycrystalline microstructure. Figure 6 shows the square of the average grain radius as a function of time during the MPF simulation with two definitions of the grain radius: volume (V )-  equivalent radius r V = [3V/4π)] 1/3 and area (A)-equivalent radius r A =(A/π) 1/2 . Here, r A is calculated for cross-sectional grains appearing on planes in the thickness direction. It is well known that the kinetics of steady-state grain growth follows the parabolic law [48,49] s

MPF simulation of grain growth
where t 0 is the initial time and K is the kinetic coefficient. The brackets represent averages over all grains. There are two stages following the parabolic law, until 0.3 μs and after 0.6 μs, in the graphs in figure 6. These regions correspond to three-and two-dimensional growth, respectively, with a transient period in between. The kinetic coefficient of three-dimensional growth is estimated to be 0.6498 from the slope of the squared average grain size as a function of volume-equivalent radius. In the same manner, the kinetic coefficients of the cross section of the three-dimensional growth and two-dimensional growth are estimated to be 0.4657 and 0.3258, respectively, from the graph of the area-equivalent radius. In our previous study [50], large-scale MPF simulations showed that the kinetic coefficients for three-and two-dimensional steady-state growth are approximately 0.48 and 0.28, respectively, and that of the cross section of three-dimensional growth is approximately 0.33. Therefore, the kinetic coefficients obtained from the present MPF simulation are larger than that of the steady-state growth in ideal grain growth at all periods. The correlation between the three-dimensional and cross-sectional characteristics was discussed in detail in our previous study [50] and therefore we do not go into detail in this study. It is expected from the MPF result that the initial configuration affects the deviation of the kinetic coefficient. Figure 7 shows the grain size distribution normalized by the average grain radius for several time steps in the region of three-dimensional growth. The grain size distribution for the steady-state growth in ideal grain growth [36,37] with parameters from the ultralarge-scale phase-field simulation [17] is shown for comparison. The frequency of small grains is much larger than that in the case of steady-state growth in the initial configuration, which is converted from the MD simulation result. Small grains easily shrink and disappear at the beginning of grain growth. Therefore, a large ratio of small grains in the Figure 7. Grain size distribution normalized by the average grain radius for several time steps in the region of three-dimensional growth. The dashed line represents the grain size distribution for the steady-state growth in ideal grain growth [36,37] with parameters obtained from the ultralarge-scale phase-field simulation [17]. microstructure increases the kinetic coefficient compared with the value under steady-state growth. However, fluctuation appears in the square of the average grain size as a function of time for the region of two-dimensional growth in figure 6 since the number of grains is no longer sufficient to quantitatively discuss the grain growth kinetics at this stage. The grain growth kinetics in a thin-film system will be separately discussed on the basis of another MPF simulation [51].

Conclusions
In this study, a micrometer-scale polycrystalline microstructure is directly obtained from a 10 billion atom MD simulation of the homogeneous nucleation and growth of crystals from an undercooled melt. Such a large-scale simulation is successfully performed on a GPU-rich supercomputer and is the largest MD simulation of a metallurgical process ever to be reported to the best of our knowledge except for several benchmark calculations using the Lennard −Jones potential [52,53]. The atomistic configuration of the polycrystalline microstructure obtained from the MD simulation is successfully converted into the phase-field profile of the PFM using our novel technique. Then, the MPF simulation of grain growth is performed, in which the grain growth kinetics basically agrees with the parabolic law in both three-and two-dimensional grain growth. A significant achievement in this study is to show that the size of large-scale MD simulations has already reached that of phase-field (and other mesoscale) simulations owing to the recent advances in high-performance computing. This will change our criterion of multi-scale modeling. That is, average values of physical properties from a smaller-scale simulation are extracted as parameters and transferred to a larger-scale simulation in conventional multi-scale modeling. For example, parameters of the interatomic potential for MD simulations (i.e. the atomic scale) are derived from ab initio simulations (i.e. the electronic scale), and then, interfacial parameters for PFM simulations (i.e. the microstructure scale) such as the solid-liquid interfacial energy and mobility are derived from MD simulations in turn. On the other hand, the direct mapping of the atomistic configuration into the phase-field profile of the MPF simulation in this study successfully achieved the sharing of information between atomistic and mesoscale simulations, which can be regarded as an evolution of multi-scale modeling. We refer to such modeling as an inter-scale modeling to differentiate it from conventional multi-scale modeling [19]. In the near future, the direct mapping in inter-scale modeling will create further new concepts. For example, an overlap of the space scale enables the 'on-the-fly' use of information from different scales in collaboration with recent techniques from data science. In particular, the data assimilation technique [54] will enable the dynamic prediction of interfacial parameters for MPF simulations from the morphological change in the atomistic configuration in MD simulations.