1Center for Research and Development of Higher Education, The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan2School of Computer Science, Tokyo University of Technology, Hachioji, Tokyo 192-0982, Japan3Development Department, Information and Science Solution Division, HULINKS Inc., Chuo, Tokyo 103-0015, Japan
Received February 6, 2018; Accepted April 12, 2018; Published May 16, 2018
A novel tight-binding method is developed, based on the extended Hückel approximation and charge self-consistency, with referring the band structure and the total energy of the local density approximation of the density functional theory. The parameters are so adjusted by computer that the result reproduces the band structure and the total energy, and the algorithm for determining parameters is established. The set of determined parameters is applicable to a variety of crystalline compounds and change of lattice constants, and, in other words, it is transferable. Examples are demonstrated for Si crystals of several crystalline structures varying lattice constants. Since the set of parameters is transferable, the present tight-binding method may be applicable also to molecular dynamics simulations of large-scale systems and long-time dynamical processes.
This article is published by the Physical Society of Japan under the terms of the Creative Commons Attribution 4.0 License. Any further distribution of this work must maintain attribution to the author(s) and the title of the article, journal citation, and DOI.
The tight-binding (TB) method becomes more and more important as a tool of quantum chemistry and physics because of its high ability of molecular dynamics (MD) simulation such as of long-time dynamics of several nano-seconds and of large-scale systems of several ten thousands atoms or more. Now the TB-MD scheme has been applied to a wide variety of materials and phenomena.
Many successful studies have been accumulated in order to improve accuracy and transferability of the extended Hückel approximation (EHA) with the Slater-type orbitals (STO).1–6) The combination of TB methods and the density functional theory (DFT) was proposed with molecular orbitals, and is called DFTB.7–9) Another approach is that of determining the TB-parameters (Slater–Koster parameters), without wavefunctions, directly from the band structures calculated by the first principles method.10–12) The charge self-consistent (CSC) scheme was also proposed and the environmental effect of the electron charge transfer among ions is accurately included.13) Further promising improvement of STO is the utilization of double-ζ STO for s- and p-orbitals,14) which reproduces the electronic band structures in crystalline solids much accurately, because it expands the freedom of the extent of wavefunctions.
We applied the TB-MD simulation successfully to phenomenon of large-scale systems or long-time dynamics of ions in condensed materials, such as the fracture propagation in Si,15) the formation mechanism and helicity of gold nanowire,16) and an ionic liquid of N-methyl-N-propylpiperidinium bis trifluoromethanesulfonyl imide (PP13-TFSI).17,18) Then the electronic structure and the dynamics of lithium ions in lithium superionic conductors (LISICON) were investigated by the long-time quantum MD simulation such as “nano-second” (ns) processes,19,20) based on the tight-binding scheme of the extended Hückel approximation (EHA). Also we developed a novel mathematical formalism of the Krylov subspace method,21,22) which was applied to calculation of Green's functions in the “extra-large scale” systems of ten million or more atoms (100 nm scale systems).
Through these investigations, we have developed a procedure of automatic determination of parameters of EHA in molecular systems with the help of the hierarchical structure of cohesive energies,19,23) and the program package is named TB-PAM (TB-Parametrization for Molecules). Even though, when we used TB-PAM for several crystalline systems with changing lattice parameters, we face the fact that the original EHA is not highly satisfactory, especially on the description of band shift and detail of the total energy depending upon environmental atomic densities. Combination of the TB parameters and the embedded-atom method (EAM)24) may be also promising.
In the present paper, we propose a novel TB method based on the total energy of the local density approximation (LDA) and establish an automatic procedure of determining a set of parameters, and this procedure overcomes the aforementioned difficulty. Our guiding principles of constructing a novel TB method are (i) reproduction of LDA results of the band structure and the total energy, (ii) no time-consumption of determining wavefunctions, and (iii) keeping the universality of the parameters. In Sect. 2, the total energy of LDA is re-examined. Then, in Sect. 3, a novel tight-binding method is formulated, which we name “the total energy assisted TB method based on LDA of DFT” (TE-TB method). We call the set of parameters of the present TE-TB method the “TE-TB parameters”. Then, the numerical procedure of the machine-determination of the parameters is described in more detail in Sect. 4. Section 5 is devoted to presentation and discussion of the calculated results for crystalline Si. The last section is a discussion and summary of the present work.
2. Total Energy of LDA
Let us discuss the total energy of LDA so as to formulate the TB method on the basis of LDA total energy expression. Total energy of LDA can be given as25) \begin{align} E_{\text{Total}}^{\text{LDA}} &= \sum_{\alpha}^{\text{occ}} \varepsilon_{\alpha} -\frac{1}{2}\int d\boldsymbol{{r}}\int\mathrm{d}\boldsymbol{{r}}^{\prime} \frac{n(\boldsymbol{{r}})n(\boldsymbol{{r}}^{\prime})}{|\boldsymbol{{r}}-\boldsymbol{{r}}^{\prime}|}+E_{\text{xc}}[n]\notag\\ &\quad - \int d\boldsymbol{{r}}\,n(\boldsymbol{{r}})\mu_{\text{xc}}(\boldsymbol{{r}})+\sum_{\{a \neq b\}}\frac{Z_{a}Z_{b}}{R_{ab}}, \end{align} (1) where \(\varepsilon_{\alpha}\) is the eigen energy of a band α (the n-th band at a \(\boldsymbol{{k}}\)-point in the Brillouin zone) and the sum is over the occupied bands over the entire Brillouin zone. The first term of Eq. (1) is called the band energy and we write it as \begin{equation} E_{\text{band}} = \sum_{\alpha}^{\text{occ}} \varepsilon_{\alpha}. \end{equation} (2) The electrostatic energy of electron–electron interaction is doubly counted in \(E_{\text{band}}\), which should be subtracted as in the second term of Eq. (1). The third and fourth terms are the exchange–correlation energy and potential, which may be separated into individual local energies in the local density approximation. The last term is the interaction between ions a and b of a distance \(R_{ab}\), where \(Z_{a}\) is the effective atomic number of an ion a (i.e., the atomic number minus the number of core electrons), and the summation \(\sum_{\{a\neq b\}}\) is over all ion pairs. The total (valence) electron density \(n(\boldsymbol{{r}})\) can be assumed to be a sum of atomic electron density as \begin{equation} n(\boldsymbol{{r}}) = \sum_{a} n_{a}(\boldsymbol{{r}}), \end{equation} (3) where the atomic electron density \(n_{a}(\boldsymbol{{r}})\) may overlap with each other.
The second term of Eq. (1) consists of two parts of the on-site and off-site contributions. We separate them and incorporate one into the on-site energy (the sum of on-site contribution) \(E_{\text{on-site}}\) and the other into the lattice energy (the off-site contribution) \(E_{\text{lattice}}\): \begin{align} E_{\text{on-site}} &= E_{\text{xc}}[n]-\int d\boldsymbol{{r}}\,n(\boldsymbol{{r}})\mu_{\text{xc}}(\boldsymbol{{r}})\notag\\ &\quad - \frac{1}{2} \sum_{a} \int d\boldsymbol{{r}}\int\mathrm{d}\boldsymbol{{r}}^{\prime} \frac{n_{a}(\boldsymbol{{r}})n_{a}(\boldsymbol{{r}}^{\prime})}{|\boldsymbol{{r}}-\boldsymbol{{r}}^{\prime}|}, \end{align} (4) \begin{align} E_{\text{lattice}} &= \sum_{\{a \neq b\}} \left\{\frac{Z_{a}Z_{b}}{R_{ab}} -\int d\boldsymbol{{r}}\int\mathrm{d}\boldsymbol{{r}}^{\prime} \frac{n_{a}(\boldsymbol{{r}})n_{b}(\boldsymbol{{r}}^{\prime})}{|\boldsymbol{{r}}-\boldsymbol{{r}}^{\prime}|} \right\}. \end{align} (5) The lattice energy vanishes when the atomic charge \(n_{a}(\boldsymbol{{r}})\) is not overlapped with each other and each atom is electrically neutral.
Now the LDA total energy can be divided into three parts; the band energy, the lattice energy and the on-site energy such as \begin{equation} E_{\text{Total}}^{\text{LDA}} = E_{\text{band}} + E_{\text{lattice}} + E_{\text{on-site}}. \end{equation} (6)
3. Total Energy Assisted Tight-binding Method Based on LDA
We intend to apply the present framework to the MD simulation of long-time dynamics. Then, the fully self-consistent calculation of wavefunctions (atomic orbitals) becomes heavily time-consuming in repeated computations and we do not use this procedure. Once one fixes atomic orbitals, the Hamiltonian and overlap matrix elements can be calculated. Then coefficients of linear combination of atomic orbitals are determined by solving the eigenvalue equation. Next the band energy and local charge density are calculated, and one can evaluate the lattice energy and the on-site energy. This procedure can be achieved during MD simulation of long-time dynamics.20)
Slater-type orbitals
We use Slater-type orbitals (STOs), simply because the overlap matrix elements can be calculated analytically.26) The STO is written, for the orbital of the principal, angular and magnetic quantum numbers \((n, l, m)\) of an atom a, with the spherical harmonics \(Y_{l}^{m}\) as; \begin{equation} \phi_{nlm}^{\text{STO}(a)}(\mathbf{r}) = R_{nl}^{(a)}(r)Y_{l}^{m}(\theta,\varphi). \end{equation} (7) The radial function of STO is given as, in case of single-ζ, \begin{equation} R_{nl}^{(a)}(r) = \sqrt{\frac{(2\zeta_{l}^{a})^{2n+1}}{(2n)!}}r^{n-1}\exp(-\zeta_{l}^{a} r), \end{equation} (8) or, in case of double-ζ,14) \begin{align} R_{nl}^{(a)}(r) &= C_{(1)} \sqrt{\frac{(2\zeta_{l(1)}^{a})^{2n+1}}{(2n)!}}r^{n-1}\exp(-\zeta_{l(1)}^{a} r)\notag\\ &\quad + C_{(2)} \sqrt{\frac{(2\zeta_{l(2)}^{a})^{2n+1}}{(2n)!}}r^{n-1}\exp(-\zeta_{l(2)}^{a} r). \end{align} (9) The coefficients \(C_{(1)}\) and \(C_{(2)}\) should be determined by the normalization condition \begin{equation} {C_{(1)}}^{2}+{C_{(2)}}^{2}+2C_{(1)}C_{(2)}\left\{\frac{4\zeta_{l(1)}^{a}\zeta_{l(2)}^{a}}{(\zeta_{l(1)}^{a}+\zeta_{l(2)}^{a})^{2}}\right\}^{n+\frac{1}{2}} {}= 1. \end{equation} (10)
Molecular orbitals (MO) or band wavefunctions are expressed as a linear combination of STOs as \begin{equation} \psi_{\alpha} = \sum_{a,nlm}C_{a;nlm}^{\alpha} \phi_{nlm}^{\text{STO}(a)}, \end{equation} (11) in which the coefficients \(C_{a;nlm}^{\alpha}\) are determined by solving a generalized eigenvalue equation \begin{equation} H\psi_{\alpha} = \varepsilon_{\alpha} S \psi_{\alpha}. \end{equation} (12) The parameters, \(\zeta_{l(1)}\) and \(\zeta_{l(2)}\) (or \(\zeta_{l}\)), \(C_{(1)}\) and \(C_{(2)}\), should be determined so as to make the calculated band structure agree with that by the LDA calculation for a certain fixed structure. After fixing them, the values of coefficients \(\{C_{a;NM}^{\alpha}\}\) are calculated again with individual atomic configuration at each time step of MD procedure.
Band energy
Hamiltonian in Extended Hückel approximation +
In the formulation of the extended Hückel approximation (EHA), the diagonal elements of the Hamiltonian should be determined empirically as \begin{equation} H_{al,al}^{0} = \varepsilon_{al}^{\text{atom}} \end{equation} (13) and the off-diagonal elements as \begin{equation} H_{al,bl^{\prime}}^{0} = K\frac{H_{al,al}^{0}+H_{bl^{\prime},bl^{\prime}}^{0}}{2}S_{al,bl^{\prime}};\quad (al) \neq (bl^{\prime}), \end{equation} (14) where \(S_{al,bl^{\prime}}\) is the overlap matrix elements. The parameter K can be a constant as in original EHA or the modified Wolfsberg–Helmholtz constant5) of a function of an inter-atomic distance. One of the advantages of this simple form is the fact that there are no particular problem even if atoms enter adjacent positions during MD simulation.
In the original EHA5,6) or DFTB,9,13) the parameter \(\varepsilon_{al}^{\text{atom}}\) is fixed to a value of an isolated atom or ion. However, in this case, the shift of the energy band structure due to changing the lattice constants could not be reproduced. Then we adopt a form of the atomic energy which depends on the local environment as12) \begin{equation} \varepsilon_{al}^{\text{atom}} = a_{al}+b_{al}\rho_{a}^{2/3}+c_{al}\rho_{a}^{4/3}+d_{al}\rho_{a}^{6/3}. \end{equation} (15) The effective atomic density \(\rho_{a}\) represents effects of local environment around an atom a; \begin{equation} \rho_{a} = \sum_{b(\neq a)}\mathrm{e}^{-\lambda^{\rho} R_{ab}}F_{C}^{\rho}(R_{ab}), \end{equation} (16) where the summation is over all neighboring atoms b, and \(\lambda^{\rho}\) is a parameter which ensures the contribution from the neighboring sites. \(F_{C}^{\rho}(R)\) is a smooth cut-off function defined as \begin{equation} F_{C}^{\rho}(R) = (1+\mathrm{e}^{(R-R_{0}^{\rho})/L^{\rho}})^{-1}. \end{equation} (17) We call the parameter set, \(\zeta_{1}\), \(\zeta_{2}\), \(C_{(1)}/C_{(2)}\), \(\lambda^{\rho}\), \(R_{0}^{\rho}\), \(L^{\rho}\), \(a_{al}\), \(b_{al}\), \(c_{al}\), and \(d_{al}\) “the STO parameters”. Now the tight-binding Hamiltonian is defined and one can evaluate the band structures and band energy.
We use values \(L^{\rho}=0.5\) (atomic unit), \(R_{0}^{\rho}=14.0\) (atomic unit).10–12) Other parameters should be determined so that the band structures with varying the lattice constant and different structures can be reproduced. The detailed procedure will be discussed in Sect. 4. The EHA Hamiltonian with the density dependent atomic energy Eq. (15) may be called “extended Hückel approximation + (plus)”.
Fine correction to Hamiltonian
After trying to fit the band structures \(\varepsilon_{\alpha}\) of several crystalline structures of Si to those by the first principles LDA calculation, one may find a small discrepancy remaining in the band energy, e.g., a discrepancy of a few hundreds or a few tens meV. This small discrepancy may be serious in the total band energy \(E_{\text{band}}\) when one compares it with the LDA results of several different structures.
We observe that the discrepancy of the band energy is almost constant with changing lattice constant in each structure. Now we consider this issue with a strategy of a rigid uniform shift of bands, and introduce a rigid shift \(\Delta(x)\) of the band structure energy \(\varepsilon_{\alpha}\) as a function of the relative difference of packing fraction η \begin{equation} x = (\eta-\eta_{\text{S}_{0}})/\eta_{\text{cp}}. \end{equation} (18) Here \(\eta_{S_{0}}\) denotes the packing fraction of a crystalline structure \(S_{0}\) in an equilibrium condition, e.g., “the principal structure” (e.g., diamond structure in the case of Si), and \(\eta_{\text{cp}}\) is that of a certain fixed structure (e.g., closed packed structure). In more local situation like surfaces, we may need to define a packing fraction in more local way. In later section in an example of Si, we fix the form of Δ as \begin{equation} \Delta (x) = a_{\Delta} [1+\tanh\{b_{\Delta} (x-x_{\Delta})\}], \end{equation} (19) since the amount of the shift of the band energy is observed to change in a rather narrow range of x.
The band energy \(E_{\text{band}}\) can be fitted fairly well, once we introduce a rigid shift \(\Delta (x)\) of energy bands. The corresponding replacement in the eigen state formulation Eq. (12) should be \begin{equation} H_{al, bl^{\prime}} \rightarrow H_{al, bl^{\prime}}+\Delta(x) S_{al, bl^{\prime}}. \end{equation} (20)
Charge self-consistent term
One serious difficulty in the conventional TB-scheme may be the fact that a Hamiltonian is independent of electron charge transfer. The CSC calculation proposed by Elstner et al.13) treats the charge transfer effect correctly by the second order perturbation theory. We adopt this scheme and the Mulliken fluctuation charge is calculated at each iteration step.19)
The charge density does not equal the starting one during iteration steps and the deviation of the electron energy can be evaluated as13) \begin{equation} E_{\text{2nd}} = \frac{1}{2} \sum_{a,b} \Delta q_{a} \Delta q_{b} \gamma_{ab}. \end{equation} (21) Here \(\Delta q_{a}\) is the deviation of the resultant Mulliken charge from the starting one. The parameter \(\gamma_{ab}\) is a function of the inter-atomic distance \(R_{ab}\). The diagonal element \(\gamma_{aa}\) is the chemical hardness or the Hubbard parameter \(U_{aa}\) and may be replaced by the absolute value of the ionization energy. The resultant matrix element of the CSC is written as \begin{equation} H_{al,bl^{\prime}}^{\text{CSC}} = \frac{1}{2}S_{al,bl^{\prime}} \sum_{cl^{\prime \prime}}(\gamma_{al, cl^{\prime \prime}}+\gamma_{bl^{\prime}, cl^{\prime \prime}})\Delta q_{cl^{\prime \prime}}. \end{equation} (22) In many cases, we would use a simplified approximation of the following form: \begin{equation} H_{al,bl^{\prime}}^{\text{CSC}} = \delta_{al,bl^{\prime}}U_{aa}\Delta q_{a}. \end{equation} (23)
Lattice energy
Since one has obtained the atomic electron density \(n_{a}(\boldsymbol{{r}})\) so far in the framework of our tight-binding scheme, one can evaluate the lattice (Coulomb) energy \(E_{\text{lattice}}\) which may be small but not identically zero. The rigorous analytical expression of the second term of Eq. (5) will be discussed in detail in Appendix A.
On-site energy
We calculate the band energy and the lattice energy in the present tight-binding framework. Now, let us determine the form of the on-site energy to reproduce the LDA total energy.
The on-site energy \(E_{\text{on-site}}\) in Eq. (4) is, as is already discussed, purely local. Then we assume that the on-site energy can be a sum of local atomic contributions and expressed with local parameters as \begin{align} E_{\text{on-site}} &= \sum_{a} E_{\text{on-site}}(r_{s}^{a},\varrho_{a}), \end{align} (24) \begin{align} E_{\text{on-site}}(r_{s}^{a},\varrho_{a}) &= f(r_{s}^{a})+g(\varrho_{a})+h(r_{s}^{a},\varrho_{a}). \end{align} (25) Here, we introduce two different parameters for the electron charge distribution around an atom a; i.e., one represents short-range circumstances and the other long-range one.
Possible definition of \(r_{s}^{a}\) is a radius of a sphere containing a single valence electron around an atom a (i.e., a local Wigner–Seitz sphere or a Volonoi polyhedron or a nearest neighbor shell per one valence electron). Now we define \(r_{s}^{a}\) as \begin{equation} \frac{1}{r_{s}^{a}} = (N_{\text{val}}^{a})^{1/3}\left[\sum_{b} \left(\frac{1}{R_{ab}}\right)^{6}\right]^{1/6}, \end{equation} (26) where \(N_{\text{val}}^{a}\) is the number of valence electrons in an atom a. The parameter \(r_{s}^{a}\) in Eq. (26) may be proportional to the nearest neighbor distance and becomes smaller slowly with larger number of neighbor atoms.
The parameter \(\varrho_{a}\) is an effective atom density around an atom a, which may be defined as \begin{equation} \varrho_{a} = \sum_{b(\neq a)}\mathrm{e}^{-\Lambda^{\varrho} R_{ab}}F_{C}^{\varrho}(R_{ab}). \end{equation} (27) The cut-off function \(F_{C}^{\varrho}(R)\) is assumed to be in the same form of \(F_{C}^{\rho}(R)\) in Eq. (17), but with different parameters. Therefore, \(r_{s}^{a}\) represents the short-range environment and \(\varrho_{a}\) the long-range one.
Let us define the functions \(f(r_{s}^{a})\), \(g(\varrho_{a})\), and \(h(r_{s}^{a},\varrho_{a})\) in the following general forms of power expansion; \begin{align} f(r_{s}^{a}) &= k_{m1}^{a} (r_{s}^{a})^{-1}+k_{m2}^{a}(r_{s}^{a})^{-2}+k_{m3}^{a}(r_{s}^{a})^{-3}, \end{align} (28) \begin{align} g(\varrho_{a}) &= A_{a} + B_{a} \varrho_{a}^{2/3} +C_{a}\varrho_{a}^{4/3}+D_{a}\varrho_{a}^{2}+ E_{a}\varrho_{a}^{8/3}\notag\\ &\quad + F \varrho_{a}^{10/3} + G\varrho_{a}^{4}, \end{align} (29) \begin{align} h(r_{s}^{a},\varrho_{a}) &= \{n_{1}^{a}(r_{s}^{a})^{-1} + n_{3}^{a}(r_{s}^{a})^{-2} + n_{5}^{a}(r_{s}^{a})^{-3}\} \varrho_{a}^{2/3}\notag\\ &\quad + \{n_{2}^{a}(r_{s}^{a})^{-1} + n_{4}^{a}(r_{s}^{a})^{-2} + n_{6}^{a}(r_{s}^{a})^{-3}\} \varrho_{a}^{4/3}. \end{align} (30) The values of the parameters \(R_{0}^{\varrho}\), \(L^{\varrho}\), and \(\Lambda^{\varrho}\) may be fixed a priori from physical view points. These parameters (\(R_{0}^{\varrho}\), \(L^{\varrho}\), \(\Lambda^{\varrho}\), \(k_{m1}\sim k_{m3}\), \(A\sim G\), and \(n_{1}\sim n_{6}\)) are named “the local energy parameters”, and should be determined for \(E_{\text{on-site}}\) to reproduce the behavior of LDA results in several bulk crystalline structures as many as possible.
The parameters which we should determine in the TE-TB method are the STO and the local energy parameters, and we call them together “the TE-TB parameters”. The expansion form Eq. (28)–(30) is certainly not unique. The most important point is that the origin of the on-site energy is made clear and that the terms can be represented by a quite general form.
Total energy assisted tight-binding method based on LDA and its transferability
We have formulated a novel tight-binding method whose framework is based on the expression of the LDA total energy Eq. (6) and name this method the total energy assisted tight-binding method based on LDA of DFT (or Total Energy-TB method or simply TE-TB method).
We have fixed first STOs with a few STO parameters in the principal structure S0 with a fixed lattice constant as described in Sect. 3.1. Then remaining STO parameters in Eqs. (15) and (16) are determined to get overall agreement of band structures in the same structure S0 with changing lattice constant. The coefficients \(\{C_{a;nlm}^{\alpha}\}\) in Eq. (11) are determined by solving the eigenvalue equation in individual atomic configuration. The local energy parameters (\(k_{m1}\sim k_{m3}\), \(A\sim G\), and \(n_{1}\sim n_{6}\)) are so determined that one gets a good description of the on-site energy in several reference structures with wide range of values of ϱ and \(r_{s}\). Therefore, the present formulation of TE-TB method can be a transferable scheme of calculating the energy bands, consistent with the total energy in several different environments or dynamical processes.
4. Numerical Procedure for Optimization of TE-TB Parameters
Crystalline materials can exist in several different phases or structures. For instance, silicon (Si) crystallizes to the diamond structure at low temperature under the atmospheric pressure and to the β-tin structure under high pressure. All TE-TB parameters would be determined to reproduce the band structures and the total energies in as many crystalline structures as possible. Our developed program package, named TB-PAC (TB-PAparametrization for Crystals), determines the STO parameter values, using the downhill simplex method.27) Then the local energy parameters are determined by Levenberg–Marquardt method of the nonlinear least-squares.28)
First, we should choose the principal structure S0 (e.g., the diamond structure in the case of Si). Then we adopt the procedure of machine-determination of the TE-TB parameters as following:
Step1
For the principal structure S0 with a fixed lattice constant (an effective atomic density \(\rho^{0}\)), determine values of STO parameters \((\zeta_{1},\zeta_{2}, C_{(1)}/C_{(2)}, a_{al}+b_{al}(\rho^{0})^{2/3}+c_{al}(\rho^{0})^{4/3}+d_{al}(\rho^{0})^{6/3})\) so as to be in a good agreement between the first principles (FP) calculation and the present TE-TB calculation for E–\(\boldsymbol{{k}}\) relation \(E_{n}(\boldsymbol{{k}})\) in this fixed \(\rho^{0}\).
Step2
For the same structure S0 with changing the lattice constant (changing ρ), calculate the E–\(\boldsymbol{{k}}\) relation and the band energy \(E_{\text{band}}\). Then determine \(\lambda^{\rho}\), \(a_{al}\), \(b_{al}\), \(c_{al}\), and \(d_{al}\), so as to get a good agreement of the band structures and the band energy between FT and TE-TB; \begin{align*} E_{n}(\boldsymbol{{k}})^{\text{TE-TB}}|_{\text{S}_{0}} &\simeq E_{n}(\boldsymbol{{k}})^{\text{FP}}|_{\text{S}_{0}}\\ E_{\text{band}}^{\text{TE-TB}}|_{\text{S}_{0}} &\simeq E_{\text{band}}^{\text{FP}}|_{\text{S}_{0}}. \end{align*} If one could not see a good agreement, return to Step 1 and choose another set of parameters.
Step3
Calculate E–\(\boldsymbol{{k}}\) relation and the band energy \(E_{\text{band}}\) for different crystalline structures S1, S2, S3, … and check the agreement between FP and TE-TB; \begin{align*} E_{n}(\boldsymbol{{k}})^{\text{TE-TB}}|_{\text{S}_{i}} &\simeq E_{n}(\boldsymbol{{k}})^{\text{FP}}|_{\text{S}_{i}}\\ E_{\text{band}}^{\text{TE-TB}}|_{\text{S}_{i}} &\simeq E_{\text{band}}^{\text{FP}}|_{\text{S}_{i}}\quad i = 1, 2,3,\ldots. \end{align*} Here we introduce \(\Delta(x)\)-term in Eq. (19) for representing characteristics of different local environment in several different crystalline structures. If the agreement is not satisfactory, choose another values of STO parameters, \(\lambda^{\rho}\), \(a_{al}\), \(b_{al}\), \(c_{al}\), and \(d_{al}\), and return to Step 2.
Step4
Calculate \(E_{\text{lattice}}^{\text{TE-TB}}\) in Eq. (5) for structures S0, S1, S2, S3, … with the resultant electron charge density obtained in Steps 1–3.
Step5
Using the calculated \(E_{\text{band}}^{\text{TE-TB}}\) and \(E_{\text{lattice}}^{\text{TE-TB}}\), evaluate the on-site term in TE-TB as \begin{equation} E_{\text{on-site}}^{\text{TE-TB}}|_{\text{S}_{i}} = \{E_{\text{total}}^{\text{FP}} - E_{\text{band}}^{\text{TE-TB}} -E_{\text{lattice}}^{\text{TE-TB}}\}|_{\text{S}_{i}}. \end{equation} (31) Then, using Eqs. (24)–(30), one should try to find a set of local energy parameters to satisfy \begin{equation} E_{\text{on-site}}(r_{s},\varrho)|_{\text{S}_{i}} = E_{\text{on-site}}^{\text{TE-TB}}|_{\text{S}_{i}}\quad (i = 0,1,2, \ldots). \end{equation} (32)
Step6
Calculate \begin{equation} E_{\text{total}}^{\text{TE-TB}} |_{\text{S}_{i}} \equiv \{E_{\text{band}}^{\text{TE-TB}} + E_{\text{lattice}}^{\text{TE-TB}} +E_{\text{on-site}}(r_{s},\varrho)\}|_{\text{S}_{i}} \end{equation} (33) for S0, S1, S2, S3, …. Because of the definition Eq. (31) in Step 5, one expects a good agreement for \(E_{\text{total}}\); \begin{equation} E_{\text{total}}^{\text{TE-TB}}|_{\text{S}_{i}} \simeq E_{\text{total}}^{\text{FP}}|_{\text{S}_{i}}\quad i = 0, 1, 2, 3,\ldots. \end{equation} (34)
Using above Steps 1–6, one obtains a good set of the TE-TB parameters. Therefore, if one chooses reference structures properly, one can expect a good estimate of electronic structures and inter-atomic forces for any arbitrary atomic configuration.
5. Results and Discussions: Example of Crystalline Si
In this section, we apply the present TE-TB method to crystalline Si with adopting ABINIT29) as the FP band structure calculation method. The diamond structure is chosen as the principal structure S0, and fcc (face centered cubic) and sc (simple cubic) structures as reference ones. The results will be tested in β-tin structure.
Adjusted band structures and band energy
For tetrahedral bonded semiconductors, one may be interested in the s, p bands. Details of s, p occupied bands are sensitive to the existence of atomic d orbitals. So we include the d-orbitals and the number of bands for adjustment is set to be up to p-bands (4 bands per atom). The off-set energy is set to be 20.0 eV.14) The \(\boldsymbol{{k}}\)-points are chosen, on 5–6 symmetry lines, totally 20–30 \(\boldsymbol{{k}}\)-points for each structure. The Hubbard parameter is chosen to be \(U_{aa}=1.0\) eV.
The resultant values of the STO parameters \(\zeta_{1}\), \(\zeta_{2}\), \(C_{(1)}/C_{(2)}\), and \(a_{al}\sim d_{al}\) are summarized in Table I. The parameters \(\zeta_{1}\), \(\zeta_{2}\) and \(C_{(1)}/C_{(2)}\) were fixed after the calculation of the diamond structure (S0) of the equilibrium lattice constant. The parameters, \(\lambda^{\rho}\), \(a_{al}\sim d_{al}\) are so determined to get overall agreement in reference structures, such as diamond, fcc and sc structures. The agreement between results of the first principle calculation and the present one is excellent, but a small discrepancy of 0.5 eV order would remain in the band structure for several different structures other than the diamond structure. This small discrepancy causes a conspicuous difference in \(E_{\text{band}}\), since the total energy difference between different crystalline structures is of order of 0.1 eV/atom. Then we introduce Δ-term in Eq. (19) \begin{align} \Delta(x) &= 0.0115[1+\tanh\{38(x-0.227)\}]\quad \notag\\ &\quad\text{(in atomic unit)}. \end{align} (35)
Table I. STO parameters of crystalline Si in atomic unit. \(a_{\Delta}=0.0115\), \(b_{\Delta}=38.0\), \(x_{\Delta}=0.227\) (in atomic unit).
The band structures of the diamond structure are shown, with changing the lattice constant, in Fig. 1 and those of fcc, and sc in Figs. 2 and 3, compared with those of FP results. Then, the band structures of β-tin structure are shown in Fig. 4, calculated by the same STO parameters. The β-tin structure is just the compressed diamond structure along c-axis, and in Si case \(c/a=0.552\).30) We assure that the agreement is excellent and the resultant STO parameter set can be said “universal” or “transferable”.
Figure 1. (Color online) Band structure of Si in the diamond structure. The red dotted lines: by the first principles calculation (ABINIT), and the black solid lines: by present TE-TB method. From the left top to right bottom, the lattice constant changes (a) 4.616 Å (15% compression), (b) 4.900 Å (10% compression), (c) 5.159 Å (5% compression), (d) 5.430 Å (equilibrium), and (e) 5.7015 Å (5% dilatation). The dotted lines: by the first principles calculation (ABINIT), and the bold solid lines: by present TE-TB method. Notice the scale difference on the ordinate and the change of the whole band width.
Figure 2. (Color online) Band structure of Si in fcc structure. The red dotted lines: by the first principles calculation (ABINIT), and the black solid lines: by present TE-TB method. The lattice constant changes; (a) 3.6 Å (5% compression), (b) 3.8 Å (equilibrium), (c) 3.9 Å (5% dilatation), and (d) 4.2 Å (10% dilatation).
Figure 3. (Color online) Band structure of Si in sc (simple cubic) structure. The red dotted lines: by the first principles calculation (ABINIT), and the black solid lines: by present TE-TB method. The lattice constant changes; (a) 2.4 Å (5% compression), (b) 2.5 Å (equilibrium), (c) 2.6 Å (5% dilatation), and (d) 2.8 Å (10% dilatation).
Figure 4. (Color online) Band structure of Si in β-tin structure (\(c/a=0.552\)). The β-tin structure reduces to the diamond structure when \(c/a=\sqrt{2}\). The red dotted lines: by the first principles calculation (ABINIT), and the black solid lines: by present TE-TB method. The lattice constant changes; (a) 4.3 Å (10% compression), (b) 4.6 Å (4% compression), (c) 4.8 Å (equilibrium), and (d) 5.1 Å (6% dilatation).
The band energies \(E_{\text{band}}^{\text{TE-TB}}\) of the present TE-TB method are shown and compared with that \(E_{\text{band}}^{\text{FP}}\) of the reference FP method in Fig. 5. The slope of the present calculation is slightly less steeper in common. This may suggest that \(R_{ab}\) dependence of K in Eq. (14) would be efficiently important (Wolfsberg–Helmholtz constant) and it remains unconcluded issue at present. Since the agreement is so good in the band structures, the agreement of the band energy is excellent.
Figure 5. (Color online) Band energies of Si of diamond, fcc, sc, and β-tin structures as a function of n.n. distance.
Lattice energy
Once a set of STO parameters is fixed and the local charge density is calculated, one can calculate the lattice energy Eq. (5) with the help of the expressions in Appendix A. The evaluated lattice energies are shown in Fig. 6. Since Si atoms in the simple crystalline structures are electrically neutral, the lattice energy is small. This term might be much larger in the case of larger charge transfer in ionic compounds.
Figure 6. (Color online) Lattice energy for Si diamond, fcc, sc, and β-tin structures.
On-site energy
There remains another term in the total energy, e.g., the on-site energy which is the sum of the exchange–correlation energy and the on-site electron–electron electrostatic energy, not included in the band energy. The on-site energy \(E_{\text{on-site}}\) is, in the present TE-TB method, defined by the value of (the total energy of FP calculation) minus (the band energy plus lattice energy of the present TE-TB method) as Eq. (31). However, it must be noticed that, when we compare the on-site energy with that of the FP-calculation, a certain term in the Coulomb interaction may be included in the “on-site” energy. This issue is discussed in more detail in Appendix B.
The evaluated local energy parameters are summarized in Table II. In order to guarantee matching of calculated values by the parameter determination using Levenberg–Marquardt method, it is necessary to avoid cancellation of significant digits in computation and we give these values in double precision in the Table II. In the ϱ–\(r_{s}\) space only a narrow range is used, so uniqueness of parameters may be lost. This may be one reason why the local energy parameters in Table II show large values.
Table II. Local energy parameters of crystalline Si (atomic unit).
The calculated on-site energy is then depicted in Figs. 7(a) and 7(b). One can conclude that the local energy parameters can be determined so as to reproduce the on-site energy. One must notice that we could not assure the applicability, in more wide range of \(\varrho^{a}\) and \(r_{s}^{a}\), of the form of the on-site energy Eqs. (28)–(30), because our scheme is a kind of interpolation but not extrapolation. As shown in Fig. 7(b), the range of the ϱ–\(r_{s}\) curve for interpolation (actually between the diamond and fcc structures) is limited and the result could not be extremely guaranteed outside. To extend the ϱ–\(r_{s}\) region, we need to prepare reference structures further.
Figure 7. (Color online) (a) On-site energy for Si diamond, fcc, sc, and β-tin structures, as a function n.n. distance. (b) On-site energy Eq. (25) as a function of ϱ and \(r_{s}\). Black points and thin black line: diamond structure, green points and chain line: fcc, blue points and dotted line: simple cubic and red points and chain line: β-tin structure. It may be noticed that the points of fcc, sc, and β-tin structures follow almost the same curves in ϱ–\(r_{s}\) space. The numbers put on contour lines in the figure are in eV unit.
Total energy
All the partial energies, the band energy, the lattice energy and the on-site energy, have been determined so far. Then we can estimate the total energies of the TE-TB method \(E_{\text{total}}^{\text{TE-TB}}\) Eq. (33), which are drawn in Fig. 8. Owing to the actual procedure, the agreement with that of the FP-calculations is excellent.
Figure 8. (Color online) Total energy for Si diamond, fcc, sc, and β-tin structures as functions of the average atomic density.
Test in β-tin structure of Si
The total energy of the β-tin structure is estimated with no adjustable parameters, and the reproduction is excellent. The origin of the apparent discrepancy in β-tin structure is the band energy in larger n.n. distance region and the on-site energy in shorter n.n. distance region. The discrepancy of the total energy of the β-tin structure is less than 0.1 eV/atom.
Though a further improvement may be necessary, one would conclude that the TE-TB parameter set is transferable over a variety of arbitrary structures.
Dynamical stability and effects of CSC-term
In the case of monotonic crystalline solids with changing the lattice constants, the resultant electron charges (Mullikan charge) for all atoms/ions in TE-TB calculation do not change. Then they should be the same as those of the starting ones for constructing the TB Hamiltonian, which are the results of the first principles calculation and the deviation of charges (\(\Delta q_{a}\)) should be zero. Therefore, the CSC-term does play no role.
When one distorts the structure or tries to ascertain the dynamical stability, the CSC-term starts to work. In some calculations, the stability of the crystalline structure cannot be kept without the CSC-terms. In the present TE-TB formulation, the band energy of the empty band is just zero and the on-site energy may, sometime, be lower when the atomic density would become higher. In this case, once an atom looses a certain amount of electrons, then the system does not have a mechanism to restore electrons back to the original ion.
In actual physical situation, when an atom looses an electron, the ionic potential becomes deeper by an affinity energy, and the atom attracts an electron back to the original ion. On the contrary, when an atom gets an excess electron, the atomic potential becomes higher by an ionization energy, and the atom releases an electron to neighboring ions. The ρ-dependence of \(\varepsilon_{al}^{\text{atom}}\) cannot be responsible for this charge restoration and the CSC-term should be used.
6. Discussion and Conclusion
We have formulated a novel TB method based on the band structure and the total energy calculated by the LDA. The total energy is divided into the band energy, the lattice energy, and the on-site energy (the exchange–correlation energy and electron–electron electro-static energy not included in the band energy). Each term is parametrized (the TE-TB parameters), based on the individual physical origin and the environmental dependence. The functional form of each term is physically meaningful, and the force acting on each atom can be calculated without time-consumption for the MD simulation.
The present formulation may be promising for further achievement of quantum molecular-dynamics simulation, which is now studying in several examples of various materials. Attempts have already been made to apply the present method to crystal GaAs and SiC, in which case the parameter values of the single substance can be good starting values for the compound. Further research on a more comprehensive approach to extend the applicable areas of the parameters ϱ and \(r_{s}\) (including selection of other parameters) is needed. Introduction of Machine Learning and other techniques that can prepare various artificial reference structures is one of future possibilities.
We have now prepared two basic tools, the solver of linear algebra (i.e., the Krylov subspace method)21,22) and the present novel TE-TB method. The program package ELSES (Extra Large Scale Electronic Structure calculations) with the present automatic parameter determination program (TB-Parametrization in Crystals = TB-PAC) is now being prepared to be open.31) The similar parameter determination procedure would be implemented also in TB-PAM.
Acknowledgements
This work is partly supported by Adaptable and Seamless Techronology transfer Program through target-driven R&D (ASTEP project) of Japan Science and Technology Agency (JST). We also thank all members of ELSES Consortium for their kind and warm support.
Appendix A:
Calculation of the Two-center Coulomb Integral
The inter-atomic electron–electron interaction in the case of \(a\neq b\) is called two-center Coulomb integral: \begin{equation} \int d\boldsymbol{{r}}\int \mathrm{d}\boldsymbol{{r}}^{\prime} \frac{n_{a}(\boldsymbol{{r}})n_{b}(\boldsymbol{{r}}^{\prime})}{|\boldsymbol{{r}}-\boldsymbol{{r}}^{\prime}|}. \end{equation} (A·1) We show, in this Appendix A, how to calculate Eq. (A·1) faithfully as much as possible. The inter-atomic electron–electron interaction with an arbitrary angular momentum component can be exactly calculated for STO but it would be time-consuming in the repeated calculation in MD simulations. Then, we set two assumptions; (1) the atomic charge distribution to be spherical and (2) the STO charge density to be replaced by a Gaussian-type charge density. Refer to Ref. 32 which gives explicit procedure and expressions for calculation of the two-center Coulomb integrals.
Assumption of spherical charge
A spherical atomic charge distribution at an atom a is assumed as \begin{align} n_{a}(\boldsymbol{{r}}_{a}) &\simeq \sqrt{\frac{1}{4\pi}} n_{a}(r_{a})Y_{00}(\theta, \phi)\notag\\ &= \frac{1}{4\pi} \{\rho_{a}^{\text{s}}(r_{a}) N_{a}^{s} + \rho_{a}^{p}(r_{a}) N_{a}^{\text{p}} + \rho_{a}^{d}(r_{a}) N_{a}^{\text{d}}\}, \end{align} (A·2) where \(\rho_{a}^{l} (r_{a}) = |R_{nl}^{(a)}(r_{a})|^{2}\) and \(\int_{0}^{\infty}\rho_{a}^{l}(r_{a}) r_{a}^{2}\,\mathrm{d}r_{a} =1\) and we have \begin{equation} \int \mathrm{d}\boldsymbol{{r}}_{a}\,n_{a}(\boldsymbol{{r}}_{a}) = N_{a}^{\text{s}} + N_{a}^{\text{p}} + N_{a}^{\text{d}}, \end{equation} (A·3) where \(N_{a}^{l}\) is the electron occupation number of atomic orbital l (l = s, p, or d).
Let us write a single-ζ STO as \begin{equation} (anlm) \equiv \phi_{nlm}^{(a)}(r) = \sqrt{\frac{(2\zeta_{l}^{a})^{2n+1}}{(2n)!}}r^{n-1}\mathrm{e}^{-\zeta_{l}^{a} r} Y_{l}^{m}(\theta, \phi), \end{equation} (A·4) and also we define \begin{equation} (n,\zeta_{l}) \equiv \phi_{n00}^{(a)}(r) = \sqrt{\frac{(2\zeta_{l})^{2n+1}}{(2n)!}} r^{n-1}\mathrm{e}^{-\zeta_{l} r}Y_{00}. \end{equation} (A·5) Since we are interested only in the spherical distribution of charges, or the charge distribution is assumed to be of spherical symmetry in Eq. (A·2), we define the following charge distribution, \begin{align} (n,\zeta)(n,\zeta^{\prime}) &= \left(\frac{2\zeta}{\zeta+\zeta'}\right)^{n+1/2}\left(\frac{2\zeta'}{\zeta+\zeta'}\right)^{n+1/2} [2n-1,\bar{\zeta}], \end{align} (A·6) \begin{align} [N,\bar{\zeta}] &= \sqrt{\frac{1}{4\pi}} \frac{(2\bar{\zeta})^{N+2}}{(N+1)!} r^{N-1}\mathrm{e}^{-2\bar{\zeta}r}Y_{00}(\theta,\phi), \end{align} (A·7) where \begin{equation*} \bar{\zeta} = (\zeta+\zeta^{\prime})/2, \end{equation*} and \([N,\bar{\zeta}]\) is normalized as \(\int\mathrm{d}\boldsymbol{{r}}[N,\bar{\zeta}] =1\).
In the case of a double-ζ STO, \((n;\zeta_{1},\zeta_{2}) =C_{(1)}(n,\zeta_{1})+C_{(2)}(n,\zeta_{2})\), we have 3-terms of normalized STO charge \([\cdots]\) as \begin{align} |(n;\zeta_{1},\zeta_{2})|^{2} &= C_{(1)}^{2}[2n-1,\zeta_{1}]+C_{(2)}^{2}[2n-1,\zeta_{2}]\notag\\ &\quad + 2C_{(1)}C_{(2)} \left(\frac{2\zeta}{\zeta+\zeta^{\prime}}\right)^{n+1/2} \left(\frac{2\zeta^{\prime}}{\zeta+\zeta^{\prime}}\right)^{n+1/2} \notag\\ &\quad \times\left[2n-1,\frac{\zeta+\zeta'}{2}\right]. \end{align} (A·8) Each of above three terms of the STO charge density is replaced by a Gaussian-type charge density.
Replacement of STO charge density by Gaussian-type spherical charge density
Let us replace the STO charge density \([2n_{S}-1,\bar{\zeta}]\) at the atom a by a Gaussian charge density \(n_{a}^{G}(n_{G},\alpha;\boldsymbol{{r}})\); \begin{align} n_{a}^{G}(n_{G},\alpha;\boldsymbol{{r}}) &= \frac{1}{4\pi}N_{n_{G}}^{G}(\alpha)^{2}\,r^{2(n_{G}-1)}\mathrm{e}^{-2\alpha r^{2}}, \end{align} (A·9) \begin{align} N_{n_{G}}^{G}(\alpha)^{2} &= \frac{2^{n_{G}+1}}{(2n_{G}-1)!!}\sqrt{\frac{(2\alpha)^{2n_{G}+1}}{\pi}},\notag \end{align} where \(\int\mathrm{d}\boldsymbol{{r}}\,n_{a}^{G}(n_{G},\alpha;\boldsymbol{{r}}) = 1\). Here we set \begin{equation} n_{S} = n_{G} \end{equation} (A·10) and α to maximize the overlap of \(\sqrt{[2n_{S}-1,\bar{\zeta}]}\) and \(\sqrt{n^{G}(n_{G};\alpha;\boldsymbol{{r}})}\), which gives rise an expression (\(n=n_{S}=n_{G}\)) \begin{equation} \alpha = \cfrac{n+\cfrac{1}{2}}{2} \cdot \frac{\displaystyle\int_{0}^{\infty} \mathrm{d}r\,r^{2n} \mathrm{e}^{-\bar{\zeta}r-\alpha r^{2}}}{\displaystyle\int_{0}^{\infty} \mathrm{d}r\,r^{2n+2} \mathrm{e}^{-\bar{\zeta}r-\alpha r^{2}}}. \end{equation} (A·11) Since we assume the spherical Gaussian charge density for local charge distribution, the electron density can be expanded as \begin{equation} n(\boldsymbol{{r}}) = \sum_{a} n_{a}(\boldsymbol{{r}}) = \sum_{a\kappa} D_{a}^{\kappa} n_{a}^{G}(n_{a},\alpha_{a};\boldsymbol{{r}}-\boldsymbol{{R}}_{a}), \end{equation} (A·12) where \(D_{a}^{\kappa}\)'s have been determined in the band structure calculation in the present TE-TB scheme.
Two-center Coulomb integrals
For the calculation of the two-center Coulomb integral, we just follow the procedure given in Ref. 32, where the Gaussian-type orbitals (GTO) are defined as \begin{equation} \chi_{nil}(\alpha,\boldsymbol{{r}}) = N^{\chi}_{nl}(\alpha)\,\mathrm{e}^{-\alpha r^{2}}\,r^{2n+l}Y_{Lem}(\hat{\boldsymbol{{r}}}) \end{equation} (A·13) and \begin{equation*} \{N^{\chi}_{nl}(\alpha)\}^{2} = \left[\frac{\{2(2n+l+1)-1\}!!}{2^{(2n+l+1)+1}} \sqrt{\frac{\pi}{(2\alpha)^{2(2n+l+1)+1}}}\right]^{-1}. \end{equation*} Our Gaussian-type charge distribution (A·9) is expressed as a product of spherical GTO's χ as \begin{equation} n_{a}^{G}(n_{a}, \alpha_{a}; \boldsymbol{{r}}) = C^{G/S}(n)\,\chi_{0,0,0}(\alpha_{a}, \boldsymbol{{r}}_{a}) \chi_{n_{a}-1,0,0}(\alpha_{a}, \boldsymbol{{r}}_{a}), \end{equation} (A·14) then the two-center Coulomb integral can be rewritten as \begin{align} &\int d\boldsymbol{{r}}_{1}\int\mathrm{d}\boldsymbol{{r}}_{2} \frac{n_{a}^{G}(n_{a},\alpha_{a};\boldsymbol{{r}}_{1})n_{b}^{G}(n_{b},\alpha_{b};\boldsymbol{{r}}_{2})}{|\boldsymbol{{r}}_{1}-\boldsymbol{{r}}_{2}|}\notag\\ &\quad = C^{G/S}(n_{a})\,C^{G/S}(n_{b}) \langle \chi_{0,0,0}(\alpha_{b}, \boldsymbol{{r}}_{1})) \chi_{0,0,0}(\alpha_{a}, \boldsymbol{{r}}_{2})/|r_{12}|\notag\\ &\qquad \times \chi_{n_{b}-1,0,0}(\alpha_{b}, \boldsymbol{{r}}_{1}) \chi_{n_{a}-1,0,0}(\alpha_{a}, \boldsymbol{{r}}_{2}) \rangle. \end{align} (A·15) Following the procedure in Ref. 32, we can get the final expression as \begin{align} &\int d\boldsymbol{{r}}\int\mathrm{d}\boldsymbol{{r}}^{\prime} \frac{n_{a}^{G}(n_{a},\alpha_{a};\boldsymbol{{r}})n_{b}^{G}(n_{b},\alpha_{b};\boldsymbol{{r}}^{\prime})}{|\boldsymbol{{r}}-\boldsymbol{{r}}^{\prime}|}\notag\\ &\quad = \frac{1}{R_{ab}} \cdot \frac{2}{\sqrt{\pi}}\int_{0}^{\sqrt{\zeta_{ab}}R_{ab}}\mathrm{e}^{-t^{2}}\,\mathrm{d}t\notag\\ &\qquad - \sqrt{\frac{2}{\pi}}\cdot\frac{\{2(n_{a}-1)\}!!\{2(n_{b}-1)\}!!}{(2n_{a}-1)!!(2n_{b}-1)!!}\notag\\ &\qquad \times \sum_{k(\neq 0)}^{(n_{a}+n_{b}-2)} C(n_{a},\alpha_{a};n_{b},\alpha_{b};k)(2\zeta_{ab})^{k+1/2}\notag\\ &\qquad \times \mathrm{e}^{-\zeta_{ab} R_{ab}^{2}} H_{2k-1}(\sqrt{2\zeta_{ab}}R_{ab})/(\sqrt{2\zeta_{ab}}R_{ab}), \end{align} (A·16) where \begin{equation*} \zeta_{ab} = 2\alpha_{a}\alpha_{b}/(\alpha_{a}+\alpha_{b}), \end{equation*} and \begin{align*} &C(n_{a},\alpha_{a};n_{b},\alpha_{b};k)\\ &\quad = \sum_{j=\textit{Max}(0,k-(n_{b}-1))}^{\textit{Min}(n_{a}-1,k)} \frac{\{1/(8\alpha_{a})\}^{j}\{1/(8\alpha_{b})\}^{k-j}}{j!(k-j)!}\\ &\qquad \cdot \frac{\Gamma(n_{a}+3/2)}{(n_{a}-j)! \Gamma(j+3/2)} \cdot \frac{\Gamma(n_{b}+3/2)}{(n_{b}-k+j)!\Gamma(k-j+3/2)}. \end{align*} The function \(H_{2k-1}(x)\) is an odd-order Hermite polynomial defined as \begin{equation*} H_{2n-1}(x) = \sum_{r=0}^{n-1} (-1)^{r}(2r-1)!! \begin{pmatrix} 2n-1\\ 2r \end{pmatrix}x^{2n-1-2r}. \end{equation*}
When we set the limit \(R_{ab}\rightarrow\infty\) in Eq. (A·16), the second term vanishes and the first term remains as \begin{equation*} \int d\boldsymbol{{r}}\int\mathrm{d}\boldsymbol{{r}}^{\prime} \frac{n_{a}^{G}(n_{a},\alpha_{a};\boldsymbol{{r}})n_{b}^{G}(n_{b},\alpha_{b};\boldsymbol{{r}}^{\prime})}{|\boldsymbol{{r}}-\boldsymbol{{r}}^{\prime}|} \rightarrow \frac{1}{R_{ab}}, \end{equation*} which is the rigorous limiting form.
Appendix B:
First Principles Pseudo Potential Scheme for the Reference Calculation
The LDA total energy of a system, with several atoms in a unit cell, can be given in the framework of the first principles pseudo-potential theory as follows:33) \begin{align} E_{\text{Total}}^{\text{LDA}} &= \sum_{i} \varepsilon_{i}(V_{\textit{Coul}}(\mathbf{0}) = 0, U_{ps}^{\mu}(\mathbf{0}) = \alpha_{1}^{\mu}) \end{align} (B·1) \begin{align} &\quad + \Omega \left(\frac{\bar{Z}}{\Omega_{c}}\right) 4\pi \beta \end{align} (B·2) \begin{align} &\quad + \left[E_{xc} -\int\mathrm{d}\boldsymbol{{r}}\,\mu_{xc}(\boldsymbol{{r}})\rho(\boldsymbol{{r}}) \right] \notag\\ &\quad- \frac{1}{2}\sum_{a} \iint \mathrm{d}\boldsymbol{{r}}\,\mathrm{d}\boldsymbol{{r}}^{\prime} \frac{\rho_{a}(\boldsymbol{{r}})\rho_{a}(\boldsymbol{{r}}^{\prime})}{|\boldsymbol{{r}}-\boldsymbol{{r}}^{\prime} |} \end{align} (B·3) \begin{align} &\quad + \left\{-\frac{1}{2}\sum_{a \neq b} \iint\mathrm{d}\boldsymbol{{r}}\,\mathrm{d}\boldsymbol{{r}}^{\prime} \frac{\rho_{a}(\boldsymbol{{r}})\rho_{b}(\boldsymbol{{r}}^{\prime})}{|\boldsymbol{{r}}-\boldsymbol{{r}}^{\prime} |} + \frac{N_{\textit{cell}}^{2}}{2\Omega}\int \frac{\bar{Z}^{{2}}}{r}\,\mathrm{d}\boldsymbol{{r}}\right\} \end{align} (B·4) \begin{align} &\quad + \left\{\frac{1}{2}\sum_{a \neq b}\frac{Z^{2}}{R_{ab}} -\frac{N_{\textit{cell}}^{2}}{2\Omega}\int \frac{\bar{Z}^{{2}}}{r}\,\mathrm{d}\boldsymbol{{r}}\right\}, \end{align} (B·5) where \(\bar{Z}\) is the sum of effective core charges of ions in the unit cell, \(\Omega_{c}\) is the volume of the unit cell, Ω is the total volume of the system, and \(N_{\textit{cell}}=\Omega/\Omega_{c}\) is the total number of unit cells in the system. Then first term (B·1) is the band energy, and \(\alpha_{1}^{\mu}\) is the remaining constant value value in the pseudopotential. The third line (B·3) is just our \(E_{\text{on-site}}\). Here the charge density in the reciprocal space is expressed as \begin{equation} \rho(\boldsymbol{{G}}) = \frac{\bar{Z}}{\Omega_{c}} +\beta G^{2} + O(G^{4}), \end{equation} (B·6) and β in (B·2) and (B·6) is defined as \begin{equation} 4\pi \beta = \frac{1}{\Omega} \iint\mathrm{d}\boldsymbol{{r}}\,\mathrm{d}\boldsymbol{{r}}^{\prime} \cdot \frac{\rho(\boldsymbol{{r}})-\bar{Z}/\Omega_{c}}{| \boldsymbol{{r}}-\boldsymbol{{r}}^{\prime} |}. \end{equation} (B·7) The sum of the forth and fifth terms, (B·4) and (B·5), are our lattice energy \(E_{\text{lattice}}\).
The term (B·5) is usually called lattice ion–ion energy and may be calculated by, for instance, the Ewald method. The sum of divergent terms with long wave-length components (\(\boldsymbol{{G}}\rightarrow 0\)) in the electron–electron Coulomb interaction, the sum of the pseudo-potential and the lattice ion–ion term results in zero. And we introduce terms, such as \(\int (1/r)\,\mathrm{d}\boldsymbol{{r}}\), to cancel the diverging components in the expression of the total energy above. The β-term appears in the Coulomb term in the band energy and the electron–electron Coulomb interaction in (B·4) and they cancel with each other. Then, they do not appear explicitly in the usual first-principles pseudopotential theory. However, in the present framework, we treat Eqs. (B·4) and (B·5) as one unit as \(E_{\text{lattice}}\). Thus the β-term (B·2) is combined with other contributions of the exchange–correlation energy in Eq. (B·3) to be an on-site energy term.
References
1 R. Hoffmann, J. Chem. Phys. 39, 1397 (1963). 10.1063/1.1734456 Crossref, Google Scholar
2 R. Hoffmann, J. Chem. Phys. 40, 2474 (1964). 10.1063/1.1725550 Crossref, Google Scholar
3 R. Hoffmann, J. Chem. Phys. 40, 2480 (1964). 10.1063/1.1725551 Crossref, Google Scholar
4 R. Hoffmann, J. Chem. Phys. 40, 2745 (1964). 10.1063/1.1725601 Crossref, Google Scholar
5 M. Wolfsberg and L. Helmholz, J. Chem. Phys. 20, 837 (1952). 10.1063/1.1700580 Crossref, Google Scholar
6 G. Calzafferri and R. Rytz, J. Phys. Chem. 100, 11122 (1996). 10.1021/jp960840t Crossref, Google Scholar
7 A. B. Anderson, J. Chem. Phys. 62, 1187 (1975). 10.1063/1.430562 Crossref, Google Scholar
8 K. Nath and A. B. Anderson, Phys. Rev. B 41, 5652 (1990). 10.1103/PhysRevB.41.5652 Crossref, Google Scholar
9 D. Porezag, Th. Fauenheim, Th. Köhler, G. Seifert, and R. Kaschner, Phys. Rev. B 51, 12947 (1995). 10.1103/PhysRevB.51.12947 Crossref, Google Scholar
10 M. J. Mehl and D. A. Papaconstantopoulos, Phys. Rev. B 54, 4519 (1996). 10.1103/PhysRevB.54.4519 Crossref, Google Scholar
11 S. H. Yang, M. J. Mehl, and D. A. Papaconstantopoulos, Phys. Rev. B 57, R2013 (1998). 10.1103/PhysRevB.57.R2013 Crossref, Google Scholar
12 D. A. Papaconstantopoulos and M. J. Mehl, J. Phys.: Condens. Matter 15, R413 (2003). 10.1088/0953-8984/15/10/201 Crossref, Google Scholar
13 M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys. Rev. B 58, 7260 (1998). 10.1103/PhysRevB.58.7260 Crossref, Google Scholar
14 J. Cerdá and F. Soria, Phys. Rev. B 61, 7965 (2000). 10.1103/PhysRevB.61.7965 Crossref, Google Scholar
15 T. Hoshi, Y. Iguchi, and T. Fujiwara, Phys. Rev. B 72, 075323 (2005). 10.1103/PhysRevB.72.075323 Crossref, Google Scholar
16 Y. Iguchi, T. Hoshi, and T. Fujiwara, Phys. Rev. Lett. 99, 125507 (2007). 10.1103/PhysRevLett.99.125507 Crossref, Google Scholar
17 S. Nishino, T. Fujiwara, H. Yamasaki, S. Yamamoto, and T. Hoshi, Solid State Ionics 225, 22 (2012). 10.1016/j.ssi.2012.01.045 Crossref, Google Scholar
18 S. Nishino, T. Fujiwara, H. Yamasaki, S. Yamamoto, S. Yamamoto, and T. Hoshi, The 62nd Annual Meeting of the International Society of Electrochemistry, 2011. Google Scholar
19 S. Nishino and T. Fujiwara, J. Mol. Modeling 19, 2363 (2013). 10.1007/s00894-013-1767-2 Crossref, Google Scholar
20 S. Nishino, T. Fujiwara, and H. Yamasaki, Phys. Rev. B 90, 024303 (2014). 10.1103/PhysRevB.90.024303 Crossref, Google Scholar
21 H. Teng, T. Fujiwara, T. Hoshi, T. Sogabe, S.-L. Zhang, and S. Yamamoto, Phys. Rev. B 83, 165103 (2011). 10.1103/PhysRevB.83.165103 Crossref, Google Scholar
22 T. Hoshi, S. Yamamoto, T. Fujiwara, T. Sogabe, and S.-L. Zhang, J. Phys.: Condens. Matter 24, 165502 (2012). 10.1088/0953-8984/24/16/165502 Crossref, Google Scholar
23 S. Nishino, T. Fujiwara, N. Watanabe, and S. Yamamoto, J. Mol. Modeling 21, 169 (2015). 10.1007/s00894-015-2694-1 Crossref, Google Scholar
24 M. S. Daw and M. I. Baskes, Phys. Rev. B 29, 6443 (1984). 10.1103/PhysRevB.29.6443 Crossref, Google Scholar
25 W. Kohn and L. S. Sham, Phys. Rev. 140, A1133 (1965). 10.1103/PhysRev.140.A1133 Crossref, Google Scholar
26 A. C. Wahl, P. E. Cade, and C. C. J. Roothann, J. Chem. Phys. 41, 2578 (1964). 10.1063/1.1726326 Crossref, Google Scholar
27 J. A. Nelder and R. Mead, Comput. J. 7, 308 (1965). 10.1093/comjnl/7.4.308 Crossref, Google Scholar
28 K. Levenberg, Q. Appl. Math. 2, 164 (1944); 10.1090/qam/10666 Crossref;, Google ScholarD. Marquardt, SIAM J. Appl. Math. 11, 431 (1963). 10.1137/0111030 Crossref, Google Scholar
Total-energy Assisted Tight-binding Method Based on Local Density Approximation of Density Functional Theory
