Nonlocal pseudopotential energy density functional for orbital-free density functional theory.

Orbital-free density functional theory (OF-DFT) is an electronic structure method with a low computational cost that scales linearly with the number of simulated atoms, making it suitable for large-scale material simulations. It is generally considered that OF-DFT strictly requires the use of local pseudopotentials, rather than orbital-dependent nonlocal pseudopotentials, for the calculation of electron-ion interaction energies, as no orbitals are available. This is unfortunate situation since the nonlocal pseudopotentials are known to give much better transferability and calculation accuracy than local ones. We report here the development of a theoretical scheme that allows the direct use of nonlocal pseudopotentials in OF-DFT. In this scheme, a nonlocal pseudopotential energy density functional is derived by the projection of nonlocal pseudopotential onto the non-interacting density matrix (instead of "orbitals") that can be approximated explicitly as a functional of electron density. Our development defies the belief that nonlocal pseudopotentials are not applicable to OF-DFT, leading to the creation for an alternate theoretical framework of OF-DFT that works superior to the traditional approach.

Introduction

Ab initio calculations using Kohn-Sham (KS) density functional theory (DFT)[1,2] can accurately describe the fundamental properties of various materials. However, its computational cost scales with the cube of the number of electrons in the simulation cell, which poses a major challenge to large-scale simulations. In contrast, orbital-free (OF) DFT is inherent of the lower computational cost that scales linearly with the number of atoms in the system, as it relies only on the electron density and the use of KS orbitals is avoided. As a result, OF-DFT is successfully applied to large-scale simulations of systems with up to millions of atoms[3-6].

The accuracy of OF-DFT simulations depends strongly on the quality of the non-interacting kinetic energy and the electron-ion (or electron-pseudocore) interaction energy employed in the simulations. Many approximate kinetic energy density functionals (KEDFs) have been proposed to evaluate the non-interacting kinetic energy in OF-DFT[7-39]. Their use in combination with local pseudopotentials[40-44] can achieve results that agree reasonably with those derived by KS-DFT, especially for main-group metals, III–V semiconductors[24,34,36,45], and even systems with inhomogeneous electron density such as metal clusters and quantum dots[37,39].

Unfortunately, the local pseudopotentials[28,31,32,40,41,43] used to evaluate the electron-ion interaction energy suffer from a lack of transferability[43], as they fail to reproduce the correct scattering behavior of the all-electron potentials[46-48]. Overcoming the transferability problem requires a reliance on either all-electron potential or nonlocal pseudopotentials (NLPPs), which are widely used in orbital-based approaches. However, it is practically unfeasible to use the all-electron potential, as an accurate all-electron KEDF for OF-DFT calculations is not yet available[49,50]. Furthermore, the use of NLPPs[46,51] runs against conventional understanding, as no orbitals are available in the traditional framework of OF-DFT[41,42,44,52,53].

A crucial nonlocal energy term with a set of angular-momentum-dependent energies has recently been added to OF-DFT calculations[54,55] in an effort to correct errors arising from the use of KEDFs and local pseudopotentials. This approach has successfully reproduced the bulk properties of several standard structures of Ti. However, special care must be taken when applying it to a wide range of practical simulations, as frozen on-site orbitals and empirically directed fitting parameters are part of the model[52]. There is substantial demand for a general approach to evaluate the electron-ion interaction energy using NLPPs in OF-DFT calculations. In this manuscript, we developed a theoretical scheme that allows the direct use of the NLPPs for the calculation of electron-ion interaction energy in OF-DFT, together with a specially designed theoretical framework of OF-DFT. This development leads to an OF-DFT calculation that gives a better transferability than the existing OF-DFT method based on local pseudopotentials.

Results and discussion

Nonlocal pseudopotential energy density functional

In general, the total energy density functional of OF-DFT can be expressed as:where ρ, T, E, E, E, {R}, and V are the electron density, KEDF, Hartree energy, exchange-correlation energy, ion-ion repulsion energy, the set of atomic positions, and local pseudopotential, respectively. To include the nonlocal electron-ion interactions, the total energy density functional of OF-DFT is reformulated as:where the total electron-ion interaction energy E[ρ] can be separated into two parts: a local part E[ρ] = ∫ρ(r)V(r)d3r and a nonlocal part E[ρ]. All of the terms in Eq. (2) except the nonlocal part of pseudopotential (E[ρ]) can be evaluated easily.

The exact NLPP energy depends on the KS orbitals or the density matrix:where f, , and represent the occupation number of the ith KS orbital ψ, the real-space representation of the nonlocal part pseudopotential, and the non-interacting density matrix, respectively. Considering that the density matrices can be used to approximate the KEDFs[7,56,57], an NLPP energy density functional (NLPPF) relying directly on the density matrix is proposed to evaluate the nonlocal electron-ion interaction energy. The nonlocal electron-ion interaction energy is then rewritten as a function of electron densityBy taking the Kleinman-Bylander form[58] of norm-conserving NLPPs, the nonlocal part pseudopotential[59] can be written aswhere and . The terms and are the atomic pseudo-wave-function and the short-range pseudopotential corresponding to the lmth angular momentum of ath atom, respectively. γ[ρ] denotes the density matrix as a function of electron distribution ρ. Although there is no exact analytic form available for the density matrix functional, a modified Gaussian (MG)[60] form derived from the second-order Taylor expansions of the density matrix[61] was employed to approximate the density matrix functional:where and . The second term in the square bracket is O(s4) correction[60], where A is an adjustable parameter. β(r) denotes the “local temperature” [62,63], where t(r) is the exact kinetic energy density defined as and ρ(r) = ∣ψ(r)∣2 is ith KS orbital’s density (see Refs. [60,61]). To remove the orbital-dependent problem in t(r) and obtain a solely density-dependent form of the density matrix functional, the kinetic energy density is obtained directly from the integrand of KEDFs to replace the exact one: t(r) ≈ t[ρ](r). The widely used Wang-Teter (WT) KEDF[28] is chosen as an exemplary case, and t[ρ](r) can be expressed aswhere is the kernel of WT functional. The Supplementary Notes give the details of the kinetic energy densities obtained from WT and Xu-Wang-Ma[38] KEDFs.

The direct numerical evaluations of and at the average position are very complicated. They are therefore approximated using q-mean “nonlocal density” and two-point average temperature for systems with slowly varying electron densities. The density matrix functional of Eq. (6) can then be reformulated asCombining Eqs. (4), (5), and (8) gives the NLPPF aswhere the integral domain Ω is the ath ionic core region. Owing to the short-range nature of , the computational cost of Eq. (9) scales linearly with the number of atoms (N), where c can be regarded as a constant derived from the double integral within the near-core region. Within the NLPPF scheme, a new theoretical framework of OF-DFT has been built and implemented in ATLAS[5,64]. The further computational details are provided in Methods Section. The parameters of pseudopotential and NLPPF are presented in Supplementary Tables 1 and 2, respectively.

Computational accuracy of NLPPF scheme

We first applied this scheme for OF-DFT calculations of Li, Mg, and Cs within hexagonal-close-packed (HCP), face-centered cubic (FCC), body-centered cubic (BCC), simple cubic (SC), and cubic diamond structures. For each structure, 13 energy-volume points were calculated by expanding and compressing the approximate equilibrium volume by up to 20%, and the bulk properties (the equilibrium cell volume V0, bulk modulus B0, and the relative energy E with respect to HCP structure) were determined by fitting the energy-volume curve against Murnaghan’s equation of state[65]. The comparison of the results obtained by OF-DFT using both local pseudopotentials (BLPS[41] and OEPP[43]) and NLPPs against those calculated by KS-DFT using the projector augmented-wave (PAW)[66] method are presented in Supplementary Table 3. For Li and Mg solids, the OF-DFT calculations within both local pseudopotentials and NLPPs give reasonable predictions of V0, B0, and E, which are comparable with the KS-DFT results. It is noteworthy that our scheme shows an improvement over the local pseudopotentials of OEPP for bulk Cs. The accurate bulk properties of Li/Mg/Cs obtained by the current scheme demonstrate its valid applicability to simple metallic solids.

Further assessment of the accuracy of our scheme was demonstrated by molecular dynamics calculations for Li-Mg alloy. The calculations used a canonical ensemble (at 1000 K) in a supercell containing 108 atoms (Li54Mg54). The calculated pair distribution functions g(r) for Li-Mg alloy are shown in Fig. 1. Overall, the predicted shapes and peak neighbors of pair distribution functions by OF-DFT within NLPPF match the results calculated by KS-DFT. Especially notable are the resulting contributions of the partial distributions (Fig. 1b–d) calculated by OF-DFT within the NLPPF being almost identical to the KS-DFT calculations, which are superior to those obtained by the local pseudopotentials.

Fig. 1

Pair distribution functions for Li54Mg54 alloy.

a Total, b Mg-Mg, c Li-Mg, and d Li-Li pair distribution functions.

Transferability of NLPPF scheme

To demonstrate the transferability of our scheme, we randomly generated 50 structures of Li systems using CALYPSO[67,68]. The total energies of these structures were calculated by OF-DFT and KS-DFT. The comparisons of total energy relative to the HCP structure (E) are shown in Fig. 2. The orderings of energy are well captured by OF-DFT within NLPPF. The relative energies of different phases are overall well reproduced and in reasonable agreement with the KS-DFT results. The least-square fitting lines of WT-NLPPF are generally closer to the KS-PAW results than those from the local pseudopotentials. For example, the mean error of E for Li systems obtained by OF-DFT within the WT-NLPPF is 45 meV/atom, which is lower than that within BLPS (73 meV/atom), or OEPP (201 meV/atom). Therefore, this framework of OF-DFT with improved transferability is superior to the traditional one.

Fig. 2

Relative energies for random structures of elemental Li.

The results are calculated by OF-DFT using BLPS, OEPP and NLPPF in comparison with that by KS-DFT using the PAW method. Blue dash-dotted and red dashed lines are the least-square fittings of WT-BLPS and WT-NLPPF results, respectively.

Previous studies have shown that OF-DFT with local pseudopotentials can be applied to most s- and p-block metals. However, OF-DFT simulation using the local pseudopotential OEPP shows unacceptable errors for various crystalline phases of Be (Fig. 3): the curves of energy with respect to volume for HCP and FCC structures show the total energy monotonically decreasing with increasing volume. In contrast, the curves with clear minima predicted by OF-DFT within NLPPF agree well with those produced by KS-DFT. These findings indicate the significant superiority of our proposed framework over the conventional one.

Fig. 3

Relative energy versus volume curves for Be systems.

a The calculated energy-volume curves of Be-HCP. b The calculated energy-volume curves of Be-FCC. The total energy shift of WT-OEPP is −33.030 eV/atom.

Note that the bulk properties (e.g., equilibrium volume, bulk modulus, and relative energy) calculated by OF-DFT within NLPPF reproduce the results of KS-DFT for Li and Mg almost exactly (Table 1), considering the maximal deviation of E is within 32 meV/atom. However, there are some discrepancies for crystalline phases of Be: in particular, the deviation of E for the SC structure is larger than 400 meV/atom. To explore the causes of these discrepancies, we estimated the errors of the kinetic energy density of the WT-KEDF with respect to the KS kinetic energy density along the [100] and [111] directions in the SC structures of Li, Mg, and Be (Fig. 4). The kinetic energy density of KS-DFT is clearly reproduced accurately by the WT-KEDF for Li and Mg with slow variations of electron densities. However, it is seriously underestimated for Be, in which the electron distribution rapidly varies in the near-core region. Therefore, we believe that errors in the kinetic energy densities for Be lead to the discrepancy in its bulk properties obtained by the framework of OF-DFT within the NLPPF. The findings are fairly consistent with our expectation that the performance of the NLPPF relies strongly on the accuracy of kinetic energy density, as manifested by Eqs. (7) and (8).

Table 1

B0 (GPa), E (eV/atom), and V0 (Å3/atom) for bulk Li, Mg, and Be by KS-DFT and OF-DFT.

		Method	HCP	FCC	BCC	SC
Li	B0	KS-PAW	13.9	13.6	13.9	12.1
		WT-NLPPF	13.5	13.5	13.7	11.0
	V0	KS-PAW	20.280	20.372	20.396	20.580
		WT-NLPPF	19.483	19.462	19.352	20.844
	ER	KS-PAW	0.000	0.000	0.001	0.120
		WT-NLPPF	0.000	0.000	0.001	0.152
Mg	B0	KS-PAW	35.8	35.5	34.8	22.7
		WT-NLPPF	33.0	31.3	31.3	21.2
	V0	KS-PAW	22.838	23.071	22.826	27.478
		WT-NLPPF	23.194	23.924	23.730	28.274
	ER	KS-PAW	0.000	0.012	0.029	0.382
		WT-NLPPF	0.000	0.011	0.031	0.372
Be	B0	KS-PAW	123.3	119.7	124.1	74.5
		WT-NLPPF	91.5	90.5	87.2	63.3
	V0	KS-PAW	7.910	7.875	7.822	10.274
		WT-NLPPF	7.690	7.942	7.798	10.160
	ER	KS-PAW	0.000	0.080	0.099	1.004
		WT-NLPPF	0.000	0.058	0.082	0.561

Fig. 4

Kinetic energy differences in SC structure for Li, Mg, and Be systems.

a The calculated differences along [100] direction. b The calculated differences along [111] direction. Note that .

Due to the existence of significant differences in kinetic energy density between WT and the exact one for Cd systems including localized d-channel electrons (see Supplementary Fig. 1), the NLPPF using WT cannot be applicable to investigate Cd systems with d-channel electrons. Therefore, the NLPP of Cd was constructed without d-channel electrons for the additional calculations. As listed in Table 2, the bulk properties predicted by the OF-DFT within NLPPF framework agree fairly well with the predictions by KS-DFT using the same NLPP, whereas WT-OEPP gives serious discrepancies compared with KS-NLPP results in all bulk properties. Although the calculations using NLPP without d-channel electrons cannot give the accurate bulk properties for Cd systems (Table 2), OF-DFT within NLPPF works superior to that within OEPP. Overall, it can be expected this NLPPF scheme using accurate KEDF and its kinetic energy density can be applied to the systems including localized electrons, such as the transition metals or covalent systems.

Table 2

B0 (GPa), E (eV/atom), and V0 (Å3/atom) for bulk Cd by KS-DFT and OF-DFT.

		Method	HCP	FCC	BCC	SC
Cd	B0	KS-PAW	40.8	41.0	34.5	29.2
		KS-NLPP	67.4	66.2	66.1	39.8
		WT-OEPP	148.2	138.8	147.1	74.8
		WT-NLPPF	67.0	64.8	65.2	39.8
	V0	KS-PAW	22.758	22.979	23.512	27.141
		KS-NLPP	15.684	15.851	15.674	18.487
		WT-OEPP	10.379	10.693	10.379	11.964
		WT-NLPPF	15.702	16.027	15.763	18.878
	ER	KS-PAW	0.000	0.004	0.053	0.121
		KS-NLPP	0.000	0.019	0.036	0.436
		WT-OEPP	0.000	0.138	0.072	0.979
		WT-NLPPF	0.000	0.036	0.056	0.451

Computational efficiency of NLPPF scheme

To assess the computational efficiency of the current scheme, static simulations of Cs BCC supercells containing 128–16,000 atoms were performed by OF-DFT within NLPPF. The total wall time of the single point energy calculations is plotted with respect to the number of atoms in Supplementary Fig. 2. The computational cost of this framework clearly scales linearly with the number of atoms in the simulation cell, in sharp contrast to the cubic scaling of KS-DFT. This shows that the OF-DFT within NLPPF is potentially applicable to the simulation of large-scale systems containing millions of atoms.

In summary, we proposed an NLPPF scheme that allows the direct use of NLPPs in OF-DFT calculations. The static and dynamic properties of s- and p-block metals calculated within this scheme agree well with KS-DFT predictions and show significant improvements in the computational accuracy and transferability over conventional OF-DFT with local pseudopotentials. With this work, we defy the conventional wisdom of orbital-dependent NLPPs being incompatible with OF-DFT, leading to the creation of an alternative framework of OF-DFT, which opens up new avenues for further development of the theory.

Methods

Pseudopotential generations

The Troullier-Martins NLPPs[59] are generated by the FHI98PP[47] code for all considered systems [see Supplementary Table 1] and the p-channel of the NLPPs is used as the local pseudopotential of V(r) in OF-DFT.

Numerical calculations

The KS-DFT calculations using the PAW[66] and NLPP are performed by VASP[69,70] and ARES packages[71], respectively. The k-point meshes are generated using the Monkhorst-Pack method[72] with the k-spacing of 0.10 Å−1. The kinetic energy cutoff is 500 eV for all the simulations using VASP. The OF-DFT calculations are carried out by ATLAS[5,64] using WT[28] as KEDF, and the corresponding kinetic energy density is used to construct the NLPPF. The generalized gradient approximation with the form of Perdew-Burke-Ernzerhof[73] is employed for both OF-DFT and KS-DFT calculations. The grid spacings of 0.18, 0.18, 0.22, 0.10, 0.15, 0.22, 0.12, and 0.15 Å are used in ATLAS/ARES for Li, Mg, Cs, Be, Cd, K, Zn, and Li-Mg alloy, respectively. The parameters of A and q in NLPPFs are presented in Supplementary Table 2 carefully tuned to yield the bulk properties, which agree with the KS-DFT (NLPPs) predictions.

Molecular dynamics

The molecular dynamic simulations of Li-Mg alloy are performed in the canonical ensemble (at 1000 K) applying the Nosé-Hoover thermostat[74,75] simulations up to 10 ps (0.5 fs/step), with the first 10,000 steps for equilibrating the system. The data for further analysis were collected from the subsequent 10,000 steps.