The underestimation of high pressure in DFT+$U$ simulation for the wide range cold-pressure of lanthanide metals

All 15 kinds of lanthanide metals investigated in this work are the face centered cubic (FCC) structure, which is a conventional close-packed structure for lanthanide metals. Since the main purpose of this work is to reveal the capability of DFT+$U$ method to simulate metals with 4$f$ electrons under high-pressure, the other probable structures are not considered in this work. To include both low and ultra-high pressure, the volumes we considered are from 9 ${\rm A^{3}}$ to 28 ${\rm A^{3}}$. PAW and LAPW simulations are performed in VASPKresse1999 and WIEN2kBlaha2020 packages respectively, and we specially refer PAW (LAPW) results to VASP (WIEN2k) simulations in this work. The exchange-correlation functional used for all calculations is generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhof formation (PBE)Perdew1996-2 . A grid of $17\times 17\times 17$ is used for the k-spacing sampling. The pseudo-potentials used for PAW simulation are recommended by VASP, which means the 5$s$, 5$p$, 5$d$, 4$f$ and 6$s$ are treated as valence for all lanthanides. And the cutoff energy of plane wave is 600 eV, which is large enough for calculation of energy and force for all conditions in this work. The augmentation radius of all 4$f$ elements in LAPW calculation is 1.11 A (2.1 bohr.), and the $R_{\rm kmax}$. The “+$U$” correction takes the isotropic form to reduce the number of parameters. The $U_{\rm eff}$ (or the $U-J$ in some context) is 5.0 eV for all lanthanide elements, and the double-counting term is “fully localized limit” (FLL), which is the rational option for many DFT+$U$ worksJiang2009 ; Moree2018 . For comparison, hybrid functional calculations using HSE (Heyd-Scuseria-Ernzerhof) formation are performed in VASPHeyd2003 . Since spin-polarization can drastically change the electronic states between different methods and makes the comparison much complicated, no spin-polarization is considered in this work. For the calculation of pressure, each energy-volume curve is fitted by Birch-Murnaghan equation with 50 sampling points in the range from 9 ${\rm A^{3}}$ to 28 ${\rm A^{3}}$.

In the “+$U$” correction, besides the correlation strength $U_{\rm eff}$ and double-counting term, one of the most important factors is the definition of local Hilbert space which the model Hamiltonian is worked onRyee2018 . From the technical point of view, a set of local orbitals is often used to define this local subspace. By acting the projection operator $\hat{P}(\mathbf{r},\mathbf{r}^{\prime})$ composed by the local orbitals on the whole Hilbert space, the local density matrix can be obtainedWang2016 . There are three kinds of most common formations of the projection operator: spherical augmented schemes, atomic orbitals and special functions (e.g. Wannier function)Shick1999 ; ORegan2010 ; Bengone2000 ; Cococcioni2005 . As the simulation tools of materials based on augmented basis sets (e.g. VASP and WIEN2k) are widely used, many works of “+$U$” correction are carried out with spherical augmented projection operator. The details of projection operators in PAW and LAPW and how they act on the single electron orbitals to obtain local density matrix are described in the appendix part.

The results of lanthanide metals’ pressure at large (upper panel, 28 $\rm A^{3}$) and small (lower panel, 9 $\rm A^{3}$) volume in different simulation conditions are shown in Fig.1. The simulation condition we considered includes DFT (PBE, green bar) and “+$U$” correction (red bar) in PAW and LAPW methods, and the HSE06 hybrid functional (blue bar) based on PAW. HSE06 is used as a benchmark, considering hybrid functional is not directly dependent on the local orbitals as “+$U$” correction, and the level of correction of correlation effect is not lower than the “+$U$” correction. Two important features are exhibited in Fig.1. Firstly, an abnormal high-pressure softening of the results from DFT+$U$ in PAW can be observed. In high pressure conditions, the pressure of some heavy lanthanide metals is underestimated by 50$\%$ to 80$\%$ in DFT+$U$ simulation from PAW comparing with that from LAPW. Meanwhile, in low pressure conditions, the DFT+$U$ results between PAW and LAPW are similar. Secondly, the pressure calculated by DFT is higher than the correlation correction methods (both “+$U$” and HSE06) in some heavy lanthanide metals (Tm, Yb and Lu). It is contrary to the general experience of these correction methods that “+$U$” and hybrid functional often give a higher pressure comparing to DFT at the same volume, as can be observed in our results of low-pressure conditionWang2016 ; Perdew1996 ; Jiang2009 ; Moree2018 . For the first feature mentioned above, it can be found that the PAW pseudo-potential is not the reason of the abnormal softening of pressure from DFT+$U$ in PAW, since both in high and low pressure the DFT results from LAPW and PAW are similar. The “+$U$” correction also should not be accused for these abnormal results, since HSE06 and “+$U$” correction in LAPW give similar results in both low and high pressure. The answer to the unreasonable results of “+$U$” correction in PAW under high pressure seems only related to the implementation of “+$U$” correction in PAW. The main difference of implementation “+$U$” in PAW and LAPW is that the overlap between augmented sphere in PAW cannot be avoided, since the radius of the sphere is fixed when the pseudo-potential is created in PAW. Other results connecting with the abnormal behavior of PAW in high-pressure condition is shown in the following, and a detailed analysis is also presented. For second feature that the pressure described by DFT is higher than that of correlation correction in high-pressure conditions, it may relate to the theorical mechanism behind correlation correction methods other than the implementation problem, since the phenomenon exhibits in both HSE06 and “+$U$” simulation. We give our interpretation of the phenomenon in the discussion.

In this subsection, the “+$U$” results of total energy, pressure and modulus changing versus volume (simplified to E-V, P-V and B-V in the work respectively) are displayed in PAW and LAPW. The E-V, P-V and B-V results are calculated in the range from 9 $\rm A^{3}$ to 26 $\rm A^{3}$. We have checked the pressure from the Hellmann-Feynman (reading from the software report) and from energy difference, and the results are the same. It should be mentioned here, the zero point of the total energy of each system is set as the energy of 9 $\rm A^{3}$. For the simulation in this subsection, we display the results of La, Gd and Yb, since these three lanthanide metals represent three kinds of 4$f$-occupation conditions, the empty occupation La, half occupation Gd and full occupation Yb. The Fig.2 shows the E-V, P-V and B-V results of La, Gd, and Yb. For La, the E-V results from PAW always change fast comparing to LAPW with volume decreasing. This is obvious in the P-V results, for pressure can be seen as the first derivative of energy with respect to volume. In the P-V results of La, the pressure of PAW is always large than that of LAPW at every volume point. However, it can be observed that the difference of pressure between PAW and LAPW is decreasing with the decreasing of volume. In the B-V results (the second derivative of energy with respect to volume), the difference between PAW and LAPW is almost vanishing at the small volume. For half occupation system Gd, with the decrease of volume, the total energy of PAW increases faster than LAPW below 100 GPa. At higher pressure, results of LAPW increase faster. In the P-V results, curves of PAW and LAPW cross at about 100 GPa in Gd. These crossings happen at a much larger volume in the B-V results of Gd. For the nearly fully occupied Yb, the variation of energy versus volume in PAW and LAPW is similar in low pressure. In P-V results, the pressure of LAPW is larger than PAW results at every volume point, and at the volume of 9 $\rm A^{3}$, the difference can be as large as 150 to 200 GPa, which has a relative difference larger than 50$\%$. For the B-V results, PAW gives smaller values comparing with LAPW, and an abnormal phenomenon is observed that modulus from PAW starts to decrease with the decreasing of volume at high pressure.

The “+$U$” correction of electronic structures at low and high pressure is investigated in this subsection. We also choose La, Gd and Yb as representatives for empty, half-filling and full occupied $f$-shell. The band structures calculated by the “+$U$” correction in PAW and LAPW are shown in Fig.3. At low pressure, it can be observed that the electronic structures from PAW and LAPW are similar, especially the itinerant bands below the Fermi level. For the local $f$-bands, the results of PAW and LAPW show the difference about 0.1 eV. The PAW results give lower $f$-bands below the Fermi level than LAPW results. At high pressure, the difference of electronic structures calculated from PAW and LAPW becomes larger. Besides all bands becoming more dispersed, PAW still makes the $f$-bands far away from the Fermi level comparing to LAPW, but the detailed features become much more complicated. However, it seems no qualitative difference between PAW and LAPW electronic structures at high pressure.

To understand the abnormal pressure feature shown in the results of “+$U$” in the PAW method, we should analyze how “+$U$” correction makes an influence on the total energy. It is because the pressure reflects the changing of total energy versus volume. The total energy of DFT+$U$ is

In this work, all formations of “$+U$” correction are shown on a diagonalized local density matrix. On one hand, it can make the analysis more clear. On the other hand, for the examples in this work, the diagonal part gives the major contribution to the results. The double counting of FLL form is adopted, which is

Without the anisotropic effect, the total energy correction in “+$U$” is

Since in most of the “+$U$” works, the $U_{\rm eff}$ is constant, $n_{m}^{\sigma}$ the occupation of local orbitals, becomes the key factor that affects the total energy.

The $n$ should have a value between 0 and 1 from the definition of the occupation of local orbitals. It is because $n_{mm}$ is

$\psi_{i}$ is KS orbital, and the summation is over all occupied orbitals. For the completeness of the KS orbitals, all the KS orbitals can construct a unit operator

Thus, it can be written out like

So,

Because $\sum_{i\notin occ}{\left\langle m\middle|\psi_{i}\right\rangle\left\langle\psi_{i}\middle|m\right\rangle}\geq 0$, $n_{mm}$ should little than $\left\langle m\middle|m\right\rangle$. If $\left|m\right\rangle$ is a well-defined atomic orbital, it should satisfy the normalization condition that $\left\langle m\middle|m\right\rangle=1$, and we can make a conclusion that $n_{mm}\leq 1.0$.

For the n with value in range 0 to 1 and $U_{\rm eff}$ $>$ 0.0 (electrons repulse with each other), the “+$U$” correction with FLL double counting term gives an energy correction larger than 0.0. This correction becomes 0 only if n equals 1 or 0. This feature of “+$U$” correction has been reported in previous works. Another well-known feature is that lower orbital energy is obtained with occupation larger than 0.5Ylvisaker .

In our PAW simulations, the occupation of some local density matrix elements exceeds 1.0, and we investigate the origin and influence of the phenomenon through the variation of local density matrix with the pressure (or volume). $\Delta N_{>1}$ is defined as the summation of all local orbitals which is more than 1.0, and has the form

$n_{mm}$ are the diagonal elements of the local density matrix, and the summation of m from -3 to 3 means all 7 orbitals in 4$f$ shell (spin polarization is not considered).

As shown in Fig.4, except for some elements with lower occupation (like, La, Ce, Pr), the occupation of most lanthanide metals can show obvious exceeding over 1.0. The systems with occupation over 1.0 can be divided into two kinds. One is that the occupation of some orbitals is over 1.0 even at low pressure (Tb, Dy, Ho). These systems correspond to the first dipping of pressure in serial of lanthanide metals in high pressure results. In the other kind occupation only exceeds 1.0 when the pressure is high enough. No matter which kind of behavior, as we have discussed above, the occupation of local orbital exceeding 1.0 happens only in two conditions: The normalization of KS orbitals is not satisfied, or the modular square of local orbitals has problem. As the DFT results of PAW and LAPW are similar, the KS-orbital should be normalized correctly. We also have tested the modular square of KS-orbitals dominated by $f$-feature, and the normalization can be achieved at $10^{-}8$ level. To analyse the local orbitals, we split the modulus square of PAW orbitals into several parts. The modular square of the PAW orbitals writes

$\left|\phi_{s}\right\rangle$ and $\left|{\widetilde{\phi}}_{s}\right\rangle$ are the all-electron partial orbital and pseudo partial orbital of the same atom respectively. $\left|\psi_{i}\right\rangle$ and $\left|{\widetilde{\psi}}_{i}\right\rangle$ are all-electron orbital and pseudo-orbital and these two terms are connected by a projector $\left|{\widetilde{p}}_{s}\right\rangle$. In usual case, $\left|\phi_{s}\right\rangle$, $\left|{\widetilde{\phi}}_{s}\right\rangle$ and $\left|{\widetilde{p}}_{s}\right\rangle$ is fixed as a certain pseudo-potential is used, while $\left|\psi_{i}\right\rangle$ and $\left|{\widetilde{\psi}}_{i}\right\rangle$ are renewed at each self-consistent cycle. The subscript s stands for a collection of orbital momentum l and m, and the site of atom. Here we split it into three parts: integral of pseudo-orbital on plane wave

integral of pseudo atomic orbital

and integral of all-electron atomic orbital

With this splitting, the integral of PAW orbital can be written as

Considering the definition of the diagonal element of local density matrix $n_{mm}$,

If $n_{mm}$ is fully contributed by one occupied KS-orbital $\psi_{i}$, $n_{mm}=\mathcal{s}_{i,atom}^{lm}$. The interval part of the integral is ${\widetilde{\mathcal{S}}}_{i,pw}-{\widetilde{\mathcal{S}}}_{i,atom}$, and the atom part of the integral is $\mathcal{S}_{i,atom}$. The interval and atom part of the integral stand for the probability to find the electron on orbital in the interval and augmented region respectively.

As we have shown, the all-electron atomic orbital has a close relation to the local orbital. And we can investigate the normalization of local orbital by the modular square of all-electron atomic orbital. From the physical meaning of orbital, $\mathcal{S}_{i,atom}$ should be less than 1.0, otherwise there is something abnormal happened to the local orbital. We analyze this problem at high pressure in PAW method in Yb as an example, since the abnormal behavior in high-pressure simulation is obvious in the heavy lanthanides. It should be clarified that the occupation defined in Eq.4 is not the exact case in augmented methods. However, if a case with only one type of radial orbital for each shell is considered, the occupation defined in Eq.4 could be equivalent with the augmented methods.

Here we consider a simplified situation, where only one kind of radial part $|R_{l}\rangle$ is used for the shell $l$ in each site. In this way, it is easier to compare the difference between occupation matrix elements obtained from augmented methods and the ones from projecting onto a set of local orbital ${|m\rangle}$, as shown in Eq.4 The occupation matrix elements in augmented methods become,

$n_{mm^{\prime}}^{\rm aug}$ stands for an occupation matrix element in PAW, and a general form of it is provided in appendix. The occupation matrix elements obtained by projecting onto local orbitals are,

$n_{mm^{\prime}}^{\rm proj}$ is the occupation matrix element from the procedure defined in Eq.4, and the radial part of local orbital $|R_{l}\rangle$ takes the same function as the all-partial orbital in PAW. It can be clearly observed that if the radial parts of the shell has the same shape, and the local orbital is well normalized, $n_{mm^{\prime}}^{\rm aug}=n_{mm^{\prime}}^{\rm proj}$.

In the table 1, we show the different parts of orbital integral of Yb at low and high pressure at $\Gamma$ point of reciprocal space. As in the “+$U$” correction of Yb, almost all $f$-feature bands are below the Fermi energy, the $\mathcal{S}_{i,atom}$ and ${\widetilde{\mathcal{S}}}_{i,atom}$ are almost contributed by $\mathcal{s}_{i,atom}^{fm}$, and the 7 orbitals listed in the table are the ones with obvious $f$-feature. At low pressure, the integral results perform normally. A finite positive value shows for every part, and the summation of all parts equals 1.0. At high pressure, though the summation of all parts is still 1.0, there shows some abnormal features. The part of the augmented region is over 1.0, and the interval region has negative values. These features are consistent with the occupation of local orbitals over 1.0 at high pressure and imply a normalization problem of the local orbital.

The drastic changes of $\mathcal{S}_{i,atom}$ and ${\widetilde{\mathcal{S}}}_{i,atom}$ between low-pressure and high-pressure conditions implies a variation of the wave function radial part with pressure in the augmented region. In VASP, there are two predefined radial functions forming the basis to describe the f-shell. We show these radial functions in Fig.5. The blue solid line represents a radial function with typical 4$f$-feature (localized around core, one peak and decaying to zero at distance), and we labeled it by $F_{1}$. The red dashed line represents a radial function more divergent, we labeled it by $F_{2}$. The black vertical dashed line labeled by $R_{\rm c}$ stands for the radius of the augmented region. The sudden drop of $F_{2}$ is a truncation, and since it is out of the $R_{\rm c}$, we believe it doesn’t affect the results directly.

In Fig.6, the average of radial charge distribution of 4$f$ orbital in Yb over a spherical surface is shown at low (blue solid line) and high pressure (red dashed line). The $R_{1/2}^{H}$ is the half of the distance between the nearby atoms at high pressure, which stands for the radial not overlap. The radial charge can reflect the shape of the local orbital. It can be observed that when R is smaller than $R_{1/2}^{H}$, the charge distributions of low and high pressure are both similar to the shape of $F_{1}$. When in the range between $R_{1/2}^{H}$ and $R_{\rm c}$, the high-pressure result shows are divergent feature like $F_{2}$. The divergent $F_{2}$ feature leading to additional electrons in the range between $R_{1/2}^{H}$ and $R_{\rm c}$. The charge distribution in high-pressure with $F_{2}$ feature implies the local orbital lost the typical characteristic of 4$f$ like $F_{1}$.

The Fig.7 shows two radial parts of local orbital functions of Yb at low and high pressure respectively. These local orbitals are from the occupied KS-orbitals close to the Fermi energy at Gamma point, which are almost f-feature. It can be clearly seen that at low pressure (blue lines), the orbital shows a typical 4$f$-feature as $F_{1}$, and almost no $F_{2}$ feature expresses. For the orbital at high pressure, an $F_{2}$ feature shows in the range $R_{1/2}^{H}$ and $R_{\rm c}$, which is consistent with the charge distribution results. Another important feature of the high-pressure orbital is that a “node” appears just at the position $R_{1/2}^{H}$. It reveals a probable explanation for the occupation over 1.0. The overlap of the augmented sphere leads to the invasion of the nearest orbital, and makes the feature far from typical 4$f$. Then, the weight of $F_{2}$ becomes larger at the projection operator. The feature of 4$f$ near the nucleus still needs to be described with $F_{1}$. At last, the excessive $F_{2}$ component together with the $F_{1}$ leads to an over 1.0 occupation.

It is demonstrated that the overlap of local orbital in PAW can cause an abnormal local density matrix, and one can explain the abnormal behavior of simulation results by PAW from this point of view. For the E-V, P-V and B-V results, the LAPW results increase faster with the decreasing of the volume than PAW. It is the results of total energy correction with FLL double counting term and local density matrix difference between PAW and LAPW. There are two different conditions. One is the occupation of every local orbital is lower than 0.5 (as in La). In this condition, the total energy grows faster when the occupation of the local density matrix grows faster, and the P-V and B-V become steeper at a small volume. The other condition is that all local orbitals are over 0.5 occupations (as in Yb). At this condition, the total energy grows slower when the occupation of the local density matrix grows faster, and the P-V, B-V becomes less steep at a small volume. The case between these two conditions (as in Gd) leads to the cross of E-V, P-V and B-V curves between PAW and LAPW results. Since the DFT results of PAW and LAPW are similar at both low and high-pressure conditions, the contribution from the DFT part to the total energy in DFT+$U$ is similar. The main difference between the DFT+$U$ results of PAW and LAPW comes from the “+$U$” correction. The “+$U$” correction based on PAW in VASP has a larger projection range, thus more electrons are counted in the local density matrix. Together with the overlap of local orbital, the local orbital occupation increases faster in PAW than it in the LAPW with the decreasing of volume. It makes the “+$U$” correction shows a more obvious effect in PAW results, which may even lead to an unphysical result.

For the abnormal behavior of Yb in PAW results, that the B-V results show a decrease with the volume decreasing, it is caused by a similar mechanism as above. In the high press condition of Yb, the occupation of local orbitals is all over 1.0 in PAW. As having been introduced at the beginning of this part, the contribution of the “+$U$” correction of energy becomes negative when the occupation is over 1.0. With the volume decreasing and the occupation of 4$f$ increasing, the negative contribution to the total energy in Yb becomes remarkable. This negative contribution from the “+$U$” correction suppresses a normal trend described by the DFT part. The P-V should also show a similar behavior in the PAW, if the volume is smaller.

The difference of electron structure results between PAW and LAPW can be explained from two aspects. Firstly, the large local orbital range in PAW has more weight in the KS orbital than LAPW dose. And it makes “+$U$” correction in PAW has more influence on the band structure. Secondly, the large projection range in PAW makes more occupation on the local orbitals, and more occupation means lower orbital energy in “+$U$” correction. It leads to a lower band energy in occupied orbitals in PAW. Though the electron structure results of PAW and LAPW show some differences, no obvious qualitative difference as in the E-V, P-V and B-V results has been captured. It may be understood from the fact that the “+$U$” correction of orbital is linear, but the correction of total energy is quadratic.

Another phenomenon observed above is that the DFT predicts higher pressure than other correction methods in small volume (high-pressure) range for Tm, Yb and Lu. This underestimating of pressure appears in all correction methods we tested. Similar results are reported in many “$+U$” and hybrid functional works at ambient pressure, where the correction methods give a larger volume and a smaller bulk modulus comparing with the DFT resultsYu2019 ; Loschen2007 .

A usual interpretation for this phenomenon is that the correction methods make the correlation electrons more concentrated around the nuclei which leads the atoms behaving like isolated atoms. The isolated behavior causes the binding between the atoms are weakened, and induce a larger volume and a smaller bulk modulus. However, in high-pressure condition with hundreds of gigapascals, the isolated atom interpretation may failed. We propose a qualitative understanding for the phenomenon in high-pressure range from the point that the capability of localizing electrons in HSE06 and “+$U$” may increase the binding between atoms. Different from the condition at ambient pressure, when the volume is small enough, there is not much space for the $f$ electrons concentrated around the atoms. The electrons may prefer to localize between atoms, and the correction methods still can stabilize these local electrons. The concentration of electrons between atoms can binding the nucleus tight, which seems like some kind of “bonding effect”. This effect makes the correction methods to predict a much smaller volume in high-pressure condition, or smaller pressure with the same volume. An illustration try to explain the picture is shown in Fig.8 schematically. However, it is just a rough interpretation qualitatively, and more detailed investigation is required for clearly understanding the results.

In the frame of LAPW, the projection operator acts on the single electron orbital $\psi_{i}(\mathbf{r})$ to obtain local density matrix, which has the formation like,

$n_{mm^{\prime}}$ is the element of matrix with quantum number $m$ and $m^{\prime}$. The summation over all electron orbitals $i$, which are occupied. In LAPW methods, the radial part of the projection operator is determined by the radial part of the augmented wave function and is truncated by the radius of augmented range. Thus ${\hat{P}}_{mm^{\prime}}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$ only extract the information of subspace by the angular part. The ${\hat{P}}_{mm^{\prime}}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$ can be written as,

$\theta_{\alpha}\left(r\right)$ is a step function which is 1 within the augmented range $\alpha$, and is 0 outside. $\delta\left(r-r^{\prime}\right)$ is a delta function. $Y_{lm}(\hat{\mathbf{r}})$ is the spherical harmonic with angular momentum $l,\ m$. Since ${\hat{P}}_{mm^{\prime}}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$ do not specify the feature of the radial part of local orbital, it is important to show the form of radial wave function if LAPW orbitals. Here shows a basis with reciprocal coordinate $\mathbf{k}+\mathbf{G}$

$\alpha$ and $I$ stand for the augmented and interval ranges respectively. $\theta_{I}\left(\mathbf{r}\right)$ and $\theta_{\alpha}\left(\mathbf{r}\right)$ are the step function, which makes corresponding restricted in the certain range. In the interval range, the plane wave is used. In augmented spherical range an atom-like wave function is used. $u_{\alpha l}\left(r_{\alpha}\right)$ and ${\dot{u}}_{\alpha l}\left(r_{\alpha}\right)$ are the solution of radial Schrodinger function at certain energy with spherical averaged KS potential and its derivative function to energy respectively. The KS single electron orbital can be described by a set of these basis.

The summation is over reciprocal lattice $\mathbf{G}$. $i$ stands for the order of KS orbital. Thus, in the augmented range, the KS orbital is

with $\mathcal{A}_{\alpha lm}^{i}=\sum_{\mathbf{G}}{C_{\mathbf{G}}^{i}A_{\alpha lm}^{i,\mathbf{G}}}$ and $\mathcal{B}_{\alpha lm}^{i}=\sum_{\mathbf{G}}{C_{\mathbf{G}}^{i}B_{\alpha lm}^{i,\mathbf{G}}}$.

Since $\left\langle u_{i}\middle|u_{i}\right\rangle=1$ and $\left\langle u_{i}\middle|{\dot{u}}_{i}\right\rangle=0$ in LAPW method, the element of local density matrix is

PAW is another important augmented method. The formation of the local density matrix and projection operator is the same as the ones in LAPW methods. The difference of the element in the local density matrix between PAW and LAPW comes from the different radial orbital of these two methods. In PAW, the KS orbital has the form,

$\left|\phi_{s}\right\rangle$ and $\left|{\widetilde{\phi}}_{s}\right\rangle$ are also atom-like orbitals with the form,

Subscript $s$ is a collection of angular momentum $lm$, atom site $\alpha$ and the order of radial function $t$. With the definition, the element of local density matrix in PAW is

and