2021 Volume 62 Issue 1 Pages 1-15
Revealing atomic-scale distributions of energy and stress in defective or complex systems, based on the behavior of electrons, should contribute much to materials science and engineering, while only few practical ab initio methods were developed for this purpose. Thus, we developed computational techniques of local-energy and local-stress calculations within the plane-wave PAW (projector augmented wave)-GGA (generalized gradient approximation) framework. This is natural extension of ab initio energy-density and stress-density schemes, while the inherent gauge dependency is removed by integrating these densities in each local region where the contained gauge-dependent terms are integrated to be zero. In this overview, we explain our scheme with some details and discuss the concepts or physical meanings of local energy and local stress via the comparison with related schemes using similar densities, LCAO (linear combination of atomic orbitals) methods, Green’s function formulation implemented by multiple scattering or TB (tight-binding) methods, or EAM (embedded-atom method) potentials. We present recent applications to metallic surfaces, grain boundaries (GBs) with and without solute segregation, tensile tests of metallic GBs, local elastic constants of microstructures in alloys, and machine-learning based GB-energy prediction, where the local-energy and local-stress analyses provide novel aspects of phenomena, deep insights into the mechanism, and effective data for novel machine-learning based modelling. We discuss unsettled issues and future applications, especially for large-scale metallic systems.
By virtue of the development of density-functional theory (DFT)1,2) implemented with LDA (local density approximation) or GGA (generalized gradient approximation),3) and the development of efficient methods for electronic-structure calculations of large supercell systems4–6) in the plane-wave pseudopotential framework,7–10) we can perform ab initio calculations of crystals, defects, surfaces, interfaces, grain boundaries (GBs), alloys and so on, where stable structures and various properties can be clarified via the behaviors of electrons and atoms. Such ab initio calculations with high accuracy and predicting power are contributing much to materials science and engineering. About the mechanical properties of materials, for example, recent review articles highlight such ab initio calculations.11–14) The development of package software for the plane-wave DFT method such as VASP, CASTEP, ABINIT and so on has also made great contributions in this research area. By such plane-wave DFT codes, we obtain a total energy and a stress tensor15,16) as quantities integrated or averaged in a supercell, while we cannot obtain any local distribution of energy or stress inside a supercell. If local distributions of energy and stress in a large supercell containing a defect were obtained, understanding on the nature of a defect should be greatly improved, and the validity of a supercell size would be clearly judged by examining the energy and stress in bulk-like regions in a supercell, while only few practical ab initio methods for this purpose were developed. Thus, we have developed practical computational techniques to obtain local energy and local stress inside a supercell17,18) in the plane-wave PAW (projector augmented wave)9,10)-GGA method as the latest plane-wave DFT framework, which were coded in the package software QMAS (Quantum Materials Simulator),19) developed in AIST (National Institute of Advanced Industrial Science and Technology), Japan.
The energy density and stress density, previously proposed in the plane-wave DFT method,20,21) have been reformulated in the PAW-GGA framework.17,19) These densities are calculated by using occupied wave functions and self-consistent charge distribution, and contain distribution information of energy and stress in a supercell, while containing gauge-dependent terms.20,21) Thus, the density values or their integrated values in arbitrary regions cannot be regarded as unique physical quantities, although the integration of these densities in the full supercell leads to the conventional total energy and stress tensor by the plane-wave DFT method. In order to settle this issue, we integrate the energy and stress densities in each partitioned local region where the gauge-dependent terms are integrated to be zero, leading to the local energy and local stress of atoms, layers or clusters as unique physical quantities.17,18) Our scheme has been applied to various systems or phenomena such as metallic surfaces,17,22,23) GBs with and without solute segregation,18,24–31) tensile tests of metallic GBs,32,33) local elastic constants of microstructures in alloys,34) and combination with machine-learning based modelling.35,36) These calculations have revealed novel aspects of materials phenomena, and provided novel insights into the mechanism and effective data for novel machine-learning based modelling. In this overview, we explain our scheme and recent applications with some details. About the local energy, a quite similar scheme within the PAW-GGA framework was also developed by Trinkle and co-workers.37) We explain the differences.
The concepts of local energy and local stress are quite important in materials science and engineering, while understanding or investigation is not enough. By using classical interatomic potentials such as EAM (embedded-atom method),38,39) we can obtain atomistic distributions of local energy and local stress. However, the definition of such local quantities is not straightforward. There is conceptual difficulty in defining such atomistic quantities unambiguously, because of intrinsic extension of wave functions, continuous electronic charge distribution, and long-range electrostatic interactions in solids. On the other hand, the quantum-mechanical stress-tensor field was discussed in the early days of quantum mechanics,40,41) and related energy or stress densities (fields) were proposed for solids,15,20,21) while there exist unsolved ambiguities, concerning uniqueness or observability. Otherwise, electronic-structure calculations of defects or surfaces in solids from local viewpoints, using LCAO (linear combination of atomic orbitals) methods or Green’s function formulation implemented by multiple scattering (KKR (Korringa-Kohn-Rostoker) type) or tight-binding (TB) methods,42–45) have a long history, and have been used to obtain related local quantities. In this overview, we discuss the history and physical meanings of local energy and local stress via the comparison with various related schemes mentioned above. About the terminology in this article, local energy and local stress mean quantities defined in some definite local spaces as local contributions to integrated or averaged quantities of the whole system.
Atomic-level stresses at surfaces, defects or amorphous structures are quite different from conventional elastic stresses by elastic strain fields, in the nature and origin. For a surface stress, surface atoms with less coordination should prefer different in-plane interatomic distances from those in the bulk, while the in-plane periodicity is restricted to the bulk one, resulting in an excess surface stress. In more detail, electron redistribution or orbital rehybridization at the surface and near-surface layers should induce excess-stress distribution there.46) Thus, the surface stress could occur even without any local strains. Atomic-level stresses in amorphous structures occur due to local deviations of coordination numbers, interatomic distances or atomic volumes, as first clarified by Egami and Vitek47–49) using classical inter-atomic potentials. Similar atomic-level stresses were also observed at core regions of metallic GBs by our ab initio local-stress scheme.18,24,25) GB-core structures have features common to amorphous or dislocation-core structures, consisting of peculiar atomic rings or clusters. Recently, Egami’s group also developed the ab initio local-stress scheme based on the Green’s function formalism,50) on which we discuss the comparison in this overview article. In contrast, a conventional elastic stress occurs by an elastic strain field in a crystal, consisting of continuous displacements of successive atomic arrays with long wave lengths. Such elastic strain and stress vary quite slowly at the atomic scale as described by a continuum model, while atomic-level stresses vary rapidly from negative to positive even for neighboring atoms, which cannot be dealt with by a continuum model. This is related to the difference between the microscopic and macroscopic stresses. As will be discussed, our scheme as well as classical inter-atomic potentials can deal with both atomic-level and conventional elastic stresses.
In this overview, we explain our local-energy and local-stress scheme17,18) with some details in Sec. 2, which supplements our previous publications. Each term in the energy and stress densities in the PAW-GGA framework is explained. The conditions for gauge-independent local regions are specified with examples of partitioning methods. In Sec. 3, we discuss the history and physical meanings of local energy and local stress via the comparison with similar or related schemes, recently37,50) or previously proposed, using similar or alternative energy or stress densities, LCAO methods, Green’s function formulation, or EAM potentials. In Sec. 4, our recent applications to various systems or phenomena are presented, including the combination with machine-learning techniques. In Sec. 5, we discuss unsettled issues in theoretical or technical aspects and future applications.
We obtain local energies and local stresses by integrating an energy density $\varepsilon _{\textit{tot}}(\vec{r})$20) and a stress density $\tau _{\alpha \beta }(\vec{r})$21) within respective local regions partitioned in a supercell. The energy and stress densities are defined as supercell-periodic real-space integrant functions, of which the integrations throughout the supercell correspond to the conventional total energy Etot and stress tensor σαβ by the plane-wave DFT method15,16) as follows:
| \begin{equation} E_{\textit{tot}} = \int_{\Omega}\varepsilon_{\textit{tot}}(\vec{r})d\vec{r} \end{equation} | (1) |
| \begin{equation} \sigma_{\alpha \beta} = \frac{1}{\Omega}\frac{\partial E_{\textit{tot}}}{\partial \varepsilon_{\alpha \beta}} = \frac{1}{\Omega}\int_{\Omega}\tau_{\alpha \beta}(\vec{r})d\vec{r}, \end{equation} | (2) |
These densities are intrinsically not unique, but gauge-dependent, because of the presence of gauge-dependent terms.20,21) Thus, density values or their integrated values in arbitrary regions cannot be regarded as unique physical quantities. In our scheme,17,18) we define partitioned local regions Ωi where the contained gauge-dependent terms are integrated to be zero, and the densities are integrated in such local regions, leading to gauge-independent local energies and local stresses as follows:
| \begin{equation} E_{\textit{tot}}^{i} = \int_{\Omega_{i}}\varepsilon_{\textit{tot}}(\vec{r})d\vec{r} \end{equation} | (3) |
| \begin{equation} \sigma_{\alpha \beta}^{i} = \frac{1}{\Omega_{i}}\int_{\Omega_{i}}\tau_{\alpha \beta}(\vec{r})d\vec{r}. \end{equation} | (4) |
| \begin{equation} E_{\textit{tot}} = \sum\nolimits_{i}E_{\textit{tot}}^{i} \end{equation} | (5) |
| \begin{equation} \sigma_{\alpha \beta} = \frac{1}{\varOmega} \sum\nolimits_{i}\varOmega_{i}\sigma_{\alpha \beta}^{i}. \end{equation} | (6) |
The energy and stress densities,20,21) originally developed within the NCPP (norm-conserving pseudopotential method7))-LDA and USPP (ultra-soft pseudopotential method8))-LDA frameworks, respectively, were reformulated in the PAW-GGA framework.17,19) The total energy in the PAW method is expressed as a sum of several terms,10) and the total-energy density $\varepsilon _{\textit{tot}}(\vec{r})$ is given as a sum of corresponding density terms as follows:
| \begin{align} \varepsilon_{\textit{tot}}(\vec{r}) &= \varepsilon_{\textit{kin}}(\vec{r}) + \varepsilon_{\textit{xc}}(\vec{r}) + \varepsilon_{H}(\vec{r}) + \varepsilon_{\textit{local}}(\vec{r}) \\ &\quad + \varepsilon_{\textit{non-local}}(\vec{r}) + \varepsilon_{\textit{ovrl-self}}(\vec{r}). \end{align} | (7) |
In eq. (7), the first four are expressed as usual density forms, while the last two terms are given as a sum of quantities with respect to each atomic site $\overrightarrow{R_{I}}$, spatially distributed by a delta function $\delta (\vec{r} - \overrightarrow{R_{I}})$, which is broadened by a Gaussian so as to have an usual density form.20) $\varepsilon _{\textit{local}}(\vec{r})$ is non-zero only inside each ionic radius,10) and the distribution itself has no physical meaning, similarly to the two terms, $\varepsilon _{\textit{non-local}}(\vec{r})$ and $\varepsilon _{\textit{ovrl-self}}(\vec{r})$. For these three terms, only the integrations inside each ionic radius are meaningful. Similar problems exist in the stress density. Thus, the whole region inside each ionic radius at each atomic site should be included in a local region for the integration of energy and stress densities. Otherwise, such ionic cores can be divided by symmetric planes in symmetric configurations.20,51)
As mentioned in Sec. 1, Trinkle and co-workers37) developed a quite similar scheme of local-energy calculations within the PAW-GGA framework. Their total-energy density consists of similar terms, while $\varepsilon _{\textit{local}}(\vec{r})$ is not included, because local pseudopotentials are constructed so as to remove it.
Details of the terms in eq. (7) are given in Ref. 17). We explain the terms concerning the gauge-dependent problem. The kinetic-energy density has two forms, asymmetric and symmetric ones as follows:
| \begin{align} \varepsilon_{\textit{kin}}(\vec{r}) & = - \frac{1}{2}\sum_{i}^{\textit{occ}}f_{i}\psi_{i}^{*}(\vec{r})\nabla^{2}\psi_{i}(\vec{r})\\ & = \frac{1}{2}\sum\nolimits_{i}^{\textit{occ}}f_{i}\nabla \psi_{i}^{*}(\vec{r}) \cdot \nabla \psi_{i}(\vec{r}) - \frac{1}{4}\nabla^{2}\rho_{e}(\vec{r}), \end{align} | (8) |
The electrostatic-energy density is expressed by the Maxwell form,20)
| \begin{equation} \varepsilon_{H}(\vec{r}) = \frac{1}{8\pi}\nabla V_{H}^{\textit{tot}}(\vec{r}) \cdot \nabla V_{H}^{\textit{tot}}(\vec{r}), \end{equation} | (9) |
| \begin{equation} \varepsilon_{H}(\vec{r}) = - \frac{1}{8\pi}V_{H}^{\textit{tot}}(\vec{r})\nabla^{2}V_{H}^{\textit{tot}}(\vec{r}). \end{equation} | (10) |
The stress density is expressed by a sum of several density terms as
| \begin{align} \tau_{\alpha \beta}(\vec{r}) & = \tau_{\alpha \beta}^{\textit{kin}}(\vec{r}) + \tau_{\alpha \beta}^{\textit{xc}}(\vec{r}) + \tau_{\alpha \beta}^{H}(\vec{r}) + \tau_{\alpha \beta}^{\textit{local}}(\vec{r}) \\ &\quad + \tau_{\alpha \beta}^{\textit{non-local}}(\vec{r}) + \tau_{\alpha \beta}^{\textit{ovrl-self}}(\vec{r})+ \tau_{\alpha \beta}^{\textit{overlap}}(\vec{r}). \end{align} | (11) |
| \begin{align} E_{\textit{tot}}^{\varepsilon} & = \int_{\det | I + \varepsilon|\Omega}\varepsilon_{\textit{tot}}((I + \varepsilon)^{-1}\vec{r})d\vec{r} = \int_{\Omega}\eta^{\varepsilon}(\overrightarrow{r'})\det | I + \varepsilon|d\overrightarrow{r'}\\ & \approx (1 + Tr(\varepsilon))\int_{\Omega}\eta^{\varepsilon}(\vec{r})d\vec{r}. \end{align} | (12) |
| \begin{equation} \begin{split} &\psi_{i}(\vec{r}) \to \psi_{i}^{\varepsilon}(\vec{r}) = \det | I + \varepsilon|^{-1/2}\psi_{i}((I + \varepsilon)^{-1}\vec{r}) \\ &\rho_{e}(\vec{r}) \to \rho_{e}^{\varepsilon}(\vec{r}) = \det | I + \varepsilon|^{-1}\rho_{e}((I + \varepsilon)^{-1}\vec{r}) \\ & V_{a}(\vec{r} - \overrightarrow{R_{a}})\to V_{a}^{\varepsilon}(\vec{r} - \overrightarrow{R_{a}}) = V_{a}(\vec{r} - (I + \varepsilon)\overrightarrow{R_{a}}) \end{split} \end{equation} | (13) |
| \begin{equation*} \Omega \to \Omega^{\varepsilon} = \det |I + \varepsilon |\Omega. \end{equation*} |
By applying the above procedure to each energy-density term in eq. (7), we obtained the corresponding stress-density term in eq. (11). Some terms or parts were derived more easily via reciprocal-space expression instead of the above real-space procedure. The compensation charges in $\rho _{e}(\vec{r})$ in the PAW method10) were treated differently from the soft-charge part treated in eq. (13). The term $\tau _{\alpha \beta }^{\textit{overlap}}(\vec{r})$ was additionally derived from the overlap operator.
Details of the stress-density terms are given in Ref. 17). The kinetic-stress density term is expressed as follows:
| \begin{align} \tau_{\alpha \beta}^{\textit{kin}}(\vec{r}) & = \sum_{i}^{\textit{occ}}f_{i}\psi_{i}^{*}(\vec{r})\nabla_{\alpha}\nabla_{\beta}\psi_{i}(\vec{r})\\ & = - \sum\nolimits_{i}^{\textit{occ}}f_{i}\nabla_{\alpha}\psi_{i}^{*}(\vec{r})\nabla_{\beta}\psi_{i}(\vec{r}) + \frac{1}{2}\nabla_{\alpha}\nabla_{\beta}\rho_{e}(\vec{r}). \end{align} | (14) |
The integration of the final forms of eqs. (7) and (11) throughout the supercell results in the same values of Etot and σαβ by the plane-wave PAW method without using the energy or stress density, which proves the validity of eqs. (7) and (11). Each two forms of the kinetic and electrostatic terms in the energy and stress densities are integrated to be the same in the full supercell, because each gauge-dependent term is integrated to be zero in the full supercell.
2.3 Partitioning into gauge-independent local regions 2.3.1 Conditions to remove the gauge dependencyBy using self-consistent charge-density distribution and occupied wave functions, we obtain the energy and stress densities as data on uniform FFT (fast Fourier transformation) mesh grids in a supercell, while the density data themselves contain ambiguities because of the gauge dependency. In the early applications, layer integration or cyclic macroscopic averaging of planar-averaged densities was performed for a surface-slab supercell,20,21) or Wigner-Seitz cell integration centered at each position was used to analyze a point defect53) without rigorous examination of the gauge-dependent terms. Otherwise, the supercell of a surface (defective) system was divided into bulk and surface (defect) parts by symmetrical planes or Voronoi polyhedra so that the integration of the density data in each part could provide a unique physical quantity.20,51,54)
In our scheme, the gauge dependency is directly removed by partitioning a supercell into each local region where the gauge-dependent term is integrated to be zero. The conditions that the gauge-dependent terms of the kinetic-energy and kinetic-stress densities in eqs. (8) and (14) are integrated to be zero in local regions Ωi and Ωj, respectively, are written as follows:
| \begin{equation} \int_{\Omega_{i}}\nabla^{2}\rho_{e}(\vec{r})d\vec{r} = \int_{S_{i}}\nabla \rho_{e}(\vec{r}) \cdot \vec{n}dS = 0 \end{equation} | (15) |
| \begin{align} \int_{\Omega_{j}}\nabla_{\alpha}\nabla_{\beta}\rho_{e}(\vec{r})d\vec{r} &= \int_{\Omega_{j}}[\nabla_{\beta}\rho_{e}(\vec{r})]_{x_{\alpha} = x_{\alpha}^{1}}^{x_{\alpha} = x_{\alpha}^{2}}dx_{\beta}dx_{\gamma} \\ &= \int_{S_{j}}\nabla_{\beta}\rho_{e}(\vec{r})\overrightarrow{e_{\alpha}} \cdot \vec{n}dS = 0. \end{align} | (16) |
Practical partitioning to satisfy the above conditions was first performed in a fcc-Al(111) surface-slab supercell, consisting of several (111) atomic layers and a vacuum region.17) Along the surface-parallel directions, denoted as x and y, the border of the two-dimensional period in the slab can be used as the partitioning border to satisfy eqs. (15) and (16) because of the presence of opposite directions of $\vec{n}$. About the direction normal to the surface, denoted as z, we set a flat x-y plane between two (111) atomic layers at a proper height z. We consider the condition of the surface integration of the gradient of $\rho _{e}(\vec{r})$ in the z direction, throughout the x-y plane area S in the two-dimensional period as follows:
| \begin{equation} \int_{S}\frac{\partial \rho_{e}(\vec{r})}{\partial z}dxdy = \frac{\partial}{\partial z}\int_{S}\rho_{e}(\vec{r})dxdy = \frac{\partial \overline{\rho_{e}}(z)}{\partial z} = 0, \end{equation} | (17) |
The Bader partitioning55) was first adopted by Trinkle and co-workers37) in similar local-energy calculations. This partitions a supercell by a region border on which the surface-normal component of the gradient of $\rho _{e}(\vec{r})$, namely $\nabla \rho _{e}(\vec{r}) \cdot \vec{n}$ in eq. (15), is zero, usually leading to each atomic region. Yu and Trinkle56) developed an efficient Bader-partitioning algorithm for the data of $\rho _{e}(\vec{r})$ on uniform grids by estimating a weight of each grid point near the border via equations of probability flux. The Bader partitioning is also applicable to a diagonal sum of a local-stress tensor, namely a hydrostatic pressure, because a diagonal sum of eq. (16) is equivalent to eq. (15). Thus, the Bader region is used for a gauge-independent local hydrostatic stress of each atom in addition to a gauge-independent atomic energy.18) We can also obtain local energies and stresses of various atomic groups by summation of atomic energies and by volume averaging of atomic stresses as eqs. (5) and (6), respectively, which is effective to remove ambiguities associated with the charge transfer between different species as discussed later in Sec. 4.5.
Trinkle and co-workers37) dealt with the gauge-dependent problem in the electrostatic-energy density, additionally. The gauge-dependent term in the electrostatic-energy density has a form of $\frac{1}{8\pi }\nabla \cdot (V_{H}^{\textit{tot}}(\vec{r})\nabla V_{H}^{\textit{tot}}(\vec{r}))$ as mentioned in Sec. 2.2.1. The condition that the integration of this term in a local region becomes zero is expressed as $\int_{S}V_{H}^{\textit{tot}}(\vec{r})\nabla V_{H}^{\textit{tot}}(\vec{r}) \cdot \vec{n}dS = 0$ for the region border, similarly to eq. (15). This can be satisfied if the normal gradient of $V_{H}^{\textit{tot}}(\vec{r})$ on the region border is zero. Such a region can be found by a method similar to the Bader partitioning via replacing $\rho _{e}(\vec{r})$ with $V_{H}^{\textit{tot}}(\vec{r})$ in eq. (15). They dealt with two kinds of local (atomic) regions to settle the two kinds of gauge-dependent problems in the energy density.
In our scheme, however, we adopt Maxwell forms for the electrostatic energy and stress densities, which are integrated in the same local regions decided for the kinetic energy and stress densities. We do not directly deal with the gauge-dependent problems in the electrostatic energy and stress, because the Maxwell forms can be regarded as a more proper way to define the distribution of electrostatic energy and stress as discussed in the original works.15,20,21) The gauge-dependent problem in the kinetic energy and stress densities is more essential, associated with the inherent quantum nature of electrons, which enforces us to define kinetic energy and stress of electrons within some definite spaces, in contrast to the gauge-dependent problem in the electrostatic energy and stress as a matter of choice for effective expressions as discussed later in Sec. 3.2.
An alternative version of the energy density was proposed by Cohen and co-workers.57) The DFT total energy is expressed as a sum of the electron kinetic energy, the electron-electron, electron-nuclei (ion), and inter-nuclei (ionic) electrostatic energies, and the exchange-correlation energy, while alternatively the total energy is given as a sum of following terms: the sum of occupied Kohn-Sham (KS) eigen energies, the double-counting corrections for the electron-electron and exchange-correlation energies, and the inter-nuclei (ionic) electrostatic energy, as discussed later in Sec. 3.4. Cohen et al. dealt with the density form of this alternative expression. The first term as the KS energy density $\varepsilon _{\textit{KS}}(\vec{r})$ is expressed as
| \begin{equation} \varepsilon_{\textit{KS}}(\vec{r}) = \int_{- \infty}^{E_{F}}ED(\vec{r},E)dE \end{equation} | (18) |
| \begin{equation} D(\vec{r},E) = \sum\nolimits_{i}f_{i}| \psi_{i}(\vec{r})|^{2}\delta (E - E_{i}). \end{equation} | (19) |
| \begin{equation} E_{\textit{band}}^{I} = \int_{\Omega_{I}}\varepsilon_{\textit{KS}}(\vec{r})d\vec{r} = \int_{- \infty}^{E_{F}}ED_{I}(E)dE. \end{equation} | (20) |
| \begin{equation} D_{I}(E) = \int_{\Omega_{I}}D(\vec{r},E)d\vec{r} = \sum\nolimits_{i}f_{i}\int_{\Omega_{I}}|\psi_{i}(\vec{r})|^{2}d\vec{r}\delta (E - E_{i}). \end{equation} | (21) |
The concepts of microscopic and macroscopic stresses of a quantum-mechanic system were introduced in the early days of quantum mechanics.40,41) The understanding on this issue was greatly advanced by Nielsen and Martin,15,16) followed by several theoretical studies.20,21,58–61) Nielsen and Martin clarified the relation among the force, stress and virial theorems in a quantum-mechanical system. By the stress theorem, they derived explicit formulation of an ab initio macroscopic stress tensor in the plane-wave DFT method via introducing scaling into the ground-state expression of the total energy as mentioned in Sec. 2.2.2. They also derived a general form of a microscopic stress-tensor field with atomic-scale variation, so that its divergence provides a force density field, resulting in kinetic (momentum-flux) and potential (configuration-dependent) terms, while it has inevitable ambiguity for the addition of the curl of an arbitrary tensor field, similarly to a classical stress field. About the potential term, the Kugler form62) was discussed as a general one, while the Maxwell form was preferred as a convenient one, due to relatively short-ranged or localized expression.15)
In successive studies,58–61) a microscopic stress-tensor field in a quantum-mechanical system was also derived in different manners, respectively, so as to generate a correct force density field by its divergence or to preserve the momentum balance, leading to essentially similar features of the stress-tensor field, consisting of kinetic and potential terms, while the inherent non-uniqueness of the stress field permits some differences. The Maxwell form was also preferred for the potential term in each study, compared to a general Kugler form.
The stress density by Filippetti and Fiorentini,21) adopted by our scheme, is a practical microscopic stress field in the framework of the plane-wave DFT method, following the discussion by Chetty and Martin20) on both the energy and stress densities. The symmetric and asymmetric forms for the kinetic-stress density were discussed, and the Maxwell form was adopted as the electrostatic-stress density. These terms as well as the exchange-correlation term have similarity to those in the other stress fields.15,58–61)
3.3 Local energy and local stress by LCAO representationFor a quantum-chemical method using LCAO representation and DFT, the energy-density analysis (EDA) was proposed as a tool to partition the total energy into the contributions of respective atoms or bonds.63,64) The exchange-correlation energy is decomposed into each atom spatially by Voronoi partitioning,65) and the other terms are dealt with by Mulliken-type atomic-orbital decomposition of the LCAO representation.
Feilbelman developed the method of ab initio stress-tensor calculations by the first-principles LCAO method using NCPP and LDA.66) A stress tensor was obtained by explicit LCAO expression of the total-energy derivative as a counterpart of the Nielsen-Martin scheme15,16) for the plane-wave DFT method, where Pulay corrections for the atomic-orbital basis were rigorously introduced. Just like local energies by the atomic-orbital decomposition of the total energy, local stresses can be obtained by the atomic-orbital decomposition of the LCAO representation of the stress tensor in addition to spatial partitioning of the exchange-correlation part.66,67)
For the Al(111) surface, layer-resolved stresses were obtained by this LCAO scheme.67) Figure 1 shows the comparison with our results by stress-density integration in each layer region explained in Sec. 2.3.2.17) The agreement is satisfactory in spite of the differences between LDA and GGA and between the basis functions. There exist large in-plane tension and compression at the top and second layers in both the results. For the third and fourth layers, our results show small in-plane tension, and the LCAO results show small in-plane compression, while the total oscillation features of in-plane layer stresses are rather similar. The oscillation of the in-plane stresses is correlated with the Friedel-type oscillation of the in-plane bond charges of surface layers, caused by remarkable redistribution of sp electrons of Al17) as will be discussed in Sec. 4.1.
Ab initio local stress of an Al (111) surface modeled by a nine-layer slab.17) Surface-parallel components of layer-resolved stresses, calculated by our scheme for both relaxed and unrelaxed surface configurations, are plotted, compared with Feibelman’s results using the first-principles LCAO method.67) Positive and negative values mean tensile and compressive stresses in units of energy without division by a local volume in eq. (4), because the local volume is not defined for the surface top layer. For thicker slabs, similar stress oscillations at the surface layers were observed.17)
Historically, local-energy changes by defects or surfaces in solids were evaluated from local electronic-structure changes through LDOS calculations based on the Green’s function formalism, implemented by the multiple-scattering or TB methods.42–45,50,68–70) We compare our formulations of local energy and local stress with those by the Green’s function formalism. Here we consider a supercell system with local ionic potentials and LDA for simplicity. The total energy is expressed as follows:
| \begin{align} E_{\textit{tot}} & = \sum_{i}^{\textit{occ}}f_{i}\int \psi_{i}^{*}(\vec{r})\left(- \frac{1}{2}\nabla^{2}\right)\psi_{i}(\vec{r})d\vec{r} + \int V_{\textit{ion}}(\vec{r})\rho_{e}(\vec{r})d\vec{r} \\ &\quad + \frac{1}{2}\int V_{H}(\vec{r})\rho_{e}(\vec{r})d\vec{r}+ \int \varepsilon_{\textit{xc}}(\vec{r})\rho_{e}(\vec{r})d\vec{r} + E_{I-I}, \end{align} | (22) |
| \begin{align} E_{\textit{tot}} & = \sum_{i}^{\textit{occ}}f_{i}E_{i} - \frac{1}{2}\int V_{H}(\vec{r})\rho_{e}(\vec{r})d\vec{r} \\ &\quad + \int (\varepsilon_{\textit{xc}}(\vec{r}) - \mu_{\textit{xc}}(\vec{r}))\rho_{e}(\vec{r})d\vec{r} + E_{I-I}\\ & = \int_{- \infty}^{E_{F}}ED(E)dE - \frac{1}{2}\int V_{H}(\vec{r})\rho_{e}(\vec{r})d\vec{r} \\ &\quad + \int (\varepsilon_{\textit{xc}}(\vec{r}) - \mu_{\textit{xc}}(\vec{r}))\rho_{e}(\vec{r})d\vec{r} + E_{I-I}, \end{align} | (23) |
We consider the decomposition of the total energy into each atomic-region contribution EI as $E_{\textit{tot}} = \sum\nolimits_{I}E_{I} $ for eqs. (22) and (23) as follows:
| \begin{align} E_{I} & = \sum_{i}^{\textit{occ}}f_{i}\int_{\Omega_{I}}\psi_{i}^{*}(\vec{r})\left(- \frac{1}{2}\nabla^{2}\right)\psi_{i}(\vec{r})d\vec{r} \\ &\quad + \int_{\Omega_{I}}V_{\textit{ion}}(\vec{r})\rho_{e}(\vec{r})d\vec{r} + \frac{1}{2}\int_{\Omega_{I}}V_{H}(\vec{r})\rho_{e}(\vec{r})d\vec{r}\\ &\quad + \int_{\Omega_{I}}\varepsilon_{\textit{xc}}(\vec{r})\rho_{e}(\vec{r})d\vec{r} + e_{I} \end{align} | (24) |
| \begin{align} E_{I} & = \int_{- \infty}^{E_{F}}ED_{I}(E)dE - \frac{1}{2}\int_{\Omega_{I}}V_{H}(\vec{r})\rho_{e}(\vec{r})d\vec{r}\\ &\quad + \int_{\Omega_{I}}(\varepsilon_{\textit{xc}}(\vec{r}) - \mu_{\textit{xc}}(\vec{r}))\rho_{e}(\vec{r})d\vec{r} + e_{I}, \end{align} | (25) |
Equation (25) is used in the Green’s function formalism, implemented by the multiple-scattering or TB methods.42–45,50,68–70) The relation between eqs. (24) and (25) indicates that our local energy possibly has the physical meaning, similar to that by the Green’s function formalism. For eq. (25), Crampin and co-workers45) performed Green’s function calculations of layer-resolved local energies of stacking faults (SFs) or twins in fcc metals as the difference from the bulk total energy within the perturbation theory. For the first term in eq. (25), KS energies and LDOSs were calculated by the accurate multiple-scattering method (layered KKR method) with the bulk potential in the atomic sphere (frozen potential approximation), and then the bulk correspondence was subtracted. As for the remaining terms in eq. (25), the difference from the bulk crystal was ignored, because of no substantial changes of atomic volumes or short-range order in SFs or twins, as validated in closed-packed metals by the force theorem.71,72) The results show substantial energy increases at layers in local hcp-like stacking with peculiar layer-energy profiles depending on kinds of SFs or twins and metal species. We also obtained similar layer-energy profiles for SFs or twins in fcc metals by our local-energy scheme, and the agreement is satisfactory.
Such perturbative treatment of the remaining terms of eq. (25) in Ref. 45) cannot be validated for general defects, and each term should be calculated. Recently, Egami’s group50) has done it by using a similar multiple-scattering Green’s function method,73) and obtained atomic contributions to both the band-structure energy (the first term) and the remaining terms in eq. (25), leading to local atomic energies. As for atomic stresses, numerical derivation of each atomic energy was performed by practical finite deformation of a system to include the effect of charge redistribution, of which the importance was pointed out in their formulation.50) Respective terms in the atomic energy and stress are obtained for each Voronoi region.
This scheme was applied to atomic hydrostatic stresses (pressures) in fcc complex concentrated alloys or high-entropy alloys (HEAs) consisting of multiple transition-metal (TM) elements,74–76) where apparent correlation between the atomic stress and charge transfer, defined by the Voronoi partitioning, was observed. Atoms gaining (losing) electrons show compressive (tensile) stresses. The authors discussed the relation between the atomic-stress fluctuation and the yield property of the alloys.76) On the other hand, Ishibashi and co-workers77) have applied our local-stress scheme to bcc HEAs consisting of refractory TM elements, and observed similar correlation between the atomic stress and charge transfer, while they observed that the results by the Bader partitioning can represent the effects of local chemical environments more reasonably, than those by the Voronoi partitioning, by incorporating the effects of realistic charge distribution between different species. The atomic charge, atomic volume and atomic stress of each constituent atom, defined by the Bader partitioning, show clear mutual correlation, which can be explained by the features of the valence electron concentration (VEC) of each atom and its first and second neighbor shells.77) The issue of the partitioning methods is discussed again in Sec. 5.
3.5 Local energy by TB or Wannier-function representationTB or related methods were also used in eq. (25) for surfaces or defects in TMs or semiconductors, already in the early days.42–44,68–70) The first term in eq. (25) is derived by the on-site TB representation70) of Eband (the first term of eq. (23)), using the atomic-orbital basis set {|χIL⟩} where I denotes an atomic site and L = lm as follows:
| \begin{align} E_{\textit{band}} & = \sum\nolimits_{i}^{\textit{occ}}f_{i}\langle \psi_{i}| H|\psi_{i} \rangle \\ &= \sum\nolimits_{i}^{\textit{occ}}f_{i}\langle \psi_{i} |\sum\nolimits_{IL}| \chi_{IL} \rangle \langle \chi_{IL}|H| \psi_{i} \rangle\\ & = \sum\nolimits_{i}^{\textit{occ}}f_{i}E_{i}\sum\nolimits_{IL}| \langle \chi_{IL}|\psi_{i} \rangle |^{2} \\ &= \sum\nolimits_{IL}\int_{- \infty}^{E_{F}}E\sum\nolimits_{i}f_{i}| \langle \chi_{IL}|\psi_{i} \rangle |^{2}\delta (E - E_{i})dE\\ & = \sum\nolimits_{IL}\int_{- \infty}^{E_{F}}ED_{IL}(E)dE = \sum\nolimits_{I}\int_{- \infty}^{E_{F}}ED_{I}(E)dE. \end{align} | (26) |
Eband can be described also by the inter-site TB representation70) via further inserting $\sum\nolimits_{I'L'}| \chi _{I'L'} \rangle \langle \chi _{I'L'} | = 1$ into eq. (26) as follows:
| \begin{align} E_{\textit{band}} & = \sum_{i}^{\textit{occ}}f_{i}\langle \psi_{i}| H|\psi_{i} \rangle\\ & = \sum\nolimits_{i}^{\textit{occ}}f_{i}\langle \psi_{i} |\sum\nolimits_{IL}| \chi_{IL} \rangle \\ &\quad \times \langle \chi_{IL}|H\sum\nolimits_{I'L'}| \chi_{I'L'} \rangle \langle \chi_{I'L'}| \psi_{i} \rangle\\ & = \sum_{ILI'L'}H_{ILI'L'}\int_{- \infty}^{E_{F}}\sum_{i}f_{i}\langle \chi_{IL}| \psi_{i} \rangle^{*}\\ &\quad \times \langle \chi_{I'L'}|\psi_{i} \rangle \delta (E - E_{i})dE\\ & = \sum\nolimits_{ILI'L'}H_{ILI'L'}\int_{- \infty}^{E_{F}}D_{ILI'L'}(E)dE. \end{align} | (27) |
For semiconductors or insulators, another local-energy expression is possible from the viewpoint of Wannier functions.80,81) The band-structure energy in eq. (23) can be changed as follows:
| \begin{align} E_{\textit{band}} & = \sum_{i}^{\textit{occ}}E_{i} = \sum_{i}^{\textit{occ}}\langle \psi_{i}| H|\psi_{i} \rangle\\ & = \sum\nolimits_{l}\langle \phi_{l}| H|\phi_{l} \rangle = \sum\nolimits_{l}E_{l}, \end{align} | (28) |
As mentioned in Sec. 1, Egami and Vitek47–49) analyzed atomic-level stresses in amorphous metals by using classical interatomic potentials such as EAM. For a system described by EAM, the energy of an atom i is given by $E_{i} = \frac{1}{2}\sum\nolimits_{j}V(r_{ij}) + F(\overline{\rho _{i}})$ as $E_{\textit{tot}} = \sum\nolimits_{i}E_{i} $, where j denotes an interacting neighbor of the atom i. V(r) is a pair potential and $F(\overline{\rho _{i}})$ is the embedding energy as a function of the density $\overline{\rho _{i}}$ at the atom i given by the atomic density function ρ(rij) as $\overline{\rho _{i}} = \sum\nolimits_{j}\rho (r_{ij}) $. The stress tensor is obtained by the strain derivative of Etot as $\sigma _{\alpha \beta } = \frac{1}{\varOmega }\frac{\partial E_{\textit{tot}}}{\partial \varepsilon _{\alpha \beta }}$ in common with the Nielsen-Martin scheme.15,16,49) By assuming eq. (6), we obtain the local stress as the contribution of the atom i,
| \begin{align} \sigma_{\alpha \beta}^{i} &= \frac{1}{\varOmega_{i}}\frac{\partial E_{i}}{\partial \varepsilon_{\alpha \beta}} = \frac{1}{2\varOmega_{i}}\sum\nolimits_{j}V'(r_{ij})\frac{r_{ij,\alpha}r_{ij,\beta}}{r_{ij}} \\ &\quad + \frac{1}{\varOmega_{i}}F'(\overline{\rho_{i}})\sum\nolimits_{j}\rho '(r_{ij})\frac{r_{ij,\alpha}r_{ij,\beta}}{r_{ij}}. \end{align} | (29) |
In Ref. 18), we made a comparison on the atomic stress and energy between our ab initio and EAM86,87) results for coincidence-site lattice (CSL) GBs in Cu and Al as shown in Figs. 2 and 3. The agreement on the atomic stress is rather good for the Cu GB, while the Al GB shows substantial differences in the GB-core region. The agreement on the atomic energy is better also in the Cu GB than in the Al GB, while both the GBs show substantial differences in the GB-core regions. The relatively poor agreement for the Al GB on both the atomic stress and energy is caused by remarkable charge redistribution at GB-core atoms, due to the nature of Al sp electrons, of which the effects cannot be properly represented by EAM. However, the sum of atomic energies at the GB-core region is properly reproduced by EAM. In both the Cu and Al GBs, there is a tendency that the agreement is better for atomic stresses than atomic energies, which is related to the fact that the local stress is an intensive quantity divided by a local volume, compared to the local energy as an extensive quantity. The latter substantially depends on the volume or border position of the local region, while the former does not. This makes the local energy more sensitive to the charge redistribution at defects via region-border changes, compared to the local stress. This should be the reason why the agreement on the atomic local energies at GB-core regions between the ab initio and EAM methods is worse.
Ab initio local stresses in the {221}Σ9 tilt GBs in Cu (a) and Al (b).18) Both the Cu and Al GBs have similar atomic configurations. Diagonal sums of local-stress tensors for respective atoms in the 42-atom supercell are plotted in units of GPa, compared with those by EAM potentials.86,87) Atoms from No. 8 (29) to No. 15 (36) belong to the GB-core region. Atoms of Nos. 10 (31) and 13 (34) are looser sites, and those of Nos. 11 (32) and 12 (33) are tighter sites. The numbers in the parentheses indicate atoms in the other symmetric interface.
The surface stress of metals dominates various phenomena such as surface reconstruction, surface growth, formation of mesoscopic structures, and so on. There are two models on the origin of surface stresses of metals;46) charge redistribution and orbital rehybridization. These two factors are not strictly separated, and depend on the nature of bonding of each metal. In our ab initio results of layer-resolved surface stresses of the Al (111) surface17) explained in Sec. 3.3 (Fig. 1), we found remarkable oscillation of in-plane stresses of the surface layers, which has direct correlation with Friedel-type oscillation of in-plane bond charges of the surface layers. The increased (decreased) in-plane bond charge induces the tensile (compressed) in-plane stress. The remarkable charge redistribution, typical of Al sp electrons, apparently causes the surface stresses of Al. We did not observe such surface-stress oscillation caused by surface-charge oscillation in fcc (111) surfaces of 4d or 5d TMs.22,23)
On the other hand, for the surface stresses of the series of 4d and 5d TMs,22,23) our ab initio local stresses agree with Friedel stresses obtained by the second-moment TB approximation44) for the surface-atom energies via eqs. (25) and (26). This indicates that the surface d-band width change due to d-orbital rehybridization, properly described by the second-moment TB approach, is the dominant factor of the surface stresses of 4d and 5d TMs, instead of charge redistribution. Our local-stress scheme has provided valuable insights into the mechanisms of surface stresses of various metals.
4.2 Local energies and stresses at GBs in metals and SiCSL tilt GBs in Cu and Al (Figs. 2 and 3)18) and those in bcc Fe24,25) were examined by our scheme. Substantial energy changes and stresses are observed for atoms in the GB-core regions, consisting of peculiar atomic rings or clusters, and become much smaller in central bulk regions in each supercell. In all the tilt GBs, except for simple twins, there exist both tighter and looser sites. The tighter (looser) sites show compressed (tensile) stresses and relatively lower (higher) energies with relatively smaller (larger) atomic volumes or shorter (longer) bond lengths. This is because atomically-stepped (not smooth) surfaces constitute usual tilt GBs. In bcc Fe, the tighter (looser) sites show reduced (enhanced) spin polarization as explained by the Stoner model. Such contrasting sites greatly affect segregation behaviors of solutes or impurities.26,27,30,31)
In Fig. 3(b), there is a negative atomic energy in the Al GB, compared to the bulk energy, at the tighter site with short bond lengths and charge accumulation.18) This represents local special stability via strong covalent-like bonding for a less-coordinated atom at the GB. As is well known, Al atoms with less coordination numbers at Al surfaces, defects or GBs tend to form such covalent-like strong bonds as typical behaviors of Al sp-valence electrons, also observed in other ab initio studies88–90) and related to the peculiar surface stresses of Al explained in Secs. 3.3 and 4.1. The occurrence of such negative atomic energies in Al tilt GBs has been predicted from the local environment via the machine-learning technique35) as explained in Sec. 4.6.
In the algorithm that we adopt, the Bader partitioning for each atom sometimes fails in fourfold semiconductors such as Si, because of charge density humps away from atomic sites. Thus, for Si GBs, we adopted the Voronoi partitioning65) to obtain atomic stresses. As discussed in Sec. 3.6, the local stress is an intensive quantity divided by a local volume, not seriously depending on the details of a region volume or a region border, and thus the results should be reliable. Figure 4 shows calculated atomic stresses in two ⟨110⟩ tilt GBs in Si.28,29) Only in the GB with a larger rotation angle ({114}Σ9), reconstructed bonds along the ⟨110⟩ axis are inevitable,83) inducing large tensile stresses in the GB-core region, in contrast to the {221}Σ9 GB without ⟨110⟩ reconstructed bonds. This explains high density of segregated oxygen atoms at the {114}Σ9 GB in Si, in contrast to no remarkable enhancement at the {221}Σ9 GB, experimentally detected by atom-probe tomography combined with transmission electron microscopy.28,29)
Ab initio local stresses in the (a) {114}Σ9 (θ = 141.06 deg.) and (b) {221}Σ9 (θ = 38.94 deg.) tilt GBs in Si.28,29) Stress of each bond, estimated by the average of atomic hydrostatic stresses of both-end atoms, is plotted on each atomic configuration obtained by ab initio calculations, viewed along the [110] axis. Only substantial tensile stresses are shown in units of GPa. The atomic configurations are in good agreement with electron microscopy observations.28,29)
There often occurs GB segregation of solute atoms in metals, leading to the reduction of doping efficiency inside grains or serious changes of GB properties such as GB embrittlement. It is desirable to understand the GB segregation energy Eseg and its origin for various solute species in each metal.91) For a substitutional solute, Eseg can be obtained as the total-energy change between the prior system of a clean GB and a solute-containing bulk crystal and the segregation system where a host atom at a GB site was exchanged with the solute atom in the bulk. We consider the local energy change of each atom between the prior and segregation systems, and group the local-energy changes of all the atoms into the four groups as T1∼T4. This is the decomposition of Eseg into the four terms T1∼T4.26,27) T1 is the local-energy change of a host atom at the GB substitution site by moving it to a bulk site. The value of −T1 means the local energy at the substitution site in a clean GB against the bulk state. T2 is the local-energy change of a solute atom from the substituted state in the bulk host to the GB-substituted state. T3 is the sum of local-energy changes of host atoms around the substitution site in the GB by replacing a host atom at the GB substitution site with a solute atom. And T4 is the sum of local-energy changes of host atoms around the substituted solute in the bulk host by replacing the solute atom with a host atom. The value of −T4 means the sum of local-energy changes of surrounding host atoms by replacing a host atom with a solute atom in the bulk host. Details of these terms were explained with figures and equations in Refs. 26, 27). Analyzing these four terms provides valuable insights into the origin of GB segregation.
Eseg and preferential segregation sites for a series of 3d-TM solutes at CSL tilt GBs in bcc Fe were examined.30) Early (late) TMs are preferentially segregated to looser (tighter) sites, while middle TMs (Cr and Mn) are segregated without fixed site preference. TMs at both ends of the 3d series show larger magnitudes of Eseg, while Mn shows a cusp at the center. This variation is similar to that of ab initio results of 3d-TM solute and screw-dislocation interactions in bcc Fe.92) We applied the local-energy decomposition of Eseg as shown in Fig. 5, which was interpreted in combination with the LDOS analysis.30) For each solute species, there is a tendency that the values of T1 and T2 counterbalance each other. For early-TM solutes, much lowered T3 values for looser sites are dominant in Eseg. 3d electrons of an early-TM solute strongly hybridize with minority-spin electrons of Fe neighbors at looser sites, leading to the stabilization of surrounding Fe hosts as represented by T3. For late-TM solutes such as Ni and Cu, much lowered T4 values are dominant in Eseg. Late-TM solutes have large numbers of d electrons, and increased d electrons show localized features in the minority-spin LDOS in bulk Fe, meaning relatively weak hybridization with Fe neighbors in bulk Fe, namely destabilization of such Fe neighbors in bulk Fe as represented by large positive values of −T4. In other words, the repulsion from bulk Fe is the main origin. For Mn, the cusp of Eseg is caused by the cusp of T4, meaning destabilization of Fe neighbors of Mn in bulk Fe. The major d electrons of a Mn solute in bulk Fe align as nearly localized high-spin states with weak Mn–Fe d-d hybridization, as explained by Hund’s rule, while a Mn solute at a Fe GB shows substantial Mn–Fe d-d hybridization. A Mn solute has contrasting effects on Fe neighbors as destabilization in bulk Fe and moderate stabilization in a Fe GB, resulting in the cusp at Mn in the profile of Eseg. This is a typical example of competitive localized and itinerant behaviors of 3d electrons.
Variations of the four terms, (a) T1, (b) T2, (c) T3 and (d) T4, in the local-energy decomposition of Eseg for the segregation at looser and tighter sites (blue and red triangles) of each 3d-TM solute in the Σ11(332) (solid triangles) and Σ3(111) (empty triangles) GBs in bcc Fe.30) The sum of the four terms is equal to Eseg. The meaning of each term is explained in the text.
As for sp-element solutes, Al, Si, P and S, in Fe GBs,31) Eseg is dominated by the difference between the total local-energy change of neighboring Fe atoms induced by a solute atom in the GB-segregation case and that in the bulk-substitution case, namely T3 + T4. Such local-energy changes of Fe atoms adjacent to an sp-element solute depend on solute-Fe distances and atomic-orbital characters of solute species, as explained by the COHP of solute-Fe bonding, mentioned in Sec. 3.5. P and S with shorter extents of atomic orbitals greatly destabilize Fe neighbors in bulk Fe, due to weak sp-d hybridization, while Fe neighbors at tighter sites of Fe GBs are stabilized or not so much destabilized, due to short solute-Fe distances, leading to larger magnitudes of T3 + T4, namely Eseg, for P and S at Fe GBs.
4.4 First-principles tensile tests of metallic GBsA first-principles tensile test (FPTT) is a powerful tool to investigate intrinsic strength and failure of intergranular or interphase boundaries or defective systems, according to the behaviors of atoms and electrons.12) For successive straining of a supercell, atomic relaxation is iterated according to atomic forces by first-principles calculations, leading to natural deformation or failure, where the evolution of local energies and stresses is quite interesting. We applied our scheme to each tensile step in the FPTTs of the CSL tilt GBs in Al and Cu,32) of which the initial configurations were analyzed in Figs. 2 and 3. The evolution of atomic energies and stresses is quite different between the Al and Cu GBs. For the Al GB, there are strengthened and stretched bonds at the tighter and looser sites, respectively. Such inhomogeneity induces complex variations of atomic energies and stresses. Initially, the back bonds of strengthened bonds at the tighter sites are mainly stretched, generating rapid increases of local stresses and energies at the tighter sites, then the strengthened bonds are suddenly weakened after some critical stretching, and finally all the interfacial bonds are broken. In each stage, local energies and stresses of the interface and back atoms are remarkably varied, associated with remarkable local charge redistribution, typical of Al sp electrons. In contrast, the evolution of local energies and stresses in the Cu GB is rather simple and homogeneous, except for the process of the SF introduction in bulk regions, only observed in Cu, because of the smaller SF energy in Cu.
This study did not deal with full behaviors of cracks or dislocations associated with GBs. For this purpose, we need much larger supercells with enough freedom of nucleation, injection, or migration of cracks or dislocations. As for the combination with fracture mechanics, energy changes of respective layers during the FPTT were obtained for an Al twist GB33) as universal embedded local-energy functions of layers against local strains. Such embedded functions independent of a bulk-region size in a supercell are provided only by our local-energy scheme.
4.5 Local elastic constants of microstructures in alloysWe can obtain a local elastic constant via the first derivative of a local stress or the second derivative of a local energy for a local region to investigate the elastic or mechanical properties of microstructures in materials. The atomic hydrostatic pressure of each Bader region by our scheme can be used to calculate the atomic bulk modulus Bi via its volumetric derivative as
| \begin{equation} B_{i} = \frac{\sigma_{i}(+ \delta) - \sigma_{i}(- \delta)}{(\Omega_{i}(+ \delta) - \Omega_{i}(- \delta))/\Omega_{i}^{0}}. \end{equation} | (30) |
This scheme was used to elucidate the mechanism of Si-concentration dependent bulk-modulus changes of Fe-rich Fe–Si alloys.34) Experimentally, the bulk modulus of a bcc-Fe alloy monotonically decreases for the increase of substituted Si concentration until about 10 at%Si, while it suddenly recovers over this concentration.93) We constructed alloy models, where SiFe8 cubic clusters as a bcc cubic unit with a substituted Si atom at a body-center site are arranged in a large bcc-Fe supercell, and the distribution and mutual contact of SiFe8 clusters change according to the Si concentration. Figure 6 shows the variations of local bulk moduli of the SiFe8-cluster and remaining pure-Fe regions, respectively, as well as the supercell-averaged one, against the Si concentration. The supercell-averaged one monotonically decreases until 11.11 at%Si, and recovers over this concentration, in good agreement with the experiments.93) As for each region, the SiFe8-cluster region becomes much weaker from 6.25 at%Si to 11.11 at%Si, and recovers over this concentration, while the variation of the remaining Fe region is small. Until 6.25 at%Si, SiFe8 clusters exist separately, while for the concentration from 6.25 at%Si to 11.11 at%Si, there occur contacts or chains of SiFe8 cubic clusters with diagonal corner Fe atoms shared, which remarkably decrease local bulk moduli of the SiFe8 clusters, especially at the corner Fe atoms shared by two SiFe8 clusters. For the concentration over 11.11 at%Si, there occur three-dimensional or two-dimensional arrangements of edge-shared (closely contacted) SiFe8 clusters, where the cluster regions show remarkably enhanced bulk moduli, because of stronger sp-d Si–Fe hybridization. The manner of arrangement or connection of SiFe8 clusters, depending on the Si concentration, dominates the elastic moduli of the SiFe8-cluster region and the whole supercell. This can be clarified only by examining local bulk moduli of respective regions using our scheme, instead of examining the averaged electronic structure.93) The Bader partitioning between different species may be involved in the charge transfer across the region border, inducing some complexity in physical interpretation of atomic quantities. This difficulty was avoided by merging atomic regions into a SiFe8-cluster region.
Calculated local bulk moduli of (a) pure-Fe and (b) SiFe8-cluster regions as well as (c) the cell-averaged ones for the Fe–Si alloy supercell models with various Si concentrations are plotted in units of GPa.34) At 12.5 at%Si, partial D03 models, I and II, contain two-dimensional and three-dimensional arrangements of edge-shared SiFe8 clusters, respectively. A number on each vertical bar indicates the volume ratio of each region in the Fe–Si alloy supercell.
Machine-learning techniques attract much attention in recent materials science and engineering.94) The combination with ab initio local-energy or local-stress calculations is promising to open new frontiers in this area. One of the purposes of machine-learning techniques is to derive novel prediction or insights, efficiently from experimental or ab initio data. Our scheme can extend the data to include ab initio local energies or stresses, which greatly improves the ability of machine-learning based prediction. This should be significant for complex systems such as surfaces, GBs, mechanical behaviors, solute behaviors, and so on, where local bonding plays a crucial role.
This was realized as machine-learning based prediction of GB energies of CSL tilt GBs in Al by Tamura and co-workers.35) First, typical CSL GBs with relatively small Σ values were selected for constructing training data of local configurations and atomic energies via ab initio local-energy calculations. Second, machine-leaning techniques were applied to find the correlation between the local environment and atomic energy via a linear regression model, LASSO (least absolute shrinkage and selection operator).95) Third, the optimized LASSO model with descriptors of local geometric features of each atom such as SOAP (smoothed overlap of atomic positions)96) was applied to test data of other CSL GBs with larger Σ values to predict atomic energies and GB energies as their sums. Enough prediction accuracy was attained if GB energies were also incorporated into the LASSO optimization.
Tamura and co-workers have extended this approach to large-scale disordered GBs in bcc Fe.36) The regression model was deduced from the data of the local environment and atomic energy, constructed by iterating a small-supercell calculation of an efficiently-selected local site, divided from a large configuration. If the size of the small supercell is enough, charge and energy densities in the central region of the small supercell should be nearly equal to those in the original large configuration, which ensures the accuracy. Each atomic energy may depend on the partitioning method, while the final GB energy is a sum of atomic energies, reducing this ambiguity. Our atomic energy is suitable to represent the effects of local structural disorder in bcc Fe within a rather short distance range, because of the screening by metallic electrons and the rather short extent of d-d hybridization changes as observed in Refs. 24, 25). By the way, finding the correlation between the local environment and atomic energy has the same meaning as the development of interatomic potentials. As for the machine-learning based interatomic potentials, a modified version97) of the ab initio local-energy scheme37) was also used for a neural-network force field of amorphous Si.98)
As for local stresses, machine-learning techniques were applied to find the correlation between an atomic stress and a local chemical environment in refractory HEAs,77) where VEC features of each atom and its neighboring shells were found to decide the atomic stress as mentioned in Sec. 3.4. Further applications of the combination of local energy or local stress with machine-learning techniques are expected for a wide range of phenomena.
Our ab initio local-energy and local-stress calculations are tools to supplement conventional plane-wave DFT calculations by retrieving the distribution information on the energy and stress in a large supercell, which is usually discarded in conventional calculations. There is no general or universal way to partition a supercell into local regions for the energy-density or stress-density integration, even within the gauge-independent conditions, and the values of local energies and stresses inevitably depend on the partitioning method. However, this does not mean that such analyses are ambiguous. Regardless of the partitioning method, the total sum of partitioned local energies and the volume average of partitioned local stresses are rigorously equal to the total energy and the cell stress of the supercell obtained by the conventional method, respectively. To understand the nature, essence or underlying mechanism of complex structures or phenomena in large-supercell DFT calculations, it is quite helpful to analyze the distributions of such local energies and stresses. The Bader partitioning is a natural and practical way to obtain such local distributions in a unified manner as shown in various applications explained above, including applications by Trinkle’s group.37,99,100) Of course, careful interpretation is required for such atomic quantities, especially for local energies, compared to local stresses, because of the difference between extensive and intensive quantities as mentioned in Sec. 3.6. Note that the quantities of atomic Bader regions can be easily extended to those of clusters or atomic groups by summation or volume averaging as mentioned in Secs. 2.3.3 and 4.5. And it is effective to combine supplementary analyses on the electronic behaviors such as LDOS or COHP as performed in Refs. 30, 31).
As explained in Sec. 3, there have been theoretical studies dealing with similar or related local quantities, especially using LCAO methods, Green’s function formulation implemented by the multiple scattering or TB methods, or EAM potentials. The comparison with calculated results by such related schemes is yet limited, while the agreement and differences are basically reasonable as explained in Secs. 3.3, 3.4 and 3.6 on metallic surfaces, SFs, alloys, and GBs. The comparison on the concept or formulation is effective to clarify the physical meanings of our quantities. Of course, there exist differences in the definitions of local energy, local stress and local regions, and in the treatment of the gauge-dependent problems, which prevent simple evaluation of the superiority among the schemes. Our scheme, however, has several advantages on the reduction of ambiguity or artificial operation in the formulation, on the simplicity, flexibility and accuracy in the computational techniques, and on the compatibility with the conventional plane-wave DFT method, sharing the same electronic structure by the plane-wave DFT method. Our scheme builds a bridge between the methods of local viewpoints such as LCAO or Green’s function methods and those of extended viewpoints.
In our scheme, the partitioning of a supercell into atomic, layer or cluster regions as well as the division into defective and bulk parts is possible within the gauge-independence conditions, according to natural valence-charge distribution. In contrast, precise distributions of local quantities by the LCAO methods depend on the atomic orbital basis, while there is no general way to define optimum atomic-orbital basis functions without any ambiguity. In the Green’s function formalism, implemented by the multiple scattering theory, the construction of atomic spheres in the atomic sphere approximation (ASA)45,50) may involve some artificial operation, especially in systems other than closed-packed structures. In the atomic energy and stress by such Green’s function formalism, the Voronoi partitioning is usually used.50) As mentioned in Sec. 3.4, the Bader partitioning, of which the border keeps zero normal gradient of the charge distribution, is more appropriate to provide physically-meaningful atomic quantities, compared to the Voronoi partitioning simply based on the atomic configuration. In HEAs, the Bader partitioning provides atomic quantities, reasonably reflecting local chemical environments.77) In elemental solids, the Bader partitioning should provide an atomic region nearer to that with the local charge neutrality.68) The superiorities of our scheme enable us to make systematic analysis of various systems containing a wide range of elements in a unified manner.
There remain several unsettled issues in theoretical or computational aspects, some of which are important for the future applications in the energy and environment materials such as catalysts and battery materials. First, the effects of charge transfer on the local energy and stress should be generally investigated, although there exists inherent ambiguity in defining a transfer value as an extensive quantity like a local energy. It is desirable to develop new analysis schemes for local energy and stress of ionic solids, compounds, or alloys with substantial charge transfers, where the sum or average of neighboring ionic pairs or clusters might be effective instead of analyzing each value of cationic or anionic atom. Second, concerning the above issue, new partitioning of bond regions instead of atom-centered regions should be investigated. As mentioned in Sec. 4.2, the usual algorithm of the atom-centered Bader partitioning is not necessarily stable for systems of four-fold semiconductors. Third, the effects of different expressions for the electrostatic terms in the energy and stress densities, such as Maxwell and Coulomb forms and Maxwell and Kugler forms, should be further investigated as mentioned in Secs. 2.2, 2.3 and 3.2. Fourth, local regions to satisfy eq. (14) for each stress-tensor component, instead of hydrostatic stresses, should be investigated. Fifth, the relation between local stresses and atomic forces within the present theoretical framework is not so clear, in contrast to general arguments on quantum-mechanical stress and force fields as mentioned in Sec. 3.2.15,20,21,58–61) Also, comparisons with classical stress fields of continuum or atomistic models84,85,101) are quite interesting.
Finally, we emphasize the future applications of local energies to large-scale systems. In Ref. 36) explained in Sec. 4.6, a total energy of a large configuration of bcc Fe was computed via machine-learning based modelling of the relation between the local environment and atomic energy, derived by the data of small-supercell calculations of atomic energies for selected local configurations, divided from the original large system. This becomes accurate if the size of small supercells is enough to reproduce atomic energies in the original large system at the center of each small supercell, which should be possible in metals, because of the screening or averaging of long-range electrostatic interactions. This approach resembles the divide and conquer (DC) type order-N scheme102,103) if all the atoms are dealt with by such small-supercell calculations of divided regions, because the atomic energy in each Bader region can seamlessly occupy the whole system. It should be also possible to perform reliable multi-scale or hybrid computations of large metallic systems by applying such DC-type computations via small-supercell calculations of local energies for greatly disordered parts and applying the machine-learning based prediction for the other parts, for example.
In conclusion, our scheme of ab initio local-energy and local-stress calculations enables us to utilize the information on the energy and stress distributions inside a large supercell, via the integration of the energy and stress densities in each gauge-independent local region within the plane-wave PAW-GGA framework. There have been several theoretical attempts to calculate similar local quantities, especially using similar densities, LCAO methods, Green’s function formulation implemented by the multiple scattering or TB methods, or EAM potentials. The comparison has made clear the physical meaning of local energy and local stress, and showed the superiority of our scheme as the reduced ambiguity, the flexibility and accuracy, and the compatibility with the conventional plane-wave DFT method. However, we must be careful in making physical interpretation on such local quantities, because of the dependence on the partitioning method, especially for local energies as extensive quantities, while the total sum of local energies is rigorously correct and meaningful. It is effective to extend local regions from atoms to clusters, depending on the situation, or to combine electron-behavior analyses such as LDOS or COHP. In our recent applications to metallic surfaces, GBs with and without solutes, and alloys, the local-energy and local-stress calculations have clarified novel aspects of phenomena, and provided novel insights into the mechanism and effective data for novel machine-learning based modelling, which could not be ever attained in the conventional DFT calculations. We have discussed unsettled issues for future wide applications. The combination with the machine-learning techniques and DC-type computations should lead to quantitative calculations of large-scale metallic systems. Our scheme and recent related schemes of local energy and local stress should create new frontiers in materials science and engineering.
This study was supported by MEXT as a social and scientific priority issue (Creation of new functional devices and high-performance materials to support next-generation industries; CDMSI) to be tackled by post-K computer (hp180220, hp190180). This study was also supported by the Elements Strategy Initiative for Structural Materials (ESISM) through MEXT, and by the Research and Development Initiative for Scientific Innovation of New Generation Batteries 2 (RISING2) project (JPNP16001) from the New Energy and Industrial Technology Development Organization (NEDO). This study was partially supported by JST Industry-Academia Collaborative Programs, ‘Materials Strength from Hamiltonian’, by the Grant-in-Aid for Scientific Research on Innovative Areas, ‘Bulk Nanostructured Metals’, and by the Grant-in-Aid for Scientific Research (Nos. 19H05177 and 20K05696). We thank Dr. Shoji Ishibashi for his contributions to the development of the QMAS code and our scheme. We thank Prof. Tomoyuki Tamura for his idea on the combination with machine-learning techniques for large-scale systems. We thank Dr. Hao Wang, Dr. Somesh Kr. Bhattacharya and Dr. Zhuo Xu for their contributions to the examination and application of our scheme. We thank Prof. Yutaka Ohno for fruitful collaboration on Si GBs. We also thank Dr. Kazuma Ito, Dr. Hideaki Sawada, Prof. Shigenobu Ogata, Prof. Ying Chen, and Prof. Tetsuo Mohri for collaborations on Fe GBs and Fe alloys.
