# Fileset

[Nature_communications_s41467-025-61246-7 (2).pdf](https://mdr.nims.go.jp/filesets/79a314d2-3a98-4549-bfd1-0030d1d9df2c/download)

## Creator

[Takumi Morino](https://orcid.org/0009-0005-0134-7713), [Machiko Ode](https://orcid.org/0000-0002-9500-5466), [Shoichi Hirosawa](https://orcid.org/0000-0002-6572-2552)

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[An explicit integration approach for predicting the microstructures of multicomponent alloys](https://mdr.nims.go.jp/datasets/e0c1141a-af07-45f1-bffa-4a18f2f9ee70)

## Fulltext

An explicit integration approach for predicting the microstructures of multicomponent alloysArticle https://doi.org/10.1038/s41467-025-61246-7An explicit integration approach forpredicting the microstructures ofmulticomponent alloysTakumi Morino 1 , Machiko Ode2 & Shoichi Hirosawa 3Predicting the complex microstructures of practical materials has been alongstanding goal since Gibbs’s pioneering work on predictions for equili-brium of heterogeneous systems. Themost promising approach for achievingthis goal is integrating Calculation of Phase Diagrams (CALPHAD) with phase-field models. This CALPHAD-coupled phase-field model requires two Gibbsfree energy minimisation conditions: equal diffusion potential and internalequilibrium, both grounded in the second law of thermodynamics. However,as implicit functions, these minimisation conditions suffer from the curse ofdimensionality when applied to multicomponent systems, which imposessignificant constraints on simulation capabilities. Here we propose anapproach that incorporates the equal diffusion potential and internal equili-brium conditions into a single explicit function in phase-field equations. Insimulations across various practical materials, our model achieved equal dif-fusion and internal equilibrium conditions. Furthermore, it overcame dimen-sionality limitations, enabling computations for systems with up to 20components. Thus, the proposed approach proves highly versatile and effi-cient, supporting a wide range of practical applications.Many materials widely used for diverse practical applications, such assteels and super-alloys, comprise numerous elements. When devel-oping newmaterials, merely optimising composition is insufficient, asmaterial properties depend significantly on both the average compo-sition and the microscopic texture. This microstructure is char-acterised by the arrangement and size of grains and phases, withdifferences in concentration, crystal structure, and other features. Itsformation process follows the second law of thermodynamics and istheoretically treated as a Gibbs free energy minimisation process.Microstructural control involves twokeymethods: the calculationof phase diagrams (CALPHAD)1,2 and phase-field methods3–5. In theformer method, phase diagrams are calculated to predict equilibriumstates. The Gibbs free energies of the potential phases are assessedusing thermodynamic models tailored to each phase, and the phasefractions and compositions are predicted as a combination thatminimises the total free energy. This paper focuses on twothermodynamic models: the quasi-regular solution model, which isprimarily for liquid and solid solution phases, and the sublatticemodelfor ordered phases, such as intermetallic compounds. In the sublatticemodel, the crystal lattice is divided into sublattices where specificelements are preferentially distributed. The site fractions in eachsublattice are determined to minimise the free energy of the orderedphase, which is referred to as the internal equilibrium2,6 hereinafter.CALPHAD requires two types of equilibrium calculations when a phaseis modelled using the sublattice approach.In contrast, the phase-field model simulates the process ofmicrostructural evolution toward the equilibrium state. The system’sfree energy is represented using a Ginzburg–Landau-type densityfunctional, which incorporates both Gibbs free energy and interphaseand/or grain boundary energy. Then, the primary governing equationis derived by taking the functional derivative. In CALPHAD, the che-mical potentials must be equal across all the equilibrium phases for aReceived: 27 November 2024Accepted: 13 June 2025Check for updates1Yokohama National University, Hodogayaku, Yokohama, Japan. 2National Institute for Materials Science, Tsukuba, Ibaraki, Japan. 3Department of MechanicalEngineering and Materials Science, Yokohama National University, Hodogayaku, Yokohama, Japan. e-mail: morino-takumi-rb@ynu.jpNature Communications |         (2025) 16:6504 11234567890():,;1234567890():,;http://orcid.org/0009-0005-0134-7713http://orcid.org/0009-0005-0134-7713http://orcid.org/0009-0005-0134-7713http://orcid.org/0009-0005-0134-7713http://orcid.org/0009-0005-0134-7713http://orcid.org/0000-0002-6572-2552http://orcid.org/0000-0002-6572-2552http://orcid.org/0000-0002-6572-2552http://orcid.org/0000-0002-6572-2552http://orcid.org/0000-0002-6572-2552http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-61246-7&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-61246-7&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-61246-7&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-61246-7&domain=pdfmailto:morino-takumi-rb@ynu.jpwww.nature.com/naturecommunicationsmultiphase system to satisfy the minimum free-energy state. Thephase-fieldmodel ensures that the diffusion potentials are equal at theinterfaces7, which is equivalent to the local equilibrium concept inthermodynamics.Adopting a Gibbs free energy function from CALPHAD is ideal forpredicting the microstructural evolution of alloys8–12. However, theequal diffusion potential conditionwith theGibbs free energy functionfrom CALPHAD imposes significant computational-time limitationsbecause it becomes an implicit function of the phase concentration,rendering it impractical to solve, even for ternary systems. PerformingCALPHAD-coupled phase-field model calculations in a straightforwardmanner requires repeating the phase-diagram calculations at everycomputational grid point and every timestep, potentially resulting inmore than a billion phase-diagram calculations.To mitigate the computational costs of CALPHAD coupling, var-ious strategies have been introduced, including the parabolicapproximation of free-energy functions13,14, incorporation of machinelearning15,16, and extrapolation of the equal diffusion potentialcondition17–19. However, these approaches face challenges, such asexponentially increasing complexity with additional components orcoarse approximations that compromise accuracy. In the grandpotential approach20,21, the equal diffusion potential condition isautomatically fulfilled; however, it is necessary to express the chemicalpotential explicitly as a function of composition. This involves linear-ising the chemical potential or using parabolic approximations of free-energy functions with respect to composition, which makes the directuse of CALPHAD functions challenging. To address this problem, wedeveloped a phase-field model called the Direct CALPHAD Coupling(DCC)model to efficiently solve the equal diffusionpotential conditionwithout approximation22. However, this model becomes increasinglyunstable as thenumber of components increases and is only applicableto quasi-regular solution-based phases.Despite the importance of ordered phases, relevant research hasbeen minimal owing to the substantial computational time required.Solving the internal equilibrium condition during the phase-fieldsimulation of the superalloy γ‘ phase using the CALPHAD softwarePyCalphad23 is estimated to take 2.3 years (see Supplementary Note 1for details). The need to solve the internal equilibrium condition, inaddition to the equal diffusion potential condition, significantlyincreases the computational time. Thus, there are a few examples inwhich the internal equilibrium condition was properly resolved duringthe phase-field simulation6,24. In most cases, either a parabolicapproximation of the free energy is used25,26 or the internal equilibriumcondition is not solved at all27–29.The phase-field model is widely used for analysing the micro-structures of materials3–5. However, there is no existing model capableof addressing the requirements for the development of new alloys inpractice. The objective of this study is to provide a state-of-the-artmethod for alloy microstructure calculations by developing a phase-field model that satisfies equal diffusion potential and internal equili-brium conditions for multicomponent alloys, including liquid, solidsolutions, and ordered phases, within a realistic computationaltimeframe.Herein, we propose a model that incorporates equal diffusionpotential and internal equilibrium conditions as an explicit function inthe phase-field model. We start by redefining both conditions in thecontext of the phase-field model and then present a formulationincorporating them into the evolution equation of the site fraction tobypass the curse of dimensionality. The proposed model performsrapid and accurate calculations of the solidification of Al, Ni, and Fesystems and the solid-state transformation of the γ‘ phase—a well-known sublattice phase—for an unprecedented number of compo-nents (> 12). This is a generalisedmodel for predictingmicrostructuresin practical materials realised over a century after Gibbs began usingthermodynamics to predict heterogeneous systems30.ResultsDefinition of variables and local minimisation conditionThe equal diffusion potential and internal equilibrium conditions aretreated separately because they are defined in the phase-field andCALPHAD methods, respectively. However, these two conditions canbe expressed by a single formula24. In this section, we first definevariables and then express the equal diffusion potential and internalequilibrium conditions as a unified equation, following Schwen’sformulations24.In the multiphase-field model, the microstructure is representedby the phase-field variable ϕα . This is a function of time and positionthat gives the local fraction of phase α (α = 1, . . . ,N phases). Forinstance,ϕα = 1 indicates the αphase, whereas 0 <ϕα < 1 represents theinterface. The interfacial region, where the phase-field variable chan-ges continuously from 0 to 1, is assumed to be a mixture of multiplephases such thatPNαϕα = 1. At the interface, the phase compositionfulfils the following mixture rule:ci =XNαϕαciα , ð1Þwhere ci represents the overall composition or mixture composition(solvent i=n and solute i= 1, . . . ,n� 1) and ciα represents the com-position of component i of phase α. The relationship between thephase composition and site fraction is defined as follows:ciα =Xa2αSayia, ð2Þwhere Sa is the stoichiometric coefficient representing the ratio of thenumber of atomic positions belonging to sublattice a to the totalnumber of atomic positions in phase α such thatPa2αSa = 1 and yia isthe site fraction of component i of sublattice a in phase α such thatPni yia = 1. a2α implies that the summation is taken only for the phaseα.In the quasi-regular solution model, the crystal lattice is not dividedinto sublattices, i.e., Sa = 1; thus, the site fraction is equal to the phasecomposition. The condition thatminimises the local Gibbs free energyof the microstructure while satisfying Eqs. (1) and (2) is derived usingLagrange’s multiplier24 (see Supplementary Note 2 for details):1Sa∂f α∂yia=1Sb∂f β∂yib, ð3Þwhere fα represents the chemical free energy density of phase α. Theright and left sides of Eq. (3) correspond to the diffusion potential ofeach phase, as follows24:1Sa∂f α∂yia=1∂cia=∂yia∂f α∂yia=∂f α∂cia: ð4ÞTherefore, Eq. (3) represents the equal diffusion potential condi-tion for the interface where α ≠ β. When α = β and α is described by thesublatticemodel, Eq. (3) becomes a condition thatminimises theGibbsfree energy of a phase, as described by the sublattice model, with theconstraint of a fixed phase composition. This is the exact definition ofthe internal equilibrium condition. In this study, Eq. (3), whichencompasses the equal diffusion potential and internal equilibriumconditions, is referred to as the localminimisation condition because itis realised by minimising the local Gibbs free energy of themicrostructure.Explicitly solvable local minimisation conditionIn conventional simulations, the site fraction is treated as a dependentvariable on the composition and phase-field variables. However, thetime evolution of the composition is essentially driven by changes inArticle https://doi.org/10.1038/s41467-025-61246-7Nature Communications |         (2025) 16:6504 2www.nature.com/naturecommunicationsNi (cave = 6.8×101)b Al (cave = 1.1×101) Co (cave = 3.6×100) Cr (cave = 2.4×100)Fe (cave = 7.0×10−1) Hf (cave = 4.7×10−2) Mo (cave = 9.4×10−1) Nb (cave = 6.9×10−1)Si (cave = 1.5×10−1) Ti (cave = 4.6×10−1) Zr (cave = 1.1×10−1)0.5cave2.0caveFe (cave = 6.7×101)a Al (cave = 7.8×10−1) Co (cave = 9.0×10−1) Cr (cave = 1.2×101) Cu (cave = 7.6×10−1)Hf (cave = 2.2×10−2) La (cave = 2.6×10−3) Mn (cave = 8.1×10−1) Mo (cave = 7.0×10−1) Nb (cave = 1.3×10−2)Ni (cave = 9.8×100) P (cave = 4.6×10−2) Pd (cave = 7.1×10−1) S (cave = 1.8×10−4) Si (cave = 8.0×10−1)Ta (cave = 7.9×10−2) Ti (cave = 3.4×10−2) V (cave = 7.1×10−2) W (cave = 7.1×10−1) Y (cave = 2.7×10−3)0.8cave1.2caveFig. 1 | Snapshots of composition distributions at the end of the two simula-tions. a Solidification of an fcc phase in a 20-component Fe system at a solidifi-cation time of 6.0 × 10−2 s. b γ’ phase precipitation and coarsening in a 12-componentNi systematanageing timeof 1.5 × 103 s. The rangeof the colourmaps isbasedon the average composition of each element (cave).Owing to the direct use ofthe CALPHAD database, the partitioning of each element varies according to thephase diagram. The CALPHAD database is obtained from the MatCalc open data-base (https://www.matcalc.at/).Article https://doi.org/10.1038/s41467-025-61246-7Nature Communications |         (2025) 16:6504 3https://www.matcalc.at/www.nature.com/naturecommunicationsthe site fraction, whichmoves toward the localminimisation conditionbecause the composition is expressed by the summation of the sitefractions, as given by Eqs. (1) and (2). Therefore, explicit treatment ofthe site fraction is a more reasonable approach. In this study, we for-mulated the evolution equation of the site fraction in the form of anexplicit function that moves toward the local minimisation conditionand is consistent with the conventional evolution equations of thephase-field variable and composition. The evolution equation pre-sented next is derived in the Supplementary Notes 3 and 4.In the context of the Eulermethod, changes in the site fraction perstep for sublattices a and b are specified as Δyia and Δyib, respectively,and their ratios are denoted as kiba. We define the diffusion potentialdifference function Hi, which corresponds to the difference betweenthe right and left sides of Eq. (3) after changes in the site fraction, asfollows:Hi Δy1a, . . . ,Δyn�1a� �=1Sa∂f α∂yiaΔy1a, . . . ,Δyn�1a� �� 1Sb∂f β∂yibk1baΔy1a, . . . , kn�1ba Δyn�1a� �: ð5ÞWe applymathematical processing to Eq. (5), including expansionto the first order, application of a condition for ensuring numericalstability, and solution of the local minimisation condition (Hi =0). Thisprocess produces the kiba that results in the Δyia and Δyib moving for-ward the local minimisation condition, as follows (see SupplementaryNotes 3 and 4 for details):kba =1Sa∂f α∂yia� 1Sb∂f β∂yib+ ðn� 1ÞΔyia 1Sa∂2f a∂yiað Þ2n� 1ð ÞΔyia 1Sb∂2f β∂yibð Þ2: ð6ÞBy combining Eq. (6) with two necessary conditions for the phase-fieldmodel, including solute diffusion in bulk phases and conservationof solute around the moving interface, we can obtain the Δyia thatmoves forward to the local minimisation condition as follows (seeSupplementary Note 3 for details):Δyia =�PNβ= 1 ϕβ +∂ϕβ∂t Δt� �Pb2βSb1Sa∂f α∂yia� 1Sb∂f β∂yib� �= n�1Sb∂2f β∂yibð Þ2� �+ ∂ci∂t �PNβ = 1∂ϕβ∂tPb2βSbyib� �ΔtPNβ= 1 ϕβ +∂ϕβ∂t Δt� �Pb2βS2bSa∂2f α∂yiað Þ2 =∂2f β∂yibð Þ2:ð7ÞSolving Eq. (7) is straightforwardbecause all variables on the right-hand side are known, making it an explicit function. The first term inthe numerator of Eq. (7) adjusts the site fraction to meet the localminimisation condition. As a result, solving Eq. (7) alone satisfies thelocal minimisation condition without the need for convergence cal-culations, which is expected to significantly reduce the computationalcost. This equation can be reduced analytically to obtain the standardevolution equation for compositions (see Supplementary Note 5 fordetails). Accordingly, an arbitrary evolution equation for compositionscan be employed. In the solidification simulations presented in thefollowing section, an anti-trapping current-infused diffusion equationnecessary for quantitative simulations31 is employed.In the present model, the local minimisation condition and freeboundary problems32, including solute diffusion in the bulk phases,conservation of solute around the moving interface, and theGibbs–Thomson effect, are solved using Eq. (7) and an evolutionequation of the phase-field variable of the standard multiphase-fieldmodel33,34. Furthermore, in the calculation of γ-γ′, the mechanicalequilibrium condition is solved to account for the effect of the elasticfield35.Case studyTwo types of microstructural evolution problems are addressed tovalidate the proposed model. The first problem involves the solidifi-cation of the solid solution phase described by the quasi-regularsolution model, and the second problem involves the solid-statetransformation of the γ′ phase using the sublattice model. Note thatthe proposedmodel directly employs CALPHAD functions without theneed for modifications, such as linearising the chemical potential orapplying parabolic approximations to free-energy functions. Thisapproach allows efficient computation regardless of temporal andspatial temperature changes, including those occurring during pro-cesses such as alloy solidification.The computational process began with (i) setting the initial con-ditions for temperature, composition, and phase fields. The initial sitefraction was input as the equilibrium values corresponding to theseconditions. Subsequently, (ii) the governing equations of the phase-field variable and composition are solved. Following this, (iii) thegoverning equation of site fraction represented by Eq. (7) is addressed.It should be noted that the phase compositions are automaticallydetermined through their relationship with the site fraction, asdescribed in Eq. (2). Finally, (iv) steps (ii) and (iii) are iteratively repe-ated. All CALPHAD data, including those for the quasi-regular solutionmodel and the sublattice model, have been implemented directly inthe phase-field simulation code. For ultra-multicomponent systems,such as those with 20 components, manually hard-coding becomesimpractical due to complexity; therefore, we developed an automatic0 5 10 15 20Number of component10−1010−810−610−4Relative error of diffusion potentialaAl dendriteNi dendriteFe dendriteNi γ-γ'0 5 10 15 20Number of component10−410−310−210−1100Error of diffusion potential (J/mol)bAl dendriteNi dendriteFe dendriteNi γ-γ'0 5 10 15 20Number of component0100200300Calculation time (s)cAl dendriteNi dendriteFe dendriteNi γ-γ'Fig. 2 | Error of the local minimisation condition and calculation time for dif-ferent numbers of components. a Average absolute relative error of the diffusionpotential. b Average absolute discrepancy of the diffusion potential. c Calculationtime. The errors were calculated every 100 timesteps in regions where multiplephases and/or sublattices coexist and subsequently averaged. Blue, orange, andgreenmarkers represent the simulation results of dendritic solidification in Al-, Ni-,and Fe- systems, respectively, as case studies using the quasi-regular solutionmodel, which requires the equal diffusion potential condition. Black markersindicate the results of γ–γ′ solid-state transformation calculated with the sublatticemodel, which requires both the equal diffusion potential condition and the internalequilibrium condition. In panels (a) and (b), error bars show the 95th and 5thpercentiles of the computed values. Source data are provided as a Source Data file.Article https://doi.org/10.1038/s41467-025-61246-7Nature Communications |         (2025) 16:6504 4www.nature.com/naturecommunicationscode generation systemasdescribed in theMethodsection. In the caseof commercial, encrypted databases, the necessary values for calcu-lations can be obtained through the software’s API.Placing a seed of the solid solution phase in the liquid phase andcooling it enabled the growth of dendrites, as shown in Fig. 1a. Whenthe γ′ phase is randomly distributed within the γ phase and subjectedto isothermal aging, it grows as depicted in Fig. 1b. In Fig. 1b, the γ′phase exhibits a cuboidal morphology owing to the contribution ofelastic energy. As these simulations directly utilise the CALPHADdatabasewithout simplifications for computational speed, Fig. 1 showsthepartitioning of eachelement intophases basedonaphasediagram,ensuring accurate and reliable results. Solidification and solid-statetransformation calculations were performed for various numbers ofcomponents. However, the figures for these systems are omitted here,as they exhibit similar patterns of concentration distribution, con-firming the robustness of our computational method across differentalloy compositions.The error in the local minimisation condition was defined as theerror between the right and left sides of Eq. (3), that is, the error in thediffusion potential, at every 100 timesteps in the regions where mul-tiple phases and/or sublattices coexist. Figure 2a depicts the averageabsolute relative error of the diffusion potential. Regardless of thenumber of components, the system, and the CALPHAD model, theerror is of the order 10−5. This is comparable to the precision of thesingle-precision floating-point number type, which is of the order 10−6.Fig. 2b shows the average absolute discrepancy in the diffusionpotential. The discrepancy is of the order 10−1 J/mol under all condi-tions. This value is negligibly small compared to the 102–103 J/moldeviation in the diffusion potential caused by the curvature effect,which was calculated using Equation (A5) from ref. 36. Therefore, itwas confirmed that the local minimisation condition is maintainedthroughout the simulation.Fig. 2c shows the calculation time for different numbers of com-ponents. The calculation time for the solid-state transformation of theγ′ phase is longer than that for the solidification of the solid solutionphase. This is attributed to differences in the thermodynamic models.In the solid solution phase described by the quasi-regular solutionmodel, thermodynamic functions, including ∂f α∂yiaand ∂2f α∂yiað Þ2, in Eq. (7)have to be computed only at the interface to satisfy the equal diffusionpotential condition,whereas in the γ′phase using the sublatticemodel,thermodynamic functions must also be computed in the bulk to fulfilthe internal equilibrium condition. Nevertheless, the calculation of thesublattice model is sufficiently fast to complete a simulation with 12components in just 293 s. The calculation time increased linearly—notexponentially—with an increase in the number of components, indi-cating that our model overcomes the curse of dimensionality. There-fore, a simulation with an unparalleled number of components (20)was completed in 63 s.Our model is superior to conventional methods with regard toboth the number of components and accuracy. ConventionalCALPHAD-coupled phase-field simulations include three componentsutilising machine learning15, three components employing parabolicapproximation37, and ten components using the extrapolationmethod18. However, all these methods use numerical techniques tocircumvent the equal diffusion potential or internal equilibriumconditions and thus inevitably compromise the accuracy of thecomputation. The Newton-based method can be used to solve theconditions directly; however, it is computationally expensive. Theexisting Newton solver for the equal diffusion potential condition—Thermo4PFM—supports only two or three components38. In contrastto these methods, our model changes the governing equations ratherthan relying on numerical techniques. Consequently, calculationsinvolving up to 20 components—significantly exceeding the numberof components for conventional methods—were achieved with highaccuracy and speed without approximation or additional computa-tional parameters.In this study, we developed a model for incorporating equal dif-fusion potential and internal equilibrium conditions into an explicitfunction that solves the problems of solute diffusion in bulk phasesand the conservation of the solute around the moving interface. Theproposed approach can be directly integrated with the CALPHADdatabase without modifying functional forms. The computationalaccuracy and efficiency of our model were evaluated by performingTable 1 | Numerical parametersLiquid-fcc Liquid-fcc Liquid-fcc γ-γ′Al system Ni system Fe system Ni systemDiffusivity in liquid (m2/s) 1.0 × 10−9 1.0 × 10−9 1.0 × 10−9 –Diffusivity in solid (m2/s) 1.0 × 10−13 1.0 × 10−13 1.0 × 10−13 1.0 × 10−16Interface energy (J/m2) 0.2 0.2 0.2 0.02Anti-phase boundary energy (J/m2) – – – 0.05Molar volume (m3/mol) 1.0 × 10−5 1.0 × 10−5 1.0 × 10−5 1.0 × 10−5Strength anisotropy (-) 0.03 0.03 0.03 –Interface mobility (s mol/J) ∞ ∞ ∞ 4.1 × 10−17Elastic constant C11 (GPa) – – – 250.8Elastic constant C12 (GPa) – – – 150.0Elastic constant C44 (GPa) – – – 123.5Lattice misfit between γ and γ′ (-) – – – 0.006Initial temperature (K) 900 1700 1700 1400Cooling rate (K/s) 50 50 50 0Grid resolution (m) 1 × 10−7 1 × 10−7 1 × 10−7 5 × 10−9Interface thickness 4 × 10−7 4 × 10−7 4 × 10−7 2 × 10−8Number of grids (-) 128 × 128 128 × 128 128 × 128 128 × 128Number of timesteps (-) 3 × 104 3 × 104 3 × 104 3 × 104Discrete time width (s) 2 × 10−6 2 × 10−6 2 × 10−6 5 × 10−2Boundary condition (-) Zero Neumann Zero Neumann Zero Neumann PeriodicArticle https://doi.org/10.1038/s41467-025-61246-7Nature Communications |         (2025) 16:6504 5www.nature.com/naturecommunicationsnumerical tests under various conditions, including quasi-regularsolution and sublattice models; Al, Ni, and Fe systems; and 3–20components. Regardless of the CALPHADmodel, system, and numberof components, the local minimisation condition, which is a unifiedexpression of the equal diffusion potential and internal equilibriumconditions, was sufficiently fulfilled. Moreover, an unprecedentedcalculation of a 20-component system was performed in just 1minusing a standard personal computer. The proposedmodel enables thesimulation of the microstructural evolution of multicomponent prac-tical materials without the curse of dimensionality.The proposedmodel is not implementable for certain specific andless frequently used CALPHAD models. Future work will focus onextending our model to encompass all CALPHAD models.MethodsSimulation conditionsThe time evolution equations are discretised using a finite differencemethod: the first-order Euler finite difference scheme is used for timeintegration, and the second-order central finite difference scheme isused for spacediscretisation. All the simulations are conducted using adouble-precisionfloating-point number. Thenumerical parameters foreach calculation system are presented in Table 1. The compositions ofthe liquid-fcc Al system, liquid-fcc Ni system, liquid-fcc Fe system, andγ-γ′ Ni system are presented in Tables 2–5.Governing equations of phase, composition, and elastic fieldsThe following multiphase-field equation is used as the governingequation for the phase-field variable:∂ϕα∂t= � 2NXNβ= 1MαβPNγ = 1W αγ �W βγ� �ϕγ +12 a2αγ � a2βγ� �∇2ϕγh i+ 8πffiffiffiffiffiffiffiffiffiffiffiffiϕαϕβpf α � f β �Pn�1i= 1∂f α∂ciαciα � ciβ� �+ΔGelasticαβ  8>>><>>>:9>>>=>>>;,ð8Þwhere ϕα is the phase-field variable, N represents the number ofphases, Mαβ represents the phase-field mobility of the phase fieldbetween phases α and β,W αβ represents the height of double-obstaclepotential between phases α and β, aαβ is the gradient coefficientbetween phases α and β, f α represents the chemical free energydensity of phase α, ciα represents the phase composition of componenti of phase α, and ΔGelasticαβ represents the elastic driving force. Thedriving force for interface motion due to chemical free energy iscalculated from the distance between the parallel tangents to the free-energy curves of the two phases, assuming that the interface is inquasi-equilibrium7,17. Therefore, the present model is not applicable tocases where the quasi-equilibrium condition is violated, such as inmixed-mode transformations. Mαβ, W αβ, and aαβ are expressed asfollows:Mαβ =π24δmαβð9Þaαβ =ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2δαβσαβqð10ÞW αβ =4σαβδαβ, ð11Þwhere mαβ, δαβ, and σαβ represent the interface mobility, interfacethickness, and interface energy between phases α and β, respectively.According to the standard multiphase-field model, the governingequation of the composition field in Eq. (7) is given as follows17,39,40:∂ci∂t=∇ �XNα = 1ϕαDiα∇ciα , ð12Þwhere Diα is the diffusion coefficient of composition i in phase α.Solidification is calculated by introducing the anti-trapping currentrequired for quantitative simulation41 to the governing equation of thecomposition field:∂ci∂t=∇ �XNα = 1ϕαDiα∇ciα +∇ �XNα >βXNβ= 1Jiαβ ð13ÞJiαβ =aαβffiffiffiffiffiffiffiffiffiffiffiffi2Wαβq ciα � ciβ� � ffiffiffiffiffiffiffiffiffiffiffiffiϕαϕβq ∂ϕα∂t∇ ϕα � ϕβ� �∇ ϕα � ϕβ� ���� ��� , ð14Þwhere Jiαβ represents the anti-trapping current of component i. Thisequation is derived assuming no diffusion in phase β. When the anti-trapping current is used, the relationship between the phase field andTable 2 | Compositions of the liquid-fcc Al systemSystem Composition (mol%)Ternary Al-0.5Cu-6.3MgSenary Al-3.8Cu-0.5Mg-0.5Cr-0.5Fe-0.5MnNonary Al-2.8Cu-0.5Mg-0.5Cr-0.5Fe-0.5Mn-0.5Ni-0.5Sc-0.5SiDuodecimal Al-3.6Cu-0.5Mg-0.5Cr-0.5Fe-0.5Mn-0.5Ni-0.5Sc-0.5Si-0.5Ti-0.5Zn-0.5ZrTable 4 | Compositions of the liquid-fcc Fe systemSystem Composition (mol%)Ternary Fe-22.6Cr-10NiSenary Fe-26Cr-10Ni-1Al-1Co-1CuNonary Fe-23.5Cr-10Ni-1Al-1Co-1Cu-0.1Hf-0.1La-1MnDuodecimal Fe-21Cr-10Ni-1Al-1Co-1Cu-0.1Hf-0.1La-1Mn-1Mo-0.1Nb-0.1PQuindecimal Fe-16.5Cr-10Ni-1Al-1Co-1Cu-0.1Hf-0.1La-1Mn-1Mo-0.1Nb-0.1P-1Pd-0.1S-1SiOctodecimal Fe-16Cr-10Ni-1Al-1Co-1Cu-0.1Hf-0.1La-1Mn-1Mo-0.1Nb-0.1P-1Pd-0.1S-1Si-0.1Ta-0.1Ti-0.1VVigesimal Fe-13.5Cr-10Ni-1Al-1Co-1Cu-0.1Hf-0.1La-1Mn-1Mo-0.1Nb-0.1P-1Pd-0.1S-1Si-0.1Ta-0.1Ti-0.1V-1W-0.1YTable 3 | Compositions of the liquid-fcc Ni systemSystem Composition (mol%)Ternary Ni-16.3Al-1CoSenary Ni-13.6Al-1Co-1Cr-0.3Cu-1FeNonary Ni-11.6Al-1Co-1Cr-0.3Cu-1Fe-0.3Hf-0.3Mn-1MoDuodecimal Ni-5.2Al-1Co-1Cr-0.3Cu-1Fe-0.3Hf-0.3Mn-1Mo-1Nb-1Si-1TiQuindecimal Ni-4.5Al-1Co-1Cr-0.3Cu-1Fe-0.3Hf-0.3Mn-1Mo-1Nb-1Si-1Ti-0.3V-1W-0.3ZrTable 5 | Compositions of the γ-γ′ Ni systemSystem Composition (mol%)Ternary Ni-18.9Al-5CoSenary Ni-16.9Al-5Co-5Cr-1Fe-0.3HfNonary Ni-17Al-5Co-5Cr-1Fe-0.3Hf-1Mo-1Nb-1SiDuodecimal Ni-14.5Al-5Co-5Cr-1Fe-0.3Hf-1Mo-1Nb-1Si-1Ti-1W-0.3ZrArticle https://doi.org/10.1038/s41467-025-61246-7Nature Communications |         (2025) 16:6504 6www.nature.com/naturecommunicationsinterface mobilities is derived under the thin-interface limitcondition42:1mαβ=σαβMαβa2αβ� π8aαβffiffiffiffiffiffiffiffiffiffiffiffi2W αβq ζαβ, ð15Þwhere dijα is the element of the inverse matrix of the interdiffusionmatrix. For calculating the solid-state transformation, the elasticdriving force for the evolution of the phase-field variable is expressedasΔGelasticαβ = � Cijklεelklε0ij , ð16Þwhere Cijkl is the effective elastic modulus tensor, εelij represents theelastic strain, and ε0ij represents the eigenstrain. The effective elasticmodulus is assumed to be uniform for the γ and γ′ phases; that is, theelastic homogeneous case is considered for simplicity. The γ and γ′phases with four different variants are represented by five phase-fieldvariables. The eigenstrain is assumed to be proportional to the localvolume fraction of the γ′ phase:ε0ij = ε0δij 1� ϕγ� �, ð17Þwhere ε0 represents the lattice misfit between phases γ and γ′, δklrepresents the Kronecker delta function, and ϕγ is the phase-fieldvariable of phase γ. The elastic strain is expressed asεelkl = �εcij + δεij + ε0ij , ð18Þwhere �εcij represents the homogeneous strain and δεij represents theheterogeneous strain. The homogeneous strain is calculated by aver-aging the eigenstrain over the volume V as follows:�εcij =1VZVε0ij dV : ð19ÞThe heterogeneous strain is calculated from the mechanicalequilibrium condition as follows:Cijkl∂2uk∂xj∂xl=Cijkl∂ε0kl∂xj, ð20Þwhere uk represents the local displacement.Computing environmentThe simulation code is developed using the Taichi programming lan-guage—a domain-specific language integrated within Python43—toincrease the computational efficiency. The Fast Fourier Transform inNumpy is used to solve the mechanical equilibrium equation. ThePython version used in this study is 3.11, with Taichi version 1.7.2 andNumPy version 1.26.4. All figures were generated using matplotlibversion 3.8.4. Simulations are conducted using a MacBook Pro laptop,Sonoma 14.5 OS, and an 11-core CPU with an Apple M3 Pro chip and 18GB of unified memory.CALPHAD coupling methodA thermodynamic database is obtained from the MatCalc open data-base (https://www.matcalc.at/). This database is provided in the TDBformat and had to be converted for utilisation in the phase-fieldsimulations. The Python function for the chemical free energy is cre-ated using Python’s regular expression operation library, Re. Thecontribution of the magnetic effects to the chemical free energy isneglected because it was very small, except at the Curie point. The firstand second derivatives of the free energy is calculated using Taichi’sdifferentiable programming and a Python function developed usingPython’s analytical differentiation library Sympy, respectively.Although the analytical differentiation of Sympy can be used todevelop a Python function for the first derivative of the free energy,differentiable programming is used because the CALPHAD functionsbecome very long when differentiated once, making the Taichi com-pilation time-consuming. However, we use Sympy for the secondderivative because CALPHAD functions become shorter when differ-entiated twice.The system for converting TDB files into Python functions hasbeen applied to TDB files released by the National Institute for Mate-rials Science. This is called the Python module project. The convertedPython functions are available from CPDDB.Reporting summaryFurther information on research design is available in the NaturePortfolio Reporting Summary linked to this article.Data availabilityThe sourcedata in the formof a Pythonfile are also available onGitHub(https://github.com/takumimorino/phase-field). Source data are pro-vided in this paper.Code availabilityAll the simulation codes, including the phase-field simulation, materialparameters, and CALPHAD coupling (TDB conversion) codes, can bedownloaded from GitHub (https://github.com/takumimorino/phase-field). The simulation codes are also accessible via Code Ocean44(https://codeocean.com/capsule/3744764/tree/v1). All codes arereleased under the MIT license.References1. Lukas, H., Fries, S. G. & Sundman, B. Computational Thermo-dynamics: The Calphad Method (Cambridge Univ., 2007).2. Hillert, M. Phase Equilibria, Phase Diagrams and Phase Transfor-mations: Their Thermodynamic Basis. Cambridge Univ: NewYork, 2008.3. Boettinger, W. J., Warren, J. A., Beckermann, C. & Karma, A. Phase-field simulation of solidification. Annu. Rev. Mater. Res. 32,163–194 (2002).4. Steinbach, I. Phase-fieldmodels inmaterials science.Modell. Simul.Mater. Sci. Eng. 17, 073001 (2009).5. Steinbach, I. & Salama, H. Lectures on Phase Field (Springer NatureSwitzerland, Cham, 2023).6. Zhang, L., Stratmann, M., Du, Y., Sundman, B. & Steinbach, I.Incorporating the CALPHAD sublattice approach of ordering intothe phase-field model with finite interface dissipation. Acta Mater.88, 156–169 (2015).7. Kim, S.G., Kim,W. T. &Suzuki, T. Phase-fieldmodel for binary alloys.60, 7186–7197 (1999).8. Grafe, U., Böttger, B., Tiaden, J. & Fries, S. G. Coupling of multi-component thermodynamic databases to a phase field model:Application to solidification and solid state transformations ofsuperalloys. Scr. Mater. 42, 1179–1186 (2000).9. Kobayashi, H., Ode, M., Gyoon Kim, S., Tae Kim, W. & Suzuki, T.Phase-field model for solidification of ternary alloys coupled withthermodynamic database. Scr. Mater. 48, 689–694 (2003).10. Steinbach, I., Böttger, B., Eiken, J., Warnken, N. & Fries, S. G. CAL-PHAD and phase-field modeling: A successful liaison. J. PhaseEquilib. Diffus. 28, 101–106 (2007).11. Nomoto, S., Kusano,M., Kitano, H. &Watanabe,M.Multi-phase fieldmethod for solidification microstructure evolution for a Ni-basedalloy in wire arc additive manufacturing.Metals 12, 1720 (2022).12. Uddagiri, M., Shchyglo, O., Steinbach, I. & Tegeler, M. Solidificationof the Ni-based superalloy CMSX-4 simulated with full complexityin 3-dimensions. Prog. Addit. Manuf. 9, 1185–1196 (2024).Article https://doi.org/10.1038/s41467-025-61246-7Nature Communications |         (2025) 16:6504 7https://www.matcalc.at/https://cpddb.nims.go.jp/https://github.com/takumimorino/phase-fieldhttps://github.com/takumimorino/phase-fieldhttps://github.com/takumimorino/phase-fieldhttps://codeocean.com/capsule/3744764/tree/v1www.nature.com/naturecommunications13. Yang, C., Wang, X., Wang, J. & Huang, H. Multiphase-field approachwith parabolic approximation scheme. Comput. Mater. Sci. 172,10–9322 (2020).14. Yang, C., Wang, J., Xing, H. & Huang, H. A parabolic approximationscheme for multi-phase-filed simulation of non-isothermal solidifi-cation. Mater. Today Commun. 28, https://doi.org/10.1016/j.mtcomm.2021.102712 (2021).15. Jiang, X., Zhang, R., Zhang, C., Yin, H. &Qu, X. Fast prediction of thequasi phase equilibrium in phase field model for multicomponentalloys based onmachine learningmethod.Calphad 66, https://doi.org/10.1016/j.calphad.2019.101644 (2019).16. Ren, L., Geng, S., Jiang, P., Gao, S. & Han, C. Numerical simulationof dendritic growth during solidification process usingmultiphase-field model aided with machine learning method.Calphad 78, https://doi.org/10.1016/j.calphad.2022.102450 (2022).17. Eiken, J., Böttger, B. & Steinbach, I. Multiphase-field approach formulticomponent alloys with extrapolation scheme for numericalapplication. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 73,066122 (2006).18. Böttger, B., Eiken, J. & Apel, M. Multi-ternary extrapolation schemefor efficient coupling of thermodynamic data to a multi-phase-fieldmodel. Comput. Mater. Sci. 108, 283–292 (2015).19. Yang, C., Xu, Q. & Liu, B. A high precision extrapolation method inmultiphase-field model for simulating dendrite growth. J. Cryst.Growth 490, 25–34 (2018).20. Plapp, M. Unified derivation of phase-field models for alloy solidi-fication from a grand-potential functional. Phys. Rev. E 84,031601 (2011).21. Choudhury, A. & Nestler, B. Grand-potential formulation for multi-component phase transformations combined with thin-interfaceasymptotics of the double-obstacle potential. Phys. Rev. E 85,021602 (2012).22. Morino, T., Ode,M. &Hirosawa, S. Direct CALPHAD coupling phase-field model: Closed-form expression for interface compositionsatisfying equal diffusion potential condition. Phys. Rev. E 109,065303 (2024).23. Otis, R. & Liu, Z.-K. Pycalphad: CALPHAD-based ComputationalThermodynamics in Python. J. Open Res. Softw. 5, https://doi.org/10.5334/jors.140 (2017).24. Schwen, D., Jiang, C. & Aagesen, L. K. A sublattice phase-fieldmodel for direct CALPHAD database coupling. Comput. Mater. Sci.195, 110466 (2021).25. Zhu, J. Z. et al. Three-dimensional phase-field simulations of coar-sening kinetics of γ′ particles in binary Ni–Al alloys. Acta Mater. 52,2837–2845 (2004).26. Osada, T. et al. Virtual heat treatment for γ-γ′ two-phase Ni-Al alloyon the materials Integration system, MInt. Mater. Des. 226,111631 (2023).27. Zhu, J. Z., Liu, Z. K., Vaithyanathan, V. & Chen, L. Q. Linking phase-fieldmodel toCALPHAD:Application toprecipitate shapeevolutionin Ni-base alloys. Scr. Mater. 46, 401–406 (2002).28. Wang, J. C., Osawa, M., Yokokawa, T., Harada, H. & Enomoto, M.Modeling the microstructural evolution of Ni-base superalloys byphase field method combined with CALPHAD and CVM. Comput.Mater. Sci. 39, 871–879 (2007).29. Kitashima, T. & Harada, H. A new phase-fieldmethod for simulatingγ′ precipitation in multicomponent nickel-base superalloys. ActaMater. 57, 2020–2028 (2009).30. Gibbs, J. W. On the equilibrium of heterogeneous substances. Am.J. Sci. s3–16, 441–458 (1878).31. Karma, A. Phase-field formulation for quantitativemodeling of alloysolidification. Phys. Rev. Lett. 87, 115701 (2001).32. Fasano, A. & Primicerio, M. (ed.). Free Boundary Problems: Theoryand Applications (Pitman, Boston, 1983).33. Steinbach, I. et al. A phase field concept for multiphase systems.Phys. Nonlinear Phenom. 94, 135–147 (1996).34. Tiaden, J., Nestler, B., Diepers, H. J. & Steinbach, I. The multiphase-field model with an integrated concept for modelling solute diffu-sion. Phys. Nonlinear Phenom. 115, 73–86 (1998).35. Chen, L.-Q. Phase-field models for microstructure evolution. Annu.Rev. Mater. Res. 32, 113–140 (2002).36. Kim, S. G. et al. Phase-field model with relaxation of the partitioncoefficient. Comput. Mater. Sci. 188, 110184 (2021).37. Yenusah, C. O. et al. Three-dimensional Phase-field simulation of γ″precipitation kinetics in Inconel 625 during heat treatment. Com-put. Mater. Sci. 187, 110123 (2021).38. Fattebert, J.-L., DeWitt, S., Perron, A. & Turner, J. Thermo4PFM:Facilitating Phase-field simulations of alloys with thermodynamicdriving forces. Comput. Phys. Commun. 288, https://doi.org/10.1016/j.cpc.2023.108739 (2023).39. Kim, S., Kim, W., Suzuki, T. & Ode, M. Phase-field modeling ofeutectic solidification. J. Cryst. Growth 261, 135–158 (2004).40. Bottger, B., Eiken, J. & Steinbach, I. Phase field simulation ofequiaxed solidification in technical alloys. Acta Mater. 54,2697–2704 (2006).41. Kim, S. G. Phase-field model with anti-trapping current for multi-component alloyswith arbitrary thermodynamic properties. ActaMater. 55, 4391–4399 (2007).42. Karma, A. & Rappel, W. J. Quantitative phase-field modeling ofdendritic growth in two and three dimensions. Phys. Rev. E 57,4323–4349 (1998).43. Hu, Y., Li, T. M., Anderson, L., Ragan-Kelley, J. & Durand, F. Taichi: Alanguage for high-performance computation on spatially sparsedata structures. ACM Trans. Graph. 38, 1–16 (2019).44. Morino, T., Ode, M. & Hirosawa, S., An explicit integration approachfor predicting the microstructures of multicomponent alloys, CodeOcean, https://doi.org/10.24433/CO.9466805.v1 (2025).AcknowledgementsWe would like to thank Niterra Materials Co., Ltd. (formerly ToshibaMaterials Co., Ltd., Kanagawa, Japan) and JSPS KAKENHI Grant No.22K04794 for financial support and Editage (www.editage.jp.) for Eng-lish language editing.Author contributionsT.M. designed the research andperformedcodingandcalculations. T.M.and M.O. wrote the paper. T.M., M.O. and S.H. contributed to the dis-cussion and revision of this paper.Competing interestsThe authors declare no competing interests.Additional informationSupplementary information The online version containssupplementary material available athttps://doi.org/10.1038/s41467-025-61246-7.Correspondence and requests for materials should be addressed toTakumi Morino.Peer review information Nature Communications thanks ChaitanyaBhave, Daniel Schwen and Ingo Steinbach for their contribution to thepeer review of this work. A peer review file is available.Reprints and permissions information is available athttp://www.nature.com/reprintsPublisher’s note Springer Nature remains neutral with regard to jur-isdictional claims in published maps and institutional affiliations.Article https://doi.org/10.1038/s41467-025-61246-7Nature Communications |         (2025) 16:6504 8https://doi.org/10.1016/j.mtcomm.2021.102712https://doi.org/10.1016/j.mtcomm.2021.102712https://doi.org/10.1016/j.calphad.2019.101644https://doi.org/10.1016/j.calphad.2019.101644https://doi.org/10.1016/j.calphad.2022.102450https://doi.org/10.1016/j.calphad.2022.102450https://doi.org/10.5334/jors.140https://doi.org/10.5334/jors.140https://doi.org/10.1016/j.cpc.2023.108739https://doi.org/10.1016/j.cpc.2023.108739https://doi.org/10.24433/CO.9466805.v1http://www.editage.jphttps://doi.org/10.1038/s41467-025-61246-7http://www.nature.com/reprintswww.nature.com/naturecommunicationsOpen Access This article is licensed under a Creative CommonsAttribution 4.0 International License, which permits use, sharing,adaptation, distribution and reproduction in any medium or format, aslong as you give appropriate credit to the original author(s) and thesource, provide a link to the Creative Commons licence, and indicate ifchanges were made. The images or other third party material in thisarticle are included in the article’s Creative Commons licence, unlessindicated otherwise in a credit line to the material. If material is notincluded in the article’s Creative Commons licence and your intendeduse is not permitted by statutory regulation or exceeds the permitteduse, you will need to obtain permission directly from the copyrightholder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.© The Author(s) 2025Article https://doi.org/10.1038/s41467-025-61246-7Nature Communications |         (2025) 16:6504 9http://creativecommons.org/licenses/by/4.0/http://creativecommons.org/licenses/by/4.0/www.nature.com/naturecommunications An explicit integration approach for predicting the microstructures of multicomponent alloys Results Definition of variables and local minimisation condition Explicitly solvable local minimisation condition Case study Methods Simulation conditions Governing equations of phase, composition, and elastic fields Computing environment CALPHAD coupling method Reporting summary Data availability Code availability References Acknowledgements Author contributions Competing interests Additional information