DFT-FE -- A massively parallel adaptive finite-element code for large-scale density functional theory calculations

TL;DR

DFT-FE is a highly scalable, accurate finite-element code for large-scale density functional theory calculations, outperforming traditional plane-wave methods in speed and efficiency for systems with thousands of electrons.

Contribution

This work introduces DFT-FE, a massively parallel finite-element DFT code that handles large systems efficiently with adaptive basis and reduced computational costs.

Findings

Achieves over 5-10x speedup compared to plane-wave codes for systems with >10,000 electrons.

Demonstrates strong parallel scalability up to 192,000 MPI tasks.

Maintains high accuracy in energies, forces, and stresses across diverse systems.

Abstract

We present an accurate, efficient and massively parallel finite-element code, DFT-FE, for large-scale ab-initio calculations (reaching $\sim 100,000$ electrons) using Kohn-Sham density functional theory (DFT). DFT-FE is based on a local real-space variational formulation of the Kohn-Sham DFT energy functional that is discretized using a higher-order adaptive spectral finite-element (FE) basis, and treats pseudopotential and all-electron calculations in the same framework, while accommodating non-periodic, semi-periodic and periodic boundary conditions. We discuss the main aspects of the code, which include, the various strategies of adaptive FE basis generation, and the different approaches employed in the numerical implementation of the solution of the discrete Kohn-Sham problem that are focused on significantly reducing the floating point operations, communication costs and latency. We demonstrate the accuracy of DFT-FE by comparing the energies, ionic forces and periodic cell stresses on a wide range of problems with popularly used DFT codes. Further, we demonstrate that DFT-FE significantly outperforms widely used plane-wave codes---both in CPU-times and wall-times, and on both non-periodic and periodic systems---at systems sizes beyond a few thousand electrons, with over $5-10$ fold speedups in systems with more than 10,000 electrons. The benchmark studies also highlight the excellent parallel scalability of DFT-FE, with strong scaling demonstrated on up to 192,000 MPI tasks.

Tables

Table 1. Table 1: Comparison between user defined adaptive mesh refinement strategy (UDAMR) and automatic adaptive mesh generation strategy (AAMR). Case Study: Pseudopotential non-periodic DFT calculation on SiF 4 molecule. SiF 4 details: Si-F bondlength = 2.91 2.91 2.91 Bohr, F-Si-F Tetrahedral angle = 109.47 0 superscript 109.47 0 {109.47}^{0} , E 0 = − 19.94369693 subscript 𝐸 0 19.94369693 E_{0}=-19.94369693 Ha/atom, 𝒇 0 subscript 𝒇 0 \bm{f}_{0} on F atom = ( 0.0288669 , 0.0 , − 0.0204124 ) 0.0288669 0.0 0.0204124 (0.0288669,0.0,-0.0204124) Ha/Bohr. β = 0.03 𝛽 0.03 \beta=0.03 in Algorithm   1 . | Δ ​ E g | = | E 0 − E h | Δ subscript 𝐸 𝑔 subscript 𝐸 0 subscript 𝐸 ℎ |\Delta E_{g}|=|E_{0}-E_{h}| and | Δ ​ f | = ‖ f 0 − f h ‖ 2 Δ 𝑓 subscript norm subscript f 0 subscript f ℎ 2 |\Delta f|=||\textbf{f}_{0}-\textbf{f}_{h}||_{2} are the discretization errors in ground-state energy and force vector, respectively.

Adaptive mesh	$h_{\text{min}}$, $h_{\text{max}}$ (a.u.), FEord ;	FE basis per atom	$\left|\Delta E_g\right|$	$\left|\Delta f\right|$
strategy	(tolq (Ha/atom))
UDAMR	0.72, 12.5, 4	24,965	1.5 $\times$ 10-3	1.5 $\times$ 10-4
UDAMR	0.36, 12.5, 4	193,108	1.1 $\times$ 10-5	1.72 $\times$ 10-5
AAMR	0.28, 20, 4 ; (1 $\times$ 10-3)	26,115	3.6 $\times$ 10-4	4.6 $\times$ 10-4
AAMR	0.28, 20, 4 ; (5 $\times$ 10-5)	47,495	6 $\times$ 10-5	4.1 $\times$ 10-5
AAMR	0.15, 20, 4 ; (2 $\times$ 10-5)	58,520	3 $\times$ 10-5	2.6 $\times$ 10-5

Table 2. Table 2: Comparison between user defined adaptive mesh refinement strategy (UDAMR) and automatic adaptive mesh generation strategy (AAMR). Case Study: All-electron non-periodic calculation on SiF 4 molecule . Details: Si-F bondlength = 2.91 Bohr, F-Si-F Tetrahedral angle = 109.47 0 , E 0 = − 137.76973324 subscript 𝐸 0 137.76973324 E_{0}=-137.76973324 Ha/atom, f 0 subscript 𝑓 0 f_{0} on 𝒇 0 subscript 𝒇 0 \bm{f}_{0} on F atom = ( 0.02595397 , 0.0 , − 0.01834901 ) absent 0.02595397 0.0 0.01834901 =(0.02595397,0.0,-0.01834901) Ha/Bohr. β = 0.03 𝛽 0.03 \beta=0.03 in Algorithm   1 . | Δ ​ E g | = | E 0 − E h | Δ subscript 𝐸 𝑔 subscript 𝐸 0 subscript 𝐸 ℎ |\Delta E_{g}|=|E_{0}-E_{h}| and | Δ ​ f | = ‖ f 0 − f h ‖ 2 Δ 𝑓 subscript norm subscript f 0 subscript f ℎ 2 |\Delta f|=||\textbf{f}_{0}-\textbf{f}_{h}||_{2} are the discretization errors in ground-state energy and force vector, respectively.

Adaptive mesh	$h_{\text{min}}$, $h_{\text{max}}$ (a.u.), FEord ;	FE basis per atom	$\left|\Delta E_g\right|$	$\left|\Delta f\right|$
strategy	(tolq (Ha/atom))
UDAMR	0.024, 12.5, 4	382,705	2 $\times$ 10-3	1.8 $\times$ 10-3
UDAMR	0.012, 12.5, 4	2,493,531	3.9 $\times$ 10-4	6.8 $\times$ 10-4
AAMR	0.019, 20, 4 ; (1 $\times$ 10-2)	91,002	1.02 $\times$ 10-3	3.3 $\times$ 10-4
AAMR	0.0097, 20, 4 ; (1 $\times$ 10-3)	125,083	1.6 $\times$ 10-4	1.2 $\times$ 10-4
AAMR	0.0024, 20, 4 ; (5 $\times$ 10-5)	187,591	9.8 $\times$ 10-6	2.15 $\times$ 10-5

Table 3. Table 3: Comparison between user defined adaptive mesh refinement strategy (UDAMR) and automatic adaptive mesh generation strategy (AAMR). Case Study: All-electron periodic calculation on Si diamond unit-cell . Si unit-cell details: Lattice constant = 10.0065 ​ B ​ o ​ h ​ r absent 10.0065 𝐵 𝑜 ℎ 𝑟 =10.0065~{}Bohr , E 0 = − 289.350382455 subscript 𝐸 0 289.350382455 E_{0}=-289.350382455 Ha/atom, σ 0 subscript 𝜎 0 \sigma_{0} on unit-cell = − 0.00077999 absent 0.00077999 =-0.00077999 Ha/Bohr 3 . β = 0.03 𝛽 0.03 \beta=0.03 in Algorithm   1 . | Δ ​ E g | = | E 0 − E h | Δ subscript 𝐸 𝑔 subscript 𝐸 0 subscript 𝐸 ℎ |\Delta E_{g}|=|E_{0}-E_{h}| and | Δ ​ σ | = | σ 0 − σ h | Δ 𝜎 subscript 𝜎 0 subscript 𝜎 ℎ |\Delta\sigma|=|\sigma_{0}-\sigma_{h}| are the discretization errors in ground-state energy and hydrostatic stress, respectively.

Adaptive mesh	$h_{\text{min}}$, $h_{\text{max}}$ (a.u.), FEord;	FE basis per atom	$\left|\Delta E_g\right|$	$\left|\Delta \sigma \right|$
strategy	(tolq (Ha/atom))
UDAMR	0.013, 1.67, 4	188,724	2.0 $\times$ 10-3	6.9 $\times$ 10-5
UDAMR	0.013, 1.67, 5	360,695	3.6 $\times$ 10-4	4.3 $\times$ 10-5
AAMR	0.0097, 1.25, 4 ; (1 $\times$ 10-2)	75,563	1.75 $\times$ 10-3	5.48 $\times$ 10-5
AAMR	0.0048, 1.25, 4 ; (1 $\times$ 10-3)	111,860	8.0 $\times$ 10-5	2.67 $\times$ 10-5
AAMR	0.0024, 1.25, 4 ; (5 $\times$ 10-5)	283,238	6.7 $\times$ 10-6	1.5 $\times$ 10-5

Table 4. Table 4: Accuracy and robustness study of mixed precision computations in CholGS orthonormalization and Rayleigh-Ritz procedure on benchmark systems. Energy difference, maximum atomic force difference magnitude ( max 1 ≤ i ≤ N a ⁡ ‖ 𝐟 dp i − 𝐟 sp i ‖ subscript 1 𝑖 subscript 𝑁 𝑎 norm subscript superscript 𝐟 𝑖 dp subscript superscript 𝐟 𝑖 sp \max\limits_{1\leq i\leq N_{a}}\left\|{\bf f}^{i}_{\rm dp}-{\bf f}^{i}_{\rm sp}\right\| ) and total number of SCFs are reported with respect to double precision calculations. 𝐟 dp i subscript superscript 𝐟 𝑖 dp {\bf f}^{i}_{\rm dp} and 𝐟 sp i subscript superscript 𝐟 𝑖 sp {\bf f}^{i}_{\rm sp} denote atomic force on i th superscript 𝑖 th {i}^{\rm th} atom for double precision and single precision calculations respectively. Discretization errors for the benchmark systems are ∼ 10 − 4 similar-to absent superscript 10 4 \sim 10^{-4} Ha/atom in ground-state energy and ∼ 10 − 4 similar-to absent superscript 10 4 \sim 10^{-4} Ha/Bohr in ionic forces. More details about the benchmark systems are given in Section  4.3 .

System	Energy difference	Maximum force difference	Total SCFs
	(Ha/atom)	magnitude (Ha/Bohr)	(Double, Mixed)
${\text{Mg}}_{\text{6x6x6}}$	$7\times {10}^{-12}$	$2\times {10}^{-6}$	(49, 49)
${\text{Cu}}_{\text{4-shell}}$	$5\times {10}^{-12}$	$3\times {10}^{-6}$	(46, 46)
${\text{Mo}}_{\text{6x6x6}}$	$3\times {10}^{-12}$	$7\times {10}^{-7}$	(49, 49)

Table 5. Table 5: Validation of DFT-FE with QE on pseudopotential benchmark systems at two different accuracy levels. FE ord , h min subscript FE ord subscript ℎ min {\rm FE}_{\rm ord},h_{\rm min} and h max subscript ℎ max h_{\rm max} denote the FE polynomial order, minimum element size and maximum element size (Bohr), respectively, in DFT-FE. E cut subscript 𝐸 cut {E}_{\rm cut} denotes the plane-wave basis cut-off used in QE (Hartree). E g subscript 𝐸 g {E}_{\rm g} denotes ground-state energy (Hartree/atom). Δ max ​ f = max 1 ≤ i ≤ N a ⁡ ‖ 𝐟 i DFT − FE − 𝐟 i QE ‖ subscript Δ max 𝑓 subscript 1 𝑖 subscript 𝑁 𝑎 norm superscript subscript 𝐟 𝑖 DFT FE superscript subscript 𝐟 𝑖 QE \Delta_{\rm max}{f}=\max\limits_{1\leq i\leq N_{a}}\left\|{\bf f}_{i}^{\rm DFT-FE}-{\bf f}_{i}^{\rm QE}\right\| (Hartree/Bohr), where 𝐟 i subscript 𝐟 𝑖 {\bf f}_{i} denotes the force on the i th superscript 𝑖 th i^{\rm th} atom. Δ ​ σ h = | σ h DFT − FE − σ h QE | Δ subscript 𝜎 ℎ superscript subscript 𝜎 ℎ DFT FE superscript subscript 𝜎 ℎ QE \Delta{\sigma_{h}}=\absolutevalue{\sigma_{h}^{\rm DFT-FE}-\sigma_{h}^{\rm QE}} (Hartree/ Bohr 3 superscript Bohr 3 {\rm Bohr}^{3} ), where σ h subscript 𝜎 ℎ \sigma_{h} denotes the hydrostatic cell stress.

System	DFT-FE	DFT-FE	QE	QE	Difference in forces
	$({\mathrm{FE}}_{\mathrm{ord}},h_{\min },h_{\max })$	$E_g$	$E_{\mathrm{cut}}$	$E_g$	& stress $({\Delta }_{\max }f,\Delta {\sigma }_h)$
${\text{Mg}}_{\text{2x2x2}}$	$4$, $0.46$, $1.92$	$-54.3195364$	$45$	$-54.3195594$	${\Delta }_{\max }f=2.1\times {10}^{-4}$
					$\Delta {\sigma }_h=3.7\times {10}^{-6}$
${\text{Mg}}_{\text{4x4x4}}$	$4$, $0.46$, $1.92$	$-54.3279442$	$45$	$-54.3279638$	${\Delta }_{\max }f=3.3\times {10}^{-4}$
					$\Delta {\sigma }_h=4.6\times {10}^{-6}$
${\text{Mo}}_{\text{2x2x2}}$	$5$, $0.74$, $1.49$	$-68.5573334$	$20$	$-68.5573613$	${\Delta }_{\max }f=1\times {10}^{-5}$
					$\Delta {\sigma }_h=1.6\times {10}^{-6}$
${\text{Mo}}_{\text{4x4x4}}$	$5$, $0.74$, $1.49$	$-68.5811483$	$20$	$-68.5811857$	${\Delta }_{\max }f=2.1\times {10}^{-5}$
					$\Delta {\sigma }_h=1.9\times {10}^{-6}$
${\text{Cu}}_{\text{3-Shell}}$	$4$, $0.39$, $12.5$	$-182.5870759$	$50$	$-182.5870221$	${\Delta }_{\max }f=7.2\times {10}^{-5}$
${\text{Cu}}_{\text{4-Shell}}$	$4$, $0.39$, $12.2$	$-182.5908621$	$50$	$-182.5908346$	${\Delta }_{\max }f=1.4\times {10}^{-4}$

Table 6. (a) Medium accuracy level comparisons.

System	DFT-FE	DFT-FE	QE	QE	Difference in forces
	$({\mathrm{FE}}_{\mathrm{ord}},h_{\min },h_{\max })$	$E_g$	$E_{\mathrm{cut}}$	$E_g$	& stress $({\Delta }_{\max }f,\Delta {\sigma }_h)$
${\text{Mg}}_{\text{2x2x2}}$	$4$, $0.46$, $1.92$	$-54.3195364$	$45$	$-54.3195594$	${\Delta }_{\max }f=2.1\times {10}^{-4}$
					$\Delta {\sigma }_h=3.7\times {10}^{-6}$
${\text{Mg}}_{\text{4x4x4}}$	$4$, $0.46$, $1.92$	$-54.3279442$	$45$	$-54.3279638$	${\Delta }_{\max }f=3.3\times {10}^{-4}$
					$\Delta {\sigma }_h=4.6\times {10}^{-6}$
${\text{Mo}}_{\text{2x2x2}}$	$5$, $0.74$, $1.49$	$-68.5573334$	$20$	$-68.5573613$	${\Delta }_{\max }f=1\times {10}^{-5}$
					$\Delta {\sigma }_h=1.6\times {10}^{-6}$
${\text{Mo}}_{\text{4x4x4}}$	$5$, $0.74$, $1.49$	$-68.5811483$	$20$	$-68.5811857$	${\Delta }_{\max }f=2.1\times {10}^{-5}$
					$\Delta {\sigma }_h=1.9\times {10}^{-6}$
${\text{Cu}}_{\text{3-Shell}}$	$4$, $0.39$, $12.5$	$-182.5870759$	$50$	$-182.5870221$	${\Delta }_{\max }f=7.2\times {10}^{-5}$
${\text{Cu}}_{\text{4-Shell}}$	$4$, $0.39$, $12.2$	$-182.5908621$	$50$	$-182.5908346$	${\Delta }_{\max }f=1.4\times {10}^{-4}$

Table 7. (b) High accuracy level comparisons.

System	DFT-FE	DFT-FE	QE	QE	Difference in forces
	$({\mathrm{FE}}_{\mathrm{ord}},h_{\min },h_{\max })$	$E_g$	$E_{\mathrm{cut}}$	$E_g$	& stress $({\Delta }_{\max }f,\Delta {\sigma }_h)$
${\text{Mg}}_{\text{2x2x2}}$	5, 0.24, 0.96	$-54.3196337$	55	$-54.3196270$	${\Delta }_{\max }f=4.6\times {10}^{-6}$
					$\Delta {\sigma }_h=2.7\times {10}^{-7}$
${\text{Mg}}_{\text{4x4x4}}$	5, 0.24, 0.96	$-54.3280448$	55	$-54.3280318$	${\Delta }_{\max }f=7.9\times {10}^{-6}$
					$\Delta {\sigma }_h=2.6\times {10}^{-7}$
${\text{Mo}}_{\text{2x2x2}}$	5, 0.37, 0.74	$-68.5574282$	50	$-68.5574315$	${\Delta }_{\max }f=4.5\times {10}^{-6}$
					$\Delta {\sigma }_h=1.1\times {10}^{-6}$
${\text{Mo}}_{\text{4x4x4}}$	5, 0.37, 0.74	$-68.5812495$	50	$-68.5812527$	${\Delta }_{\max }f=6\times {10}^{-6}$
					$\Delta {\sigma }_h=9.8\times {10}^{-7}$
${\text{Cu}}_{\text{3-Shell}}$	5, 0.18, 12.5	$-182.5872871$	70	$-182.5872868$	${\Delta }_{\max }f=2.5\times {10}^{-5}$
${\text{Cu}}_{\text{4-Shell}}$	5, 0.18, 12.2	$-182.5910308$	70	$-182.5910298$	${\Delta }_{\max }f=3.5\times {10}^{-5}$

Table 8. Table 6: Validation of DFT-FE with exciting— and NWChem on all-electron benchmark systems. FE ord , h min subscript FE ord subscript ℎ min {\rm FE}_{\rm ord},h_{\rm min} and h max subscript ℎ max h_{\rm max} denote the FE polynomial order, minimum element size and maximum element size (Bohr), respectively, in DFT-FE. The rgkmax— parameter in exciting— determines the number of LAPW basis functions. The polarization consistent Gaussian basis sets in NWChem with increasing levels of accuracy are denoted by: pc2, pc3. E g subscript 𝐸 g {E}_{\rm g} denotes ground-state energy (Hartree/atom). Δ max ​ f = max 1 ≤ i ≤ N a ⁡ ‖ 𝐟 i DFT − FE − 𝐟 i QE ‖ subscript Δ max 𝑓 subscript 1 𝑖 subscript 𝑁 𝑎 norm superscript subscript 𝐟 𝑖 DFT FE superscript subscript 𝐟 𝑖 QE \Delta_{\rm max}{f}=\max\limits_{1\leq i\leq N_{a}}\left\|{\bf f}_{i}^{\rm DFT-FE}-{\bf f}_{i}^{\rm QE}\right\| (Hartree/Bohr), where 𝐟 i subscript 𝐟 𝑖 {\bf f}_{i} denotes the force on the i th superscript 𝑖 th i^{\rm th} atom. Δ ​ σ h = | σ h DFT − FE − σ h QE | Δ subscript 𝜎 ℎ superscript subscript 𝜎 ℎ DFT FE superscript subscript 𝜎 ℎ QE \Delta{\sigma_{h}}=\absolutevalue{\sigma_{h}^{\rm DFT-FE}-\sigma_{h}^{\rm QE}} (Hartree/ Bohr 3 superscript Bohr 3 {\rm Bohr}^{3} ), where σ h subscript 𝜎 ℎ \sigma_{h} denotes the hydrostatic cell stress.

System	DFT-FE	DFT-FE	exciting	exciting	Difference in forces
	$({\mathrm{FE}}_{\mathrm{ord}},h_{\min },h_{\max })$	$E_g$	rgkmax	$E_g$	& stress $({\Delta }_{\max }f,\Delta {\sigma }_h)$
${\text{Si}}_{\mathrm{uc}}^{\mathrm{Gamma}}$	5, 0.013, 1.67	$-289.3500196$	7	$-289.3501763$	$\Delta {\sigma }_h=\mathrm{\hspace{0.17em}6.1}\times {10}^{-6}$
${\text{Si}}_{\mathrm{uc}}^{\mathrm{Multi}-\mathrm{kpt}}$	5, 0.013, 1.67	$-289.4001025$	7	$-289.4001864$	$\Delta {\sigma }_h=\mathrm{\hspace{0.17em}4.0}\times {10}^{-6}$
${\text{SiVac}}_{2\times 2\times 2}$	5, 0.013, 1.67	$-289.3943335$	7	$-289.3944222$	${\Delta }_{\max }f=7.2\times {10}^{-4}$
					$\Delta {\sigma }_h=2.5\times {10}^{-6}$

Table 9. (a) Validation of periodic calculations with exciting—.

System	DFT-FE	DFT-FE	exciting	exciting	Difference in forces
	$({\mathrm{FE}}_{\mathrm{ord}},h_{\min },h_{\max })$	$E_g$	rgkmax	$E_g$	& stress $({\Delta }_{\max }f,\Delta {\sigma }_h)$
${\text{Si}}_{\mathrm{uc}}^{\mathrm{Gamma}}$	5, 0.013, 1.67	$-289.3500196$	7	$-289.3501763$	$\Delta {\sigma }_h=\mathrm{\hspace{0.17em}6.1}\times {10}^{-6}$
${\text{Si}}_{\mathrm{uc}}^{\mathrm{Multi}-\mathrm{kpt}}$	5, 0.013, 1.67	$-289.4001025$	7	$-289.4001864$	$\Delta {\sigma }_h=\mathrm{\hspace{0.17em}4.0}\times {10}^{-6}$
${\text{SiVac}}_{2\times 2\times 2}$	5, 0.013, 1.67	$-289.3943335$	7	$-289.3944222$	${\Delta }_{\max }f=7.2\times {10}^{-4}$
					$\Delta {\sigma }_h=2.5\times {10}^{-6}$

Table 10. (b) Validation of non-periodic calculations with NWChem.

System	DFT-FE (AAMR)	NWChem
	$\left({\mathrm{FE}}_{\mathrm{ord}},h_{\min },h_{\max }:E_g\right)$	(Gaussian Basis: $E_g$)
${\text{C}}_{\text{60}}$	4, 0.037, 20.0: $-38.0701278$	pc2: $-38.0693824$
	4, 0.0185, 20.0: $-38.0709949$	pc3: $-38.0710146$
	4, 0.0094, 20.0: $-38.0710439$	-

Table 11. Table 7: CPU-time comparison of DFT-FE with QE and ABINIT: Average time per SCF iteration step in Node-Hrs.

System	Number of atoms	FE basis	DFT-FE	Plane-wave	QE	ABINIT
	(Number of electrons)			basis
${\text{Mg}}_{\text{4x4x4}}$	255 (2,550)	6,673,513	0.3	530,051	0.1	0.3
${\text{Mg}}_{\text{6x6x6}}$	863 (8,630)	19,852,441	3.3	1,788,771	4.4	20.2
${\text{Mg}}_{\text{8x8x8}}$	2,047 (20,470)	45,954,505	21.6	4,240,071	123.5	-
${\text{Mg}}_{\text{10x10x10}}$	3,999 (39,990)	93,972,153	103.4	-	-	-

Table 12. (a) Benchmark system (i): hcp Mg periodic supercells with a mono-vacancy.

System	Number of atoms	FE basis	DFT-FE	Plane-wave	QE	ABINIT
	(Number of electrons)			basis
${\text{Mg}}_{\text{4x4x4}}$	255 (2,550)	6,673,513	0.3	530,051	0.1	0.3
${\text{Mg}}_{\text{6x6x6}}$	863 (8,630)	19,852,441	3.3	1,788,771	4.4	20.2
${\text{Mg}}_{\text{8x8x8}}$	2,047 (20,470)	45,954,505	21.6	4,240,071	123.5	-
${\text{Mg}}_{\text{10x10x10}}$	3,999 (39,990)	93,972,153	103.4	-	-	-

Table 13. (b) Benchmark system (ii): bcc Mo periodic supercells with a mono-vacancy.

System	Number of atoms	FE basis	DFT-FE	Plane-wave	QE	ABINIT
	(Number of electrons)			basis
${\text{Mo}}_{\text{6x6x6}}$	431 (6034)	5,475,843	0.5	194,310	0.56	0.7
${\text{Mo}}_{\text{8x8x8}}$	1,023 (14,322)	12,942,743	4.2	460,725	22.1	115.7
${\text{Mo}}_{\text{10x10x10}}$	1,999 (27,986)	25,229,995	17.7	899,849	219.5	-

Table 14. (c) Benchmark system (iii): Cu Icosahedron nano-particles of varying sizes.

System	Number of atoms	FE basis	DFT-FE	Plane-wave	QE	ABINIT
	(Number of electrons)			basis
${\text{Cu}}_{\text{3-Shell}}$	147 (2793)	6,584,861	0.3	1,080,751	0.2	0.8
${\text{Cu}}_{\text{4-Shell}}$	309 (5,871)	13,974,767	1.7	2,110,867	5.5	10.7
${\text{Cu}}_{\text{5-Shell}}$	561 (10,659)	26,060,299	5.3	3,647,655	63.4	-
${\text{Cu}}_{\text{6-Shell}}$	923 (17,537)	41,775,101	12.7	5,792,547	-	-

Table 15. Table 8: Minimum wall-time comparison of DFT-FE with QE: Average time per SCF iteration step in seconds (rounded to the nearest whole number).

System	Number of atoms	DFT-FE	QE
	(Number of electrons)
${\text{Mg}}_{\text{4x4x4}}$	255 (2,550)	19	29
${\text{Mg}}_{\text{6x6x6}}$	863 (8,630)	38	165
${\text{Mg}}_{\text{8x8x8}}$	2,047 (20,470)	91	816
${\text{Mg}}_{\text{10x10x10}}$	3,999 (39,990)	203	-

Table 16. (a) Benchmark system (i): hcp Mg periodic supercells with a mono-vacancy.

System	Number of atoms	DFT-FE	QE
	(Number of electrons)
${\text{Mg}}_{\text{4x4x4}}$	255 (2,550)	19	29
${\text{Mg}}_{\text{6x6x6}}$	863 (8,630)	38	165
${\text{Mg}}_{\text{8x8x8}}$	2,047 (20,470)	91	816
${\text{Mg}}_{\text{10x10x10}}$	3,999 (39,990)	203	-

Table 17. (b) Benchmark system (ii): bcc Mo periodic supercells with a mono-vacancy.

System	Number of atoms	DFT-FE	QE
	(Number of electrons)
${\text{Mo}}_{\text{6x6x6}}$	431 (6,034)	23	173
${\text{Mo}}_{\text{8x8x8}}$	1,023 (14,322)	52	549
${\text{Mo}}_{\text{10x10x10}}$	1,999 (27,986)	117	1883
${\text{Mo}}_{\text{13x13x13}}$	4,393 (61,502)	277	-

Table 18. (c) Benchmark system (iii): Cu Icosahedron nano-particles of varying sizes.

System	Number of atoms	DFT-FE	QE
	(Number of electrons)
${\text{Cu}}_{\text{3-Shell}}$	147 (2,793)	15	50
${\text{Cu}}_{\text{4-Shell}}$	309 (5,871)	25	183
${\text{Cu}}_{\text{5-Shell}}$	561 (10,659)	44	304
${\text{Cu}}_{\text{6-Shell}}$	923 (17,537)	70	-

Table 19. Table 9: Breakdown of average wall-time per SCF iteration step (in seconds, rounded to the nearest whole number) using DFT-FE for large systems into the following computational steps: a) ES (Total electrostatic potential solve), b) CF (Chebyshev filtering), c) CholGS (Cholesky-Gram-Schimdt Orthogonalization), d) RR (Rayleigh-Ritz procedure), e) DC (Electron-density computation) and f) O (other—Discrete Hamiltonian computation, electron density mixing, and computation of Fermi energy). “NDP” denotes number of domain decomposition MPI tasks, and “NBP” denotes number of band parallelization groups. Total number of MPI tasks is NDP times NBP.

System	No. atoms	No. electrons	DOF’s	NDP (NBP)	Total
			per atom		MPI tasks
${\text{Mg}}_{\text{8x8x8}}$	2,047	20,470	22,450	32,000 (2)	64,000
${\text{Mg}}_{\text{10x10x10}}$	3,999	39,990	23,499	51,200 (3)	153,600
${\text{Mo}}_{\text{10x10x10}}$	1,999	27,986	12,621	16,000 (3)	48,000
${\text{Mo}}_{\text{13x13x13}}$	4,393	61,502	12,594	48,000 (4)	192,000

Table 20. (a) Setup of the benchmark simulations.

System	No. atoms	No. electrons	DOF’s	NDP (NBP)	Total
			per atom		MPI tasks
${\text{Mg}}_{\text{8x8x8}}$	2,047	20,470	22,450	32,000 (2)	64,000
${\text{Mg}}_{\text{10x10x10}}$	3,999	39,990	23,499	51,200 (3)	153,600
${\text{Mo}}_{\text{10x10x10}}$	1,999	27,986	12,621	16,000 (3)	48,000
${\text{Mo}}_{\text{13x13x13}}$	4,393	61,502	12,594	48,000 (4)	192,000

Table 21. (b) Breakdown of average wall-time per SCF iteration step.

System	ES	CF	CholGS	RR	DC	O	Total
							time
${\text{Mg}}_{\text{8x8x8}}$	11	33	15	24	6	2	91
${\text{Mg}}_{\text{10x10x10}}$	12	63	53	63	9	3	203
${\text{Mo}}_{\text{10x10x10}}$	5	49	23	30	6	4	117
${\text{Mo}}_{\text{13x13x13}}$	7	80	80	97	9	4	277

Table 22. Table 10: Demonstration of geometry optimization in DFT-FE. Case study: tris (bipyridine) ruthenium (TBR) . FE ord , h min subscript FE ord subscript ℎ min {\rm FE}_{\rm ord},h_{\rm min} and h max subscript ℎ max h_{\rm max} denote the FE polynomial order, minimum element size and maximum element size (Bohr), respectively. E cut subscript E cut {\rm E}_{\rm cut} denotes plane-wave basis cut-off (Hartree), and E g subscript 𝐸 g {E}_{\rm g} denotes ground-state energy per atom (Hartree). f max = max 1 ≤ i ≤ N a ⁡ ‖ 𝐟 i ‖ subscript 𝑓 max subscript 1 𝑖 subscript 𝑁 𝑎 norm subscript 𝐟 𝑖 f_{\rm max}=\max\limits_{1\leq i\leq N_{a}}\left\|{\bf f}_{i}\right\| (Hartree/Bohr), where 𝐟 i subscript 𝐟 𝑖 {\bf f}_{i} denotes the force on the i th superscript 𝑖 th i^{\rm th} atom, and maximum difference in forces is measured by Δ max ​ f = max 1 ≤ i ≤ N a ⁡ ‖ 𝐟 i DFT − FE − 𝐟 i QE ‖ subscript Δ max 𝑓 subscript 1 𝑖 subscript 𝑁 𝑎 norm superscript subscript 𝐟 𝑖 DFT FE superscript subscript 𝐟 𝑖 QE \Delta_{\rm max}{f}=\max\limits_{1\leq i\leq N_{a}}\left\|{\bf f}_{i}^{\rm DFT-FE}-{\bf f}_{i}^{\rm QE}\right\| (Hartree/Bohr). Δ ​ E relax = E g , relax − E g , un − relaxed Δ subscript E relax subscript E g relax subscript E g un relaxed \Delta{\rm E}_{\rm relax}={\rm E}_{\rm g,relax}-{\rm E}_{\rm g,un-relaxed} denotes the difference in the ground-state energy per atom between the relaxed and un-relaxed TBR structures (Hartree). Maximum difference between relaxed and the un-relaxed TBR coordinates in DFT-FE is measured by Δ max ​ R relax = max 1 ≤ i ≤ N a ⁡ ‖ 𝐑 i relaxed − 𝐑 i un − relaxed ‖ subscript Δ max superscript 𝑅 relax subscript 1 𝑖 subscript 𝑁 𝑎 norm superscript subscript 𝐑 𝑖 relaxed superscript subscript 𝐑 𝑖 un relaxed \Delta_{\rm max}{R^{\rm relax}}=\max\limits_{1\leq i\leq N_{a}}\left\|{\bf R}_{i}^{\rm relaxed}-{\bf R}_{i}^{\rm un-relaxed}\right\| (Bohr).

DFT-FE	DFT-FE	DFT-FE	QE	QE	QE	${\Delta }_{\max }f$
$({\mathrm{FE}}_{\mathrm{ord}},h_{\min },h_{\max })$	$E_g$	$f_{\max }$	$E_{\mathrm{cut}}$	$E_g$	$f_{\max }$
4, 0.42, 13.33	-5.5566947	0.1626628	50	-5.5566265	0.1626344	$1.1\times {10}^{-4}$

Table 23. (a) Comparison of energies and forces in the un-relaxed TBR structure between DFT-FE and QE.

DFT-FE	DFT-FE	DFT-FE	QE	QE	QE	${\Delta }_{\max }f$
$({\mathrm{FE}}_{\mathrm{ord}},h_{\min },h_{\max })$	$E_g$	$f_{\max }$	$E_{\mathrm{cut}}$	$E_g$	$f_{\max }$
4, 0.42, 13.33	-5.5566947	0.1626628	50	-5.5566265	0.1626344	$1.1\times {10}^{-4}$

Table 24. (b) Comparison between relaxed and un-relaxed TBR structrues.

DFT-FE	QE	${\Delta }_{\max }{R}^{\mathrm{relax}}$
$\Delta E_{\mathrm{relax}}$	$\Delta E_{\mathrm{relax}}$
-0.0015992	-0.0016009	0.276

Equations

$$\left(-\frac{1}{2}{\nabla }^2+V_{\text{eff}}(\rho ,\text{R})\right){\psi }_i={ϵ}_i{\psi }_i,i=1,2,\cdots N\text{with}N>\frac{N_e}{2},$$

$$\rho (\text{x})=2\sum _{i=1}^{N}f({ϵ}_i,\mu )|{\psi }_i(\text{x}){|}^2.$$

$$f(ϵ,\mu )=\frac{1}{1+\exp \left(\frac{ϵ-\mu }{\sigma }\right)}.$$

$$\int \rho (\text{x})d\text{x}=2\sum _{i}f({ϵ}_i,\mu )=N_e.$$

$$V_{\text{eff}}(\rho ,\text{R})=V_{\text{xc}}(\rho )+V_{\text{el}}(\rho ,\text{R})=\frac{\delta E_{\text{xc}}}{\delta \rho }+\frac{\delta E_{\text{el}}}{\delta \rho }$$

$$E_{\text{xc}}(\rho )=\int {ϵ}_{\text{xc}}(\rho ,\nabla \rho )\rho (\text{x})d\text{x}.$$

$$E_{\text{el}}(\rho ,\text{R})=E_{\text{H}}(\rho )+E_{\text{ext}}(\rho ,\text{R})+E_{\text{zz}}(\text{R}).$$

$$E_{\text{H}}=\frac{1}{2}\int \int \frac{\rho (\text{x})\rho (\text{y})}{|\text{x}-\text{y}|}d\text{x}d\text{y},E_{\text{ext}}=-\sum _J\int \rho (\text{x})\frac{Z_J}{|\text{x}-{\text{R}}_J|}d\text{x},E_{\text{zz}}=\frac{1}{2}\sum _{\begin{array}{c}I,J\ne I\end{array}}\frac{Z_I{Z}_J}{|{\text{R}}_I-{\text{R}}_J|},$$

$$\begin{array}{rl}{E}_{zz}& =\frac{1}{2}\int \int \frac{b(\text{x},\text{R})b(\text{y},\text{R})}{|\text{x}-\text{y}|}d\text{x}d\text{y}-E_{\text{self}}\\ & =\frac{1}{2}\int \int \left(\frac{b(\text{x},\text{R})b(\text{y},\text{R})}{|\text{x}-\text{y}|}-\sum _I\frac{Z_I^2\tilde{\delta }(|\text{x}-{\text{R}}_I|)\tilde{\delta }(|\text{y}-{\text{R}}_I|)}{|\text{x}-\text{y}|}\right)d\text{x}d\text{y},\end{array}$$

$$E_{\text{ext}}=2\sum _{i=1}^N\int \int f_i{\psi }_i^{\ast }(\text{x})V_{\text{PS}}(\text{x},\text{y},\text{R}){\psi }_i(\text{y})d\text{y}d\text{x}.$$

$$V_{\text{nl}}{\psi }_i:=\int V_{\text{nl}}(\text{x},\text{y},\text{R}){\psi }_i(\text{y})d\text{y}$$

$$\text{with}{C}_{lpm}^{J,i}=\int {\chi }_{lpm}^J(\text{x},{\text{R}}_J){\psi }_i(\text{x})d\text{x},\frac{1}{h_{lp}^J}=\text{innerproduct}{\xi }_{lm}^J{\chi }_{lpm}^J.$$

$${\bar{V}}_{\tilde{\delta }}^I(\text{x})=\int \frac{-Z_I\tilde{\delta }(|\text{y}-{\text{R}}_I|)}{|\text{x}-\text{y}|}d\text{y},\phi (\text{x},\text{R})=\int \frac{\rho (\text{y})+b(\text{y},\text{R})}{|\text{x}-\text{y}|}d\text{y}.$$

$$\begin{array}{rl}{E}_{\text{el}}& =\frac{1}{2}\int (\rho (\text{x})+b(\text{x},\text{R}))\phi (\text{x},\text{R})d\text{x}-\frac{1}{2}\sum _I\int -Z_I\tilde{\delta }(|\text{x}-{\text{R}}_I|){\bar{V}}_{\tilde{\delta }}^I(\text{x})d\text{x}\\ & +\sum _J\int \left(V_{\text{loc}}^J(|\text{x}-{\text{R}}_J|)-{\bar{V}}_{\tilde{\delta }}^J(|\text{x}-{\text{R}}_J|)\right)\rho (\text{x})d\text{x}+2\sum _{i=1}^N\int \int f_i{\psi }_i^{\ast }(\text{x})V_{\text{nl}}(\text{x},\text{y}){\psi }_i(\text{y})d\text{y}d\text{x}.\end{array}$$

$$\left(-\frac{1}{2}{\nabla }^2+V_{\text{xc}}+\phi +\sum _J(V_{\text{loc}}^J-{\bar{V}}_{\tilde{\delta }}^J)+V_{\text{nl}}\right){\psi }_i={ϵ}_i{\psi }_i,$$

$$-\frac{1}{4\pi }{\nabla }^2\phi (\text{x},\text{R})=\rho (\text{x})+b(\text{x},\text{R}),-\frac{1}{4\pi }{\nabla }^2{\bar{V}}_{\tilde{\delta }}^I(\text{x},{\text{R}}_I)=-Z_I\tilde{\delta }(|\text{x}-{\text{R}}_I|),$$

$$2\sum _{i}f({ϵ}_i,\mu )=N_e,\rho (\text{x})=2\sum _{i}f({ϵ}_i,\mu )|{\psi }_i(\text{x}){|}^2.$$

$${\psi }_{n,\mathbf{k}}(\text{x})=e^{i\mathbf{k}\cdot \text{x}}{u}_{n,\mathbf{k}}(\text{x}),$$

$$\rho (\text{x})=2\sum _{n=1}^N\text{fint}{\text{displaylimits}}_{BZ}{f}_{n,\mathbf{k}}|u_{n,\mathbf{k}}(\text{x}){|}^{2}d\mathbf{k},$$

$$E_{\text{ext}}^{\text{nl}}=2\sum _{n=1}^N\text{fint}{\text{displaylimits}}_{BZ}\int {\text{displaylimits}}_{{\Omega }_{\text{p}}}\int {\text{displaylimits}}_{{\mathbb{R}}^3}{f}_{n,\mathbf{k}}{u}_{n,\mathbf{k}}^{\ast }(\text{x})e^{-i\mathbf{k}\cdot \text{x}}{V}_{\text{nl}}(\text{x},\text{y},\text{R})e^{i\mathbf{k}\cdot \text{y}}{u}_{n,\mathbf{k}}(\text{y})d\text{y}d\text{x}d\mathbf{k}.$$

$$V_{\text{nl}}{u}_{n,\mathbf{k}}:=e^{-i\mathbf{k}\cdot \text{x}}\int {\text{displaylimits}}_{{\mathbb{R}}^3}{V}_{\text{nl}}(\text{x},\text{y},\text{R})e^{i\mathbf{k}\cdot \text{y}}{u}_{n,k}(\text{y})d\text{y}=\sum _a\sum _{lpm}\sum _r{e}^{-i\mathbf{k}.(\text{x}-{\text{L}}_r)}{C}_{lpm}^{a,n\mathbf{k}}{h}_{lp}^a{\chi }_{lpm}^a(\text{x},{\text{R}}_a+{\text{L}}_r),$$

$$C_{lpm}^{a,n\mathbf{k}}=\int {\text{displaylimits}}_{{\Omega }_{\text{p}}}\sum _r{e}^{i\mathbf{k}.(\text{x}-{\text{L}}_r)}{\chi }_{lpm}^a(\text{x},{\text{R}}_a+{\text{L}}_r)u_{n,\mathbf{k}}(\text{x})d\text{x}.$$

$$\left(-\frac{1}{2}{\nabla }^2-i\mathbf{k}.\nabla +\frac{1}{2}|\mathbf{k}{|}^2+V_{\text{xc}}+\phi +\sum _J(V_{\text{loc}}^J-{\bar{V}}_{\tilde{\delta }}^J)+V_{\text{nl}}\right)u_{n,\mathbf{k}}={ϵ}_{n,\mathbf{k}}{u}_{n,\mathbf{k}}\text{on}{\Omega }_p,$$

$$-\frac{1}{4\pi }{\nabla }^2\phi (\text{x},\text{R})=\rho (\text{x})+b(\text{x},\text{R})\text{on}{\Omega }_p,-\frac{1}{4\pi }{\nabla }^2{\bar{V}}_{\tilde{\delta }}^I(\text{x},{\text{R}}_I)=-Z_I\tilde{\delta }(|\text{x}-{\text{R}}_I|)\text{on}{\mathbb{R}}^3,$$

$$2\sum _{n=1}^N\text{fint}{\text{displaylimits}}_{BZ}{f}_{n,\mathbf{k}}d\mathbf{k}=N_e,\rho (\text{x})=2\sum _{n=1}^N\text{fint}{\text{displaylimits}}_{BZ}{f}_{n,\mathbf{k}}|u_{n,\mathbf{k}}(\text{x}){|}^{2}d\mathbf{k}.$$

$${\mathcal{F}}_0(\text{R})=\underset{\mathbf{f}\in [0,1{]}^N}{\min }\underset{\mathbf{\Psi }\in (\mathcal{Y}{)}^N}{\min }\underset{\phi \in \mathcal{Y}}{\max }\mathcal{L}(\mathbf{f},\mathbf{\Psi },\phi ;\text{R})\text{such that}\int {\psi }_i^{\ast }{\psi }_{j}d\text{x}={\delta }_{ij},2\sum _i{f}_i=N_e,$$

$$\text{where}\mathcal{L}(\mathbf{f},\mathbf{\Psi },\phi ;\text{R})=\tilde{\mathcal{L}}(\mathbf{f},\mathbf{\Psi })+\underset{\mathcal{V}\in (H^1({\mathbb{R}}^3){)}^{N_a}}{\min }{\mathcal{L}}_{\text{el}}(\mathbf{f},\mathbf{\Psi },\phi ,\mathcal{V};\text{R}),$$

$$\text{with}\tilde{\mathcal{L}}(\mathbf{f},\mathbf{\Psi })=T_{\text{s}}(\mathbf{f},\mathbf{\Psi })+E_{\text{xc}}(\rho )+E_{\text{ent}}(\mathbf{f}).$$

$$T_{\text{s}}(\mathbf{f},\mathbf{\Psi })=2\sum _{i=1}^N\int f_i{\psi }_i^{\ast }(\text{x})\left(-\frac{1}{2}{\nabla }^2\right){\psi }_i(\text{x})d\text{x},$$

$$E_{\text{ent}}=-2\sigma \sum _{i=1}^N\left[f_i\ln f_i+(1-f_i)\ln (1-f_i)\right].$$

$$\begin{array}{rl}& {\mathcal{L}}_{\text{el}}(\mathbf{f},\mathbf{\Psi },\phi ,\mathcal{V},\text{R})=\int \left[-\frac{1}{8\pi }|\nabla \phi (\text{x}){|}^2+(\rho (\text{x})+b(\text{x},\text{R}))\phi (\text{x})\right]d\text{x}\\ & +\sum _I\int \left[\frac{1}{8\pi }|\nabla V^I(\text{x}){|}^2+Z_I\tilde{\delta }(|\text{x}-{\text{R}}_I|)V^I(\text{x})\right]d\text{x}+\sum _J\int \left(V_{\text{loc}}^J(|\text{x}-{\text{R}}_J|)-{\bar{V}}_{\tilde{\delta }}^J(|\text{x}-{\text{R}}_J|)\right)\rho (\text{x})d\text{x}\\ & +2\sum _{i=1}^N\int \int f_i{\psi }_i^{\ast }(\text{x})V_{\text{nl}}(\text{x},\text{y},\text{R}){\psi }_i(\text{y})d\text{y}d\text{x},\end{array}$$

$$\frac{d\mathcal{F}_{0}(\,\bm{\tau}^{\varepsilon})}{d\varepsilon}\bigg{|}_{\varepsilon=0}=\int\displaylimits_{\Omega}\bm{\textbf{E}}:\nabla\bm{\Upsilon}(\bm{\textbf{x}})\,d\bm{\textbf{x}}+\sum_{I}\int\displaylimits_{{\mathbb{R}^{3}}}\bm{\textbf{E}}^{\prime I}:\nabla\bm{\Upsilon}(\bm{\textbf{x}})\,d\bm{\textbf{x}}\,+\,\text{F}^{\text{PSP}}\,,$$

$$\displaystyle\bm{\textbf{E}}=\Biggl{(}\sum_{i=1}^{N}\Bigl{(}\bar{f}_{i}\nabla\bar{\psi}_{i}^{*}(\bm{\textbf{x}})\cdot\nabla\bar{\psi}_{i}(\bm{\textbf{x}})-2\,\bar{f}_{i}\,\bar{\epsilon}_{i}\,\bar{\psi}^{*}_{i}(\bm{\textbf{x}})\;\bar{\psi}_{i}(\bm{\textbf{x}})\Bigr{)}+\varepsilon_{\text{xc}}(\bar{\rho},\nabla\bar{\rho})\bar{\rho}(\bm{\textbf{x}})-\frac{1}{8\pi}|\nabla\bar{\varphi}(\bm{\textbf{x}})|^{2}+\bar{\rho}(\bm{\textbf{x}})\bar{\varphi}(\bm{\textbf{x}})\,$$

Full text

††journal: Computer Physics Communications

Program summary

Program Title: DFT-FE (https://github.com/dftfeDevelopers/dftfe)

Journal Reference:

Catalogue identifier:

Licensing provisions: LGPL v3

Programming language: C/C++

Computer: Any system with C/C++ compiler and MPI library.

Operating system: Linux

RAM: Ranges from few GBs for small problem sizes to around 50,000 GB for a system with 61,502 electrons.

Number of processors used: Range from 64 to 192,000 MPI tasks (demonstrated).

Keywords: Electronic structure, Real-space, Large-scale, Spectral finite-elements, Pseudopotential, All-electron.

Classification: 0.7

External routines/libraries: p4est (http://www.p4est.org/), deal.II (https://www.dealii.org/),

BLAS (http://www.netlib.org/blas/), LAPACK (http://www.netlib.org/lapack/),

ELPA (https://elpa.mpcdf.mpg.de/), ScaLAPACK (http://www.netlib.org/scalapack/),

Spglib (https://atztogo.github.io/spglib/), ALGLIB (http://www.alglib.net/),

LIBXC (http://www.tddft.org/programs/libxc/), PETSc (https://www.mcs.anl.gov/petsc),

SLEPc (http://slepc.upv.es).

Nature of problem: Density functional theory calculations.

Solution method: We employ a local real-space variational formulation of Kohn-Sham density functional theory that is applicable for both pseudopotential and all-electron calculations with arbitrary boundary conditions. Higher-order adaptive spectral finite-element basis is used to discretize the Kohn-Sham equations. Chebyshev polynomial filtered subspace iteration procedure (ChFSI) is employed to solve the nonlinear Kohn-Sham eigenvalue problem self-consistently. ChFSI in DFT-FE employs Cholesky factorization based orthonormalization, and spectrum splitting based Rayleigh-Ritz procedure in conjunction with mixed precision arithmetic. Configurational force approach is used to compute ionic forces and periodic cell stresses for geometry optimization.

Restrictions: Exchange correlation functionals are restricted to Local Density Approximation (LDA) and Generalized Gradient Approximation (GGA), with and without spin. The pseudopotentials available are optimized norm conserving Vanderbilt (ONCV) pseudopotentials and Troullier–Martins (TM) pseudopotentials. Calculations are non-relativistic.

Unusual features: DFT-FE handles all-electron and pseudopotential calculations in the same framework, while accommodating periodic, non-periodic and semi-periodic boundary conditions.

Running time: This is dependent on problem type and computational resources used. Timing results for benchmark systems are provided in the paper.

1 Introduction

Quantum mechanical calculations based on Kohn-Sham density functional theory (DFT) occupy a substantial fraction of world’s computational resources today, and have provided many important insights into a wide range of materials properties over the past decade. The Kohn-Sham approach [1; 2] to DFT provides an efficient formulation to compute the ground-state properties of materials systems by reducing the many-body Schrödinger problem of interacting electrons into an equivalent problem of non-interacting electrons in an effective mean field that is governed by the electron-density. While this effective single-electron formulation has no approximations, and is exact, in principle, the quantum mechanical interactions between electrons manifest in the form of an unknown exchange-correlation functional, which is modeled in practice. While improving the accuracy of the exchange-correlation functional is still an active area of research, the widely used models for the exchange-correlation functional [3] have been shown to predict a range of materials properties across various materials systems with good accuracy.

Large-scale DFT calculations are crucial to improve the predictive capability of modeling materials systems, and enable computation based design of new materials in a variety of application areas. For example, accurate determination of dislocation core properties [4; 5; 6; 7; 8; 9; 10; 11; 12; 13] in metals and semiconductors, understanding ion conduction mechanisms and computing diffusivities in solid-state electrolytes [14; 15], and large-scale bio-molecular simulations [16], all require the ability to conduct accurate and computationally efficient DFT calculations involving many thousands of atoms on both metallic and insulating systems. However, the solution to the governing equations in DFT demands significant computational resources, and accurate DFT calculations are routinely limited to materials systems with at most a few thousands of electrons restricting the system sizes to a few hundred atoms. This computational complexity is a serious bottleneck in the context of ab-initio molecular dynamics with long time-scales and atomic relaxations that require a large number of relaxation steps.

Traditionally, the widely used DFT codes employ either plane-waves [17; 18; 19; 20] or atomic-orbital type basis sets [21; 22; 23; 24; 25] for DFT calculations. However, the use of plane-wave basis restricts simulation domains to be periodic. Further, these basis sets do not exhibit good parallel scalability, severely limiting the range of materials systems that can be studied. On the other hand, atomic orbital type basis sets are not systematically convergent for generic materials systems. Thus, to overcome the above limitations, there has been an increasing thrust in the development of systematically improvable and scalable real-space discretization techniques like finite-elements [26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39], finite-difference [40; 41; 42; 43; 44; 45], wavelets [46], psinc functions [47], and other reduced order basis techniques  [48; 49].

We note that the traditional self-consistent approaches to solving the discretized nonlinear Kohn-Sham eigenvalue problem involves the diagonalization of the Kohn-Sham Hamiltonian to obtain the orthonormal eigenvectors, the computational complexity of which typically scales cubically with number of electrons. Hence, the computational cost of a Kohn-Sham DFT calculations becomes prohibitively expensive as the system size grows larger, and is another limiting factor in realizing large-scale DFT calculations. To circumvent this limitation, numerous efforts have focused on reducing the computational complexity by improving the system-size scaling of DFT calculations. These methods [50; 51; 47; 52; 53] rely on the exponential decay of the density matrix in real space, and the computational cost has been demonstrated to linearly scale with number of atoms for systems with non-vanishing band gaps. Though reduced scaling approaches for metallic systems at moderate electronic temperatures have been developed [54; 55; 56], they usually have a high computational prefactor. There are also other efforts in the literature which have focused on reducing the prefactor associated with the cubic computational complexity of the Kohn-Sham DFT calculations [48; 57; 58] that are suited for both insulating and metallic systems.

In this work, we focus on real-space adaptive spectral finite-element (FE) discretization of Kohn-Sham DFT which affords excellent parallel scalability. We employ Chebyshev filtering approach  [59] in conjunction with Cholesky factorization based orthonormalization procedure and spectrum splitting based Rayleigh-Ritz technique, along with mixed precision strategies, to reduce the computational prefactor of the cubic scaling steps. These strategies enabled systematically convergent, computationally efficient and massively parallel DFT calculations (demonstrated up to $192,000$ MPI tasks) on material systems with tens of thousands of electrons for both metallic and insulating systems. The choice of FE discretization among other real space discretizations for DFT in this work is motivated by some key advantages it offers for electronic structure calculations. In particular, the FE basis naturally allows for arbitrary boundary conditions, provides good scalability on parallel computing platforms due to locality of the basis, and is amenable to adaptive spatial resolution, which can effectively be exploited for efficient solution of all electron DFT calculations as well as the development of coarse-graining techniques that seamlessly bridge electronic structure calculations with continuum.

We present here DFT-FE, a massively parallel adaptive finite-element code for large-scale density functional theory calculations. The current endeavour extends the previous work  [34], which has demonstrated the advantage of higher-order spectral finite-elements in significantly improving the computational efficiency of DFT calculations. DFT-FE is based on local real-space variational formulation of Kohn-Sham DFT [34; 60], which handles all-electron and pseudopotential calculations in the same framework while accommodating periodic, non-periodic and semi-periodic boundary conditions. This unified local variational formulation is briefly discussed and the usefulness of this formulation to compute configurational forces [60] for geometry optimization is also highlighted. These configurational forces correspond to generalized variational forces computed as the derivative of the Kohn-Sham energy functional with respect to the position of a material point x. These generalized forces that result from the inner variations of the Kohn-Sham energy functional inherently account for the Pulay corrections, and provide a unified framework to compute ionic forces as well as stress tensor for geometry optimization. We then introduce the FE discretization of Kohn-Sham DFT, and the resulting discretized Kohn-Sham equations are cast into a standard eigenvalue problem with modest computational cost by using a spectral FE basis in conjunction with the Gauss-Lobatto-Legendre (GLL) quadrature rules  [34]. We subsequently discuss various adaptive mesh refinement procedures implemented in DFT-FE. Firstly, we employ an user defined adaptive mesh refinement procedure, involving user inputs for mesh sizes at different regions of interest in the simulation domain, which in turn can be estimated from a mesh size distribution $h(\bm{\textbf{x}})$ obtained from a constrained minimization of the discrete ground-state energy error estimate. Further, we also present an automatic mesh refinement procedure to construct an a priori adaptive mesh, by employing numerically constructed single atom Kohn-Sham DFT wavefunctions. This approach is based on local error indicators obtained from the energy error estimates involving the discrete electronic wavefunctions.

We employ Chebyshev filtered subspace iteration technique (ChFSI) [59; 34] to self-consistently solve the discretized Kohn-Sham DFT problem. ChFSI technique involves three key computational steps: (i) Chebyshev filtering to compute a subspace that closely approximates the relevant eigensubspace; (ii) orthonormalization procedure to construct an orthonormal basis spanning the Chebyshev filtered subspace; and (iii) Rayleigh-Ritz procedure that involves projecting the Kohn-Sham Hamiltonian into the Chebyshev filtered subspace, followed by diagonalization of the projected Hamiltonian to compute the associated eigenpairs and a subspace rotation step. The computational complexity of Chebyshev filtering scales quadratically with system size, and is the dominant cost for small- to medium-scale system sizes. However, for large-scale systems, the cubic-scaling computational cost of orthonormalization and Rayleigh-Ritz procedure dominates. To this end, the numerical implementation in DFT-FE focuses on reducing the computational prefactor by using efficient numerical strategies, which include: (i) optimized FE cell level matrix operations; (ii) Cholesky factorization based Gram-Schmidt orthonormalization; (iii) mixed precision arithmetic, and (iv) spectrum splitting based Rayleigh-Ritz procedure. We note that the use of these techniques delays the onset of cubic computational complexity in DFT-FE to very large system sizes. Using benchmark problems, we demonstrate that DFT-FE exhibits close to quadratic scaling in computational complexity till $30,000$–$40,000$ electrons.

In our numerical implementation of DFT-FE, we use two levels of parallelization: (i) domain-decomposition based on partitioning of the adaptive FE mesh, and (ii) parallelization over wavefunctions (band parallelization). Our numerical investigations show that best parallel scaling efficiency is achieved using domain-decomposition combined with moderate band parallelization. We assess the parallel scaling performance of DFT-FE on a range of system sizes, and observe excellent scalability on all systems. Notably, we obtain a scaling efficiency of $\sim 50\%$ on 102,400 MPI tasks for a system containing $\sim 40,000$ electrons. We note that DFT-FE’s massive parallel scalability is a result of the locality of the FE basis as well as an effective parallel implementation of the various algorithms that reduce communication costs and latency.

In order to assess the accuracy, computational efficiency and scalability of DFT-FE, we conduct a comprehensive comparison study with widely used plane-wave codes—QUANTUM ESPRESSO (QE) [17; 61] and ABINIT [18]—on various pseudopotential DFT benchmark problems, which include periodic bulk metallic systems with a defect and non-periodic nano-particles. In the validation studies, we compare DFT-FE ground-state energies, ionic forces and cell stresses against QE results at different accuracy levels. We also consider all-electron periodic and non-periodic benchmark systems for assessing the accuracy of DFT-FE with exciting [20] and NWChem [25] codes. Next, we assess the computational efficiency and scalability of pseudopotential DFT-FE calculations with respect to QE and ABINIT. The FE and plane-wave discretizations in these benchmark calculations are chosen to be commensurate with chemical accuracy ($\sim 10^{-4}$ Ha/atom and $\sim 10^{-4}$ Ha/Bohr in ground-state energies and forces). The CPU-times on these benchmark systems show that DFT-FE is computationally more efficient than QE, even for periodic systems, beyond system sizes containing $3,000$ electrons, and significantly more efficient for larger system sizes by $4.5-12\times$. Comparing the minimum wall times, we find DFT-FE to be significantly faster than QE for all the systems considered here, with speedups reaching up to $16\times$ for larger systems. Finally, we also demonstrate the capability of DFT-FE to simulate very large systems reaching up to $\sim 60,000$ electrons, which are computationally prohibitive using plane-wave codes. For the largest benchmark system containing $\sim 60,000$ electrons, which is solved using a FE discretization commensurate with chemical accuracy, we scale up to $192,000$ MPI tasks at $42\%$ efficiency and consume an average of $4.6$ minutes per SCF iteration.

Finally, we demonstrate the following capabilities in DFT-FE: i) ionic relaxations using an organometallic complex as the benchmark system, which is validated using QE ii) NVE ab-initio molecular dynamics simulations using a representative bulk-metallic system containing $108$ atoms, which demonstrates energy conservation to very stringent accuracy, and iii) band-structure calculations of bulk Si and Si supercell with a vacancy, which are validated using QE. The accuracy, efficiency and scalability of the DFT-FE code demonstrated in this work, presents this as a useful code to conduct systematically convergent large-scale DFT calculations in an efficient manner, which can enable studies on complex materials systems that have not been possible heretofore. We note that the present release of DFT-FE is only a first step in a longer term effort, which includes future planned activities of incorporating advanced exchange-correlation functionals, Hamiltonians with spin-orbit coupling, dielectric calculations, electron-phonon coupling, and others capabilities into DFT-FE.

The remainder of the paper is organized as follows. Section 2 discusses the governing equations in DFT, describes the local real space formulation implemented in DFT-FE that treats pseudopotential and all-electron calculations in the single framework, followed by a brief discussion on configurational forces used to evaluate ionic forces and stresses for geometry optimization. FE discretization of the Kohn-Sham DFT problem is subsequently introduced with a discussion on the discrete Kohn-Sham equations. Section 3 describes the various aspects of numerical implementation in DFT-FE, which include adaptive mesh refinement strategies implemented in DFT-FE, and numerical implementation of the various algorithms employed in DFT-FE for the solution of the Kohn-Sham problem. Section 4 presents the results on various benchmark systems, demonstrating the accuracy, computational efficiency and parallel scaling performance of the developed DFT-FE code in comparison with other widely used DFT codes. Section 5 highlights other capabilities of the DFT-FE code, and we conclude with an outlook in Section 6.

2 Real space Kohn-Sham DFT formulation

2.1 Governing equations in DFT

We consider a materials system with $N_{e}$ electrons and $N_{a}$ atoms whose position vectors are denoted by $\bm{\textbf{R}}=\{\bm{\textbf{R}}_{1},\,\bm{\textbf{R}}_{2},\,\cdots\bm{\textbf{R}}_{N_{a}}\}$. Neglecting spin, the variational problem of evaluating the ground-state properties in density functional theory is equivalent to solving the $N$ lowest eigenvalues of the following non-linear eigenvalue problem [1]:

[TABLE]

where $\epsilon_{i}$ and $\psi_{i}$ denote the eigenvalues and corresponding eigenfunctions (also referred to as the canonical wavefunctions) of the Hamiltonian, respectively. For clarity and notational convenience, the case of spin independent Hamiltonian is discussed here. However, extension to spin-dependent Hamiltonians [62] is straightforward, and incorporated in DFT-FE. The electron density $\rho$ in equation  (1) can be expressed in terms of the orbital occupancy function $f(\epsilon,\mu)$ and the canonical wavefunctions as

[TABLE]

The range of $f(\epsilon_{i},\mu)$ lies in the interval $\left[0,1\right]$, and $\mu$ represents the Fermi-energy. In material systems with large number of eigenstates around the Fermi energy, the numerical instabilities that may arise in the solution of the non-linear Kohn-Sham eigenvalue problem are avoided by using a smooth orbital occupancy function. In DFT-FE, $f$ is represented by the Fermi-Dirac distribution [19; 50] given by

[TABLE]

In the above, $\sigma=k_{B}T$ denotes the regularization parameter with $k_{B}$ denoting the Boltzmann constant and $T$ representing the finite temperature. We note that as $\sigma\to 0$, the Fermi-Dirac distribution tends to the Heaviside function. The constraint on the total number of electrons in the system ($N_{e}$) determines the Fermi-energy $\mu$, and is given by

[TABLE]

We note that $f(\epsilon_{i},\mu)$ is denoted as $f_{i}$ in the remainder of the manuscript. The effective single-electron potential, $V_{\text{eff}}(\rho,\bm{\textbf{R}})$, in the Hamiltonian in equation (1) is given by

[TABLE]

In the above, $V_{\text{xc}}(\rho)$ denotes the exchange-correlation potential that accounts for quantum-mechanical interactions between electrons [3], and is given by the first variational derivative of the exchange-correlation energy $E_{\text{xc}}$. We adopt the generalized gradient approximation (GGA) [62; 63] for the exchange correlation functional description throughout the manuscript. However other forms of functionals involving local density (LDA, LSDA) are also incorporated in DFT-FE. In the case of GGA, the exchange-correlation energy is given by

[TABLE]

Numerous forms for $\epsilon_{\text{xc}}(\rho,\nabla\rho)$ have been proposed, and the three widely used forms are Becke (B88)  [64], Perdew and Wang (PW91)  [65] and Perdew, Burke and Enzerhof (PBE)  [66].

The term $V_{\text{el}}(\rho)$, in the effective single-electron potential (equation (5)), accounts for the electrostatic interactions. In particular, it is the variational derivative of the classical electrostatic interaction energy between electrons and nuclei, $E_{\text{el}}$, which can further be decomposed as

[TABLE]

In the above, $E_{\text{H}}$, $E_{\text{ext}}$ and $E_{\text{zz}}$ denote the electrostatic interaction energy between electrons (Hartree energy), interaction energy between nuclei and electrons, and repulsive energy between nuclei, respectively. These are given by

[TABLE]

where $Z_{I}$ denotes the charge on the $I^{th}$ nucleus. In the case of non-periodic boundary conditions representing an isolated atomic system, all integrals in equations (8) are over ${\mathbb{R}^{3}}$, and the summations include all the atoms in the system. In the case of periodic boundary conditions representing an infinite periodic crystal, all integrals involving x in equation (8) are over the periodic domain (supercell), whereas the integrals involving y are over ${\mathbb{R}^{3}}$. Further, the summation over $I$ is on atoms in the periodic domain, and the summation over $J$ extends over all the lattice sites. Henceforth, unless otherwise specified, we will adopt this convention. Next, we define the nuclear charge distribution $b(\bm{\textbf{x}},\bm{\textbf{R}})=-\sum_{I}Z_{I}\tilde{\delta}(|\bm{\textbf{x}}-\bm{\textbf{R}}_{I}|)$ with $\tilde{\delta}(\bm{\textbf{x}}-\bm{\textbf{R}}_{I})$ denoting a regularized Dirac distribution centered at $\bm{\textbf{R}}_{I}$ (and similarly $b(\bm{\textbf{y}},\bm{\textbf{R}})=-\sum_{J}Z_{J}\tilde{\delta}(|\bm{\textbf{y}}-\bm{\textbf{R}}_{J}|)$) to reformulate the repulsive energy $E_{\text{zz}}(\bm{\textbf{R}})$ as

[TABLE]

where $E_{\text{self}}$ denotes the self energy of the nuclear charges which depends only on the nuclear charge distribution.

The tightly bound core electrons close to the nucleus of an atom do not influence chemical bonding in many materials systems, and, thus, may not play a significant role in governing many materials properties. Hence, it is a common practice to adopt the pseudopotential approach, where valence electronic wavefunctions are computed in an effective potential generated by the the nucleus and the core electrons. The pseudopotential is often defined by the operator $\mathcal{V}_{\text{PS}}=\mathcal{V}_{\text{loc}}+\mathcal{V}_{\text{nl}}$ = $\sum_{J}(\mathcal{V}_{\text{loc}}^{J}+\mathcal{V}_{\text{nl}}^{J})$, where $\mathcal{V}_{\text{loc}}^{J}$ and $\mathcal{V}_{\text{nl}}^{J}$ denote the local and non-local part of the pseudopotential operator for an atom $J$, respectively. Further, in the case of norm-conserving pseudopotentials, $\mathcal{V}_{\text{nl}}^{J}$ can be constructed as a separable pseudopotential operator  [62; 67] of the form $\sum_{lpm}\ket{\chi_{lpm}}h_{lp}\bra{\chi_{lpm}}$ , with $\ket{\chi_{lpm}}$ denoting the pseudopotential projector. Here $l$ denotes the azimuthal quantum number, $p$ denotes the index corresponding to the projector component for a given $l$ while $m$ denotes the magnetic quantum number with $h_{lp}$ denoting the pseudopotential constant. Using the representation of the operator $\mathcal{V}_{PS}$ in the x basis, $E_{\text{ext}}$ in pseudopotential Kohn-Sham DFT can be expressed as

[TABLE]

Norm conserving pseudopotentials are employed in DFT-FE, where the action of the nonlocal psuedopotential operator on a wavefunction is given by

[TABLE]

[TABLE]

In the above, $\ket{\xi_{lm}^{J}}$ denotes the single atom pseudo-wavefunction. Note that $h^{J}_{lp}$ does not depend on the magnetic quantum number $m$ as the spherical harmonics associated with angular variables in the inner product  (12) are normalized to unity. We remark that equation  (11) reduces to Troullier-Martins (TM) pseudopotential [68] in the Kleinman-Bylander form [67] for one projector component, i.e. $p=1$ for every $l$, while in the case of optimized norm conserving Vanderbilt pseudopotential (ONCV)  [69] there are two projector components ($p=1,\,2$) for every $l$. Both TM and ONCV norm-conserving pseudopotentials are implemented in DFT-FE. We further note that the accuracy of ONCV pseudopotentials are shown to be on par with PAW approaches widely employed in DFT codes [70].

We note that the various components of the electrostatic interaction energy in  (8) and  (9) are non-local in real-space, and these extended interactions are reformulated as local variational problems as discussed in  [71; 34]. To this end, we define the electrostatic potential corresponding to the $I^{th}$ nuclear charge $Z_{I}\tilde{\delta}(|\bm{\textbf{x}}-\bm{\textbf{R}}_{I}|)$ to be $\bar{V}^{I}_{\tilde{\delta}}(\bm{\textbf{x}})$ and the electrostatic potential corresponding to the total charge distribution $(\rho+b)$ to be $\varphi(\bm{\textbf{x}},\bm{\textbf{R}})$, and these potentials are given by:

[TABLE]

Noting that the kernel corresponding to these extended interactions is the Green’s function of the Laplace operator, these potentials can be efficiently computed by taking recourse to the solution of a Poisson problem. Using the potentials defined in  (13) and the expressions for different components of electrostatic energy in  (8)-(10), we can rewrite the electrostatic energy $E_{\text{el}}=E_{H}+E_{\text{ext}}+E_{zz}$ as [72; 60]

[TABLE]

Finally, for given positions of nuclei, the reformulated governing equations for the Kohn-Sham DFT problem are:

[TABLE]

Though, the above equations  (14) and  (15) represent a pseudopotential treatment, we note that an all-electron treatment can be realized by setting $V^{J}_{\text{loc}}=\bar{V}^{J}_{\tilde{\delta}}$ and $V_{\text{nl}}=0$. We further remark that the equations  (14) and  (15) are equally valid for both periodic and non-periodic systems with appropriate boundary conditions. In a non-periodic setting, the simulation domain corresponds to a large enough domain, containing the compact support of the wavefunctions, with Dirichlet boundary conditions. In periodic calculations, it corresponds to a supercell with periodic boundary conditions.

For periodic systems, we now discuss the reduced Kohn-Sham equations on a periodic unit-cell. The Kohn-Sham eigenfunctions for infinite periodic crystals are given by the Bloch theorem  [62; 73], and the Bloch-periodic Kohn-Sham problem on an infinite crystal reduces to a periodic problem on a unit-cell. In numerical simulations involving such periodic calculations, it is computationally efficient to deal with unit-cells that are much smaller than the supercells, and the computation of electron-density, kinetic energy and the electrostatic interaction energy involving the non-local pseudopotentials has an additional integration over the Brillouin zone (BZ). Using Bloch theorem, the Kohn-Sham eigenfunction $\psi_{n,\bm{k}}(\bm{\textbf{x}})$ can be expressed as

[TABLE]

where $i=\sqrt{-1}$ and $u_{n,\bm{k}}(\bm{\textbf{x}})$ is a function that is periodic on the unit-cell while $\bm{k}$ denotes a reciprocal space point in the Brillouin zone. Using  (16), the computation of electron-density in equation  (2) is given by

[TABLE]

where $\fint\displaylimits_{BZ}$ denotes the volume average of the integral over the Brillouin zone corresponding to the periodic unit-cell $\Omega_{p}$ and $f_{n,\bm{k}}$ denotes the orbital occupancy function corresponding to $u_{n,\bm{k}}$. The contribution to the electrostatic interaction energy arising from the non-local pseudopotential using the Bloch theorem can be expressed as

[TABLE]

Using the separable form of the nonlocal pseudopotential operator, we have

[TABLE]

where the summation over $r$ runs on all lattice points in the periodic crystal, and $a$ runs on all the $N_{a}$ atoms in the unit-cell. We further note that $C_{lpm}^{a,n\bm{k}}$ in the above equation is given by

[TABLE]

Finally, using the local formulation of the extended electrostatic interaction energy discussed previously, the computation of the electronic ground-state, for a given position of atoms, in the context of unit-cell periodic DFT calculations is given by the following equations:

[TABLE]

In the above, periodic boundary conditions are imposed on $\Omega_{p}$ for the fields $u_{n}(\bm{\textbf{x}},\bm{k})$ and $\varphi(\bm{\textbf{x}})$. Furthermore, we remark that the equations  (21) involve an additional integration over the Brillouin zone (BZ) and is evaluated using numerical quadratures. These numerical quadratures replace the integration over BZ by a weighted sum over $\bm{k}$ points in the first BZ, and this sampling of the BZ is done using the Monkhorst-pack (MP) grid  [74] in DFT-FE. The symmetry operations of the underlying Bravais lattice is used to extract the reduced number of k-points sampling the irreducible Brillouin zone (IBZ). The origin of this MP grid could either coincide with the origin of the BZ or can be shifted by half of the grid spacing in order to maximize the benefit from crystal symmetry mediated BZ reduction  [62]. Finally, the set of equations  (21) is solved self-consistently, with the eigenvalue problem  (21a) being solved for $N$ lowest bands for every $\bm{k}$-point in IBZ. We remark that, as noted previously, periodic all-electron calculations can be realized by setting $V^{J}_{\text{loc}}=\bar{V}^{J}_{\tilde{\delta}}$ and $V_{\text{nl}}=0$.

2.2 Variational formulation

The Kohn-Sham governing equations discussed in the previous subsection are the Euler-Lagrange equations of a local variational Kohn-Sham problem that corresponds to the computation of the electronic ground-state free energy for a given position of atoms. The variational problem can be formulated in terms of wavefunctions, fractional occupancies and the electrostatic potentials as given by [60]:

[TABLE]

[TABLE]

[TABLE]

We note that $\bm{\Psi}=\{\psi_{1}(\bm{\textbf{x}}),\psi_{2}(\bm{\textbf{x}}),\psi_{3}(\bm{\textbf{x}}),\cdots,\psi_{N}(\bm{\textbf{x}})\}$, and $\bm{f}=\{f_{1},f_{2},f_{3}\cdots f_{N}\}$ denotes the vector of orbital occupancy factors, while $\mathcal{V}=\{V^{1},V^{2},\cdots,V^{N_{a}}\}$ denotes the vector containing the trial electrostatic potentials corresponding to all nuclear charges in the simulation domain. Here, $T_{\text{s}}(\bm{f},\bm{\Psi})$ denotes the kinetic energy of non-interacting electrons and $E_{\text{ent}}(\bm{f})$ denotes the electronic entropy contribution, and the corresponding expressions are given by

[TABLE]

The energy functional corresponding to electrostatic energy, $\mathcal{L}_{\text{el}}$, can be expressed in the local form as [60]

[TABLE]

where $\bar{V}^{J}_{\tilde{\delta}}$ denotes the electrostatic potential corresponding to the $J^{th}$ nuclear charge (see equation (13)), or analogously $\bar{\mathcal{V}}_{\delta}=\{\bar{V}^{1}_{\tilde{\delta}},\bar{V}^{2}_{\tilde{\delta}},\cdots,\bar{V}^{N_{a}}_{\tilde{\delta}}\}=arg\,\min_{\mathcal{V}\in(H^{1}({\mathbb{R}^{3}}))^{N_{a}}}\mathcal{L}_{\text{el}}(\bm{f},\bm{\Psi},\varphi,\mathcal{V};\bm{\textbf{R}})$. Further, we note that, $\mathcal{Y}$ in equation  (22) denotes a suitable function space that guarantees the existence of minimizers. We remark that numerical computations involve the use of bounded domains, which in non-periodic calculations correspond to a large enough domain containing the compact support of the wavefunctions, and, in periodic calculations, correspond to the super-cell333while the variational problem in equation (22) is presented for super-cells in the case of periodic calculations, it can be extended to periodic unit-cells using the Bloch Ansatz as discussed in  [60].. Denoting such an appropriate bounded domain by $\Omega$ subsequently, $\mathcal{Y}=H^{1}_{0}(\Omega)$ in the case of non-periodic problems, and $\mathcal{Y}=H^{1}_{per}(\Omega)$ in the case of periodic problems.

2.3 Configurational Forces

We employ configurational forces approach [60] to evaluate ionic forces and periodic unit-cell stresses in DFT-FE. These configurational forces correspond to the generalized variational force computed as the derivative of the Kohn-Sham energy functional  (22) with respect to the position of a material point x. This approach provides a unified framework to compute ionic forces as well as stress tensor for geometry optimization, and inherently accounts for the Pulay corrections owing to the variational nature of the formulation. For the sake of completeness, we present here the expressions for the configurational forces, in terms of the Kohn-Sham eigenfunctions, as the derivative of the energy functional  (22) with respect to x. We refer to [60] for the derivation of configurational forces, and a comprehensive discussion on this topic.

Let $\,\bm{\tau}^{\varepsilon}:{\mathbb{R}^{3}}\rightarrow{\mathbb{R}^{3}}$ represent the infinitesimal perturbation of the underlying space, mapping a material point x to a new point $\bm{\textbf{x}}^{\prime}$ such that ${\bm{\tau}}^{0}=\textbf{I}$. Further, let the generator of this mapping be $\bm{\Upsilon}=\frac{d}{d\varepsilon}\,\bm{\tau}^{\varepsilon}(\bm{\textbf{x}})\,|_{\varepsilon=0}$ such that $\,\bm{\tau}^{\varepsilon}$ is constrained to rigid body deformations in the compact support of the regularized nuclear charge distribution $\bm{b}(\bm{\textbf{x}})$ in order to preserve the integral constraint $\int\tilde{\delta}(\bm{\textbf{x}}-\bm{\textbf{R}}_{I})\,d\bm{\textbf{x}}=1$. Denoting $\mathcal{F}_{0}(\,\bm{\tau}^{\varepsilon})$ to be the ground-state free energy in the perturbed space, the configurational force is evaluated by computing the Gâteaux derivative of $\mathcal{F}_{0}(\,\bm{\tau}^{\varepsilon})$ given by

[TABLE]

where E and $\bm{\textbf{E}}^{\prime I}$ denote the Eshelby tensors whose expressions in terms of the solutions of the saddle point problem  (22) on the original space are provided below. If $\bar{\bm{\Psi}}=\{\bar{\psi}_{1},\bar{\psi}_{2},\bar{\psi}_{3}\cdots\bar{\psi}_{N}\}$ denote the Kohn-Sham eigenfunctions corresponding to the lowest $N$ eigenvalues $\bar{\bm{\epsilon}}=\{\bar{\epsilon}_{1},\bar{\epsilon}_{2},\cdots,\bar{\epsilon}_{N}\}$ with occupancies $\bar{\bm{f}}=\{\bar{f}_{1},\bar{f}_{2}\cdots\bar{f}_{N}\}$, and $\bar{\varphi}$ denotes the electrostatic potential (all solutions of the saddle point problem  (22)), the expressions for the Eshelby tensor E and $\text{F}^{\text{PSP}}$ in equation  (27) are given by

[TABLE]

where

[TABLE]

Further,

[TABLE]

where

[TABLE]

Though equation (2.3) represents the pseudopotential case, the expression for an all-electron case is realized by setting $E^{\text{nl}}=0$, $\text{F}^{\text{PSP}}=0$ and $V^{J}_{\text{loc}}(|\bm{\textbf{x}}-\bm{\textbf{R}}_{J}|)=\bar{V}^{J}_{\tilde{\delta}}(|\bm{\textbf{x}}-\bm{\textbf{R}}_{J}|)$. Finally, we remark that the configurational forces provide a unified expression for computing both the ionic forces as well as periodic unit-cell stress by using an appropriate choice of the generator $\bm{\Upsilon}(\bm{\textbf{x}})$. In particular, the force on any given atom is computed via equation  (27) by choosing the compact support of $\bm{\Upsilon}(\bm{\textbf{x}})$ to contain the atom of interest, and the periodic unit-cell stress tensor is evaluated by choosing $\Upsilon_{i}=C_{ij}x_{j}$ in  (27) as explained in  [60].

2.4 Discrete Kohn-Sham DFT equations

We introduce here the finite-element (FE) discretization of the Kohn-Sham DFT problem by representing various electronic fields in the FE basis, a piece-wise polynomial basis generated from the FE discretization [75]. In particular, we employ $C^{0}$ continuous Lagrange polynomial basis interpolated over Gauss-Lobatto-Legendre nodal points. The FE discretization of the Kohn-Sham DFT problem described here is along the lines of our prior work [34] and is briefly discussed here, highlighting the important differences in this work. We specifically note here that the real-space formulation of Kohn-Sham DFT as presented in equation  (22) results in a saddle point problem (min-max problem) in the electronic fields. Thus, it is possible that the electronic ground-state energy obtained from a single FE discretization of all the solution fields in the Kohn-Sham DFT problem can be non-variational. To address this, we seek to solve the electrostatic problem to a more stringent accuracy than the Kohn-Sham eigenvalue problem. To this end, we consider two FE triangulations for representing the wavefunctions and the electrostatic potentials, namely $\mathcal{T}^{h}$ and $\mathcal{T}^{h_{el}}$ with the characteristic mesh-sizes denoted by $h$ and $h_{el}$, respectively. We consider $\mathcal{T}^{h_{el}}$ to be a uniform subdivision of $\mathcal{T}^{h}$. Denoting the subspaces spanned by the FE basis corresponding to triangulations $\mathcal{T}^{h}$ and $\mathcal{T}^{h_{el}}$ to be $\mathbb{V}_{h}^{M}$ (with dimension $M$) and $\mathbb{V}_{h_{el}}^{M_{el}}$ (with dimension $M_{el}>M$), we note that $\mathbb{V}_{h}^{M}\subset\mathbb{V}_{h_{el}}^{M_{el}}$. Finally, the representation of the various fields in the Kohn-Sham problem (15)—the wavefunctions and the electrostatic potentials—in the FE basis is given by

[TABLE]

where $N^{h}_{j}:1\leq j\leq M$ denotes the FE basis spanning $\mathbb{V}_{h}^{M}$ and $N^{h_{el}}_{j}:1\leq j\leq M_{el}$ denotes the FE basis spanning $\mathbb{V}_{h_{el}}^{M_{el}}$ . We note that $\psi^{h}_{i}$, $\varphi^{h_{el}}$ and $\bar{V}^{J^{h_{el}}}_{\tilde{\delta}}$ denote the FE discretized fields, with $\psi^{j}_{i}$, $\varphi^{j}$ and $\bar{V}^{J^{j}}_{\tilde{\delta}}$ denoting the coefficients in the expansion of the $i^{th}$ discretized wavefunction and the electrostatic potentials, which also correspond to the nodal values of the respective fields at the $j^{th}$ node on the FE mesh.

The FE discretization of the Kohn-Sham eigenvalue problem (15) results in a generalized eigenvalue problem given by $\bm{\textbf{H}}\hat{\bm{\varPsi}}_{i}=\epsilon^{h}_{i}\bm{\textbf{M}}\hat{\bm{\varPsi}}_{i}$ where H denotes the discrete Hamiltonian matrix with matrix elements $\text{H}_{jk}$, M denotes the overlap matrix (or commonly referred to as the mass matrix in finite element literature) with matrix elements $\text{M}_{jk}$, and $\epsilon^{h}_{i}$ denotes the $i^{th}$ eigenvalue corresponding to the discrete eigenvector $\hat{\bm{\varPsi}}_{i}$. The expression for the discrete Hamiltonian matrix, $\text{H}_{jk}=\text{H}_{jk}^{\text{loc}}+\text{H}_{jk}^{\text{nl}}$, is given in terms of

[TABLE]

In the above, $V_{\text{eff,loc}}^{h}$ denotes the local part of the effective single-electron potential computed in the FE basis as the sum of discretized exchange-correlation potential $V_{\text{xc}}^{h}$, total electrostatic potential $\varphi^{h_{el}}(\bm{\textbf{x}})$ and the local pseudopotential term as follows:

[TABLE]

In the case of all-electron calculations, $V_{\text{eff,loc}}^{h}(\bm{\textbf{x}},\bm{\textbf{R}})=V_{\text{xc}}^{h}(\bm{\textbf{x}})+\varphi^{h_{el}}(\bm{\textbf{x}})$ and $\text{H}_{jk}^{\text{nl}}$ is zero. In the case of pseudopotential calculations, $H_{jk}^{\text{nl}}$ is given by

[TABLE]

Finally, the matrix elements of the overlap matrix M are given by $\text{M}_{jk}=\int_{\Omega}N^{h}_{j}(\bm{\textbf{x}})N^{h}_{k}(\bm{\textbf{x}})\,d\bm{\textbf{x}}$. We note that the matrices $\bm{\textbf{H}}^{\text{loc}}$ and M are sparse as the FE basis functions are local in real space and have a compact support (a finite region where the function is non-zero). Further, the vectors $C^{J}_{lpm,j}$ in $\bm{\textbf{H}}^{\text{nl}}$ are also sparse since the projectors $\chi^{J}_{lpm}(\bm{\textbf{x}},\bm{\textbf{R}}_{J})$ have a compact support, thus rendering a sparse structure to the discrete Hamiltonian H.

In order to explore efficient solution strategies, it is desirable to transform the generalized eigenvalue problem into a standard eigenvalue problem. Since the matrix M is positive definite symmetric, there exists a unique positive definite symmetric square root of M, and is denoted by $\bm{\textbf{M}}^{1/2}$. Hence, the following holds true:

[TABLE]

[TABLE]

We note that $\bm{\widetilde{\textbf{H}}}$ is a Hermitian matrix, and (32) represents a standard Hermitian eigenvalue problem. The actual eigenvectors are recovered by the transformation $\hat{\bm{\varPsi}}_{i}=\bm{\textbf{M}}^{-1/2}\widetilde{\bm{\varPsi}}_{i}$. Furthermore, we note that the matrix $\bm{\textbf{M}}^{-1/2}$ can be evaluated with modest computational cost by using a spectral FE basis in conjunction with the use of Gauss-Lobatto-Legendre (GLL) quadrature for the evaluation of integrals in the overlap matrix, that renders the overlap matrix diagonal [34]. This renders the matrix $\bm{\widetilde{\textbf{H}}}$ the same sparsity structure as the matrix H.

Finally, for the given positions of nuclei, the discrete Kohn-Sham eigenvalue problem along with the discretized Poisson equations for the electrostatic potentials ($\varphi^{h_{el}}$ and $\bar{V}^{J^{h_{el}}}_{\tilde{\delta}}$) are to be solved self-consistently, and are given by:

[TABLE]

We note that the nuclear charges in DFT-FE implementation are located on the nodes of the FE triangulation, and are treated as point charges. Thus, the nuclear charge distribution in the discrete setting $b^{h_{el}}(\bm{\textbf{x}},\bm{\textbf{R}})$ in equation (34b) is given by $b^{h_{el}}(\bm{\textbf{x}},\bm{\textbf{R}})=\sum_{I}b_{I}^{h_{el}}(|\bm{\textbf{x}}-\bm{\textbf{R}}_{I}|)$ with $b_{I}^{h_{el}}=-Z_{I}{\delta}(|\bm{\textbf{x}}-\bm{\textbf{R}}_{I}|)$ where $\delta(|\bm{\textbf{x}}-\bm{\textbf{R}}_{I}|)$ denotes the Dirac-delta distribution centered at the position of the atom $\bm{\textbf{R}}_{I}$. The boundary conditions used for the computation of the discrete potential field $\varphi^{h_{el}}(\bm{\textbf{x}})$ in equation  (34b) are either homogeneous Dirichlet boundary conditions or periodic boundary conditions depending on whether the problem is non-periodic or periodic. Further, the discrete self potential $\bar{V}^{J^{h_{el}}}_{\tilde{\delta}}$ associated with individual nuclear charge $J$ is solved using the discrete Poisson equation  (34c) subject to Dirichlet boundary conditions with prescribed Coulomb potential applied on a domain $\Omega_{J}$ enclosing the atom $J$. After obtaining the electronic ground-state from the solution of the discrete Kohn-Sham problem (equations (34)), we compute the discrete total ground-state energy $E^{h}$ in terms of the discrete solution fields ($\bar{\epsilon}^{h}_{i},\bar{\rho}^{h},\bar{\varphi}^{h_{el}},\bar{V}^{I^{h_{el}}}_{\tilde{\delta}}$) as follows:

[TABLE]

where,

[TABLE]

3 Numerical implementation

DFT-FE is built over the deal.II open-source finite-element library [76], and uses its underlying finite element constructs, adaptive mesh refinement architecture and efficient parallel vector objects. In this section, we discuss various aspects of numerical implementation of the Kohn-Sham DFT problem within the framework of spectral finite-element discretization in DFT-FE. We first begin with a discussion on the strategies implemented for adaptive mesh refinement in DFT-FE, followed by a detailed discussion on the various steps involved in the implementation of the Kohn-Sham self-consistent field iteration procedure.

3.1 Adaptive mesh refinement

One of the significant strengths of the FE basis is that it can accommodate adaptive spatial resolution, which in turn can be effectively exploited for the efficient solution of DFT calculations  [34; 77; 78; 27; 37] as well as the development of coarse-graining techniques that seamlessly bridge electronic structure calculations with continuum [79; 80; 81]. In DFT-FE, adaptive mesh refinement is carried out using octree-based hexahedral mesh generator based on the ‘Parallel AMR on Forests of Octrees’ (p4est) library  [82] via deal.II open-source finite-element library [76]. Spatial discretization via this octree refinement produces a non-conforming mesh (cf. Fig. 1), resulting in a potentially discontinuous function approximation. However, the continuity of a function across refinement interfaces is enforced by algebraic constraints that require the value of the function on the ‘hanging node’ to be consistent with the approximation along the neighboring element face or edge. We discuss the two adaptive mesh refinement strategies implemented in DFT-FE: (i) user-defined adaptive mesh refinement (UDAMR) procedure, involving user inputs for mesh sizes at different regions of the simulation domain, and (ii) automatic adaptive mesh refinement procedure (AAMR) guided by a priori error estimates. We note that the simulation domain encloses the given materials system for a non-periodic problem, while for a periodic problem, the simulation domain is a periodic unit-cell.

3.1.1 User defined adaptive mesh refinement (UDAMR)

The user defined adaptive mesh in DFT-FE is constructed to have a refined mesh around the atom up to a certain radius $r_{\text{atom}}$ and coarsens away. This is accomplished by discretizing the given simulation domain using an initial coarse uniform mesh with mesh size $h_{\text{base}}$. Subsequently, the FE cells whose centroids are within a distance of $r_{\text{atom}}$ units from each of the atomic positions are marked for refinement. These marked cells are refined until a target mesh size $h_{\text{atom}}$ is achieved. Hence, the mesh size parameters $h_{\text{base}}$, $h_{\text{atom}}$ and the parameter $r_{\text{atom}}$ form the user defined input to generate the adaptive mesh. Though these parameters are adequate for generating adaptive meshes in the case of pseudopotential DFT calculations involving smooth solution fields, all-electron DFT calculations require much finer meshes in the close vicinity of the nuclear position due to the highly oscillatory nature of the wavefunctions near the nuclei. Hence, we introduce another mesh parameter $h_{\text{fine}}$ which prescribes the mesh-size of the elements that share the vertex at the nuclear position. A user defined adaptive mesh depicting the above parameters in the case of SiF4 molecule for all-electron DFT calculation is shown in Fig. 1 .

We note that the mesh parameters $h_{\text{base}}$, $h_{\text{atom}}$, $h_{\text{fine}}$ for a given materials system can be estimated from the optimal mesh size distribution $h(\bm{\textbf{x}})$ obtained by minimizing the discretization error in the ground-state energy for a fixed number of FE cells $N_{E}$. To this end, we recall the procedure described in  [34] to determine this mesh-size distribution. Let $E$ be the ground-state energy for the continuous problem  (15) with $\bar{\bm{\Psi}}=\{\bar{\psi}_{1}\,,\,\bar{\psi}_{2}\,\cdots\,\bar{\psi}_{N}\}$ representing the electronic wavefunction solutions of the continuous problem  (15), and, further, let $E_{h}$ denote the discretized ground-state energy with $\bar{\bm{\Psi}}^{h}=\{\bar{\psi}_{1}^{h}\,,\,\bar{\psi}_{2}^{h}\,\cdots\,\bar{\psi}_{N}^{h}\}$ representing the electronic wavefunction solutions of the discrete problem (34). As discussed in the previous section, we note that the triangulation $\mathcal{T}^{h_{e}}$ employed in the solution of the electrostatics problems  (34b) and  (34c) is more resolved than the triangulation $\mathcal{T}^{h}$ used for solving Kohn-Sham eigenvalue problem  (34a), and is usually chosen to be the uniform subdivision of the triangulation $\mathcal{T}^{h}$. Thus, the dominant discretization error will correspond to the error in the wavefunctions, and as derived in equation (47) of  [34], the discretization error $|E-E_{h}|$ (retaining only the leading order terms) is given by

[TABLE]

where $e$ denotes an element in the FE discretization, with mesh size $h_{e}$ covering the domain $\Omega_{e}$, and $|\,.\,|_{k+1,\Omega_{e}}$ denotes the $(k+1)$ semi-norm over $\Omega_{e}$. Using the definition of semi-norms, equation (36) can be written as

[TABLE]

where $h(\bm{\textbf{x}})$ denotes the element size distribution function defining the target element size at point x in the simulation domain. The mesh distribution $h(\bm{\textbf{x}})$ can be estimated by minimizing the approximation error in energy in  (37) subject to fixed number of elements, which is given by

[TABLE]

The solution to this variational problem is given by

[TABLE]

where the constant $A$ is computed from the constraint that the total number of elements is $N_{E}$. We note that the mesh size distribution $h(\bm{\textbf{x}})$ in equation  (39) involves the knowledge of $\bar{\psi}_{i}(\bm{\textbf{x}})$, which is not known a priori. However, from a practical standpoint, single atom DFT wavefunctions constructed numerically from the solution of 1D radial Kohn-Sham DFT problem can be used. Thus, the mesh parameters $h_{\text{base}}$, $h_{\text{atom}}$, $h_{\text{fine}}$ can be estimated from $h(\bm{\textbf{x}})$, which provides a systematic procedure to prescribe these quantities.

We note that the adaptive mesh refinement procedure via the octrees employed in DFT-FE produces a non-conforming mesh such that the ratio of edge lengths of two neighboring cells is at most 2:1. Due to this constraint the target mesh sizes of $h_{\text{atom}}$ and $h_{\text{fine}}$ can only be approximately realized, and may also lead to more degrees of freedom than required. This increase in degrees of freedom can be more pronounced in the case of all-electron DFT calculations. Thus, DFT-FE also provides an automatic adaptive mesh generation strategy, guided by local error indicators involving the FE discretized solution fields, and is discussed subsequently.

3.1.2 Automatic adaptive mesh refinement (AAMR)

A number of recent works [83; 84; 85; 86; 37; 87] have been devoted to adaptive mesh refinement strategies employing local error indicators expressed in terms of the solution fields of the discrete Kohn-Sham problem. Local error indicators relying on eigenvalue problem residuals as well as jump in the derivative of solution fields across the face of FE cells (Kelly error indicators) have been employed in many of the recent works  [83; 84; 86; 37]. Further, error indicators based on $H^{1}$ semi-norms of the wavefunctions [85] and those based on coarsening mesh approaches [87] have also been employed. Most of the above methods start with an initial coarse triangulation on which the discrete Kohn-Sham problem is solved, and a local error indicator in terms of the discrete solution fields is subsequently employed to mark the cells to be refined (a posteriori mesh adaption). This procedure is usually repeated till convergence, and the process generates a sequence of increasingly refined adaptive FE approximations.

We note that many of the aforementioned a posteriori adaptive mesh refinement strategies require the solution of the Kohn-Sham problem during the course of the adaptive refinement procedure, which can be very expensive to employ for large scale calculations 444computational complexity of the Kohn-Sham problem is cubic-scaling with number of electrons, thus making the repeated solution of a large-scale problem very expensive during the course of adaptive mesh refinement procedure.. To this end, we present an efficient strategy to construct an adaptive mesh a priori, before beginning the SCF procedure, by making use of numerically computed single-atom Kohn-Sham DFT wavefunctions. The approach is based on a local error indicator obtained from an error estimate on the energy involving the discrete wavefunctions, in contrast to the error estimate in equation (36) that involves the continuous wavefunctions. We first derive this energy error estimate following the mesh adaption ideas in  [88], and subsequently present the algorithm implemented in DFT-FE for automatic adaptive mesh refinement.

Let $\widetilde{\psi}^{h}(\bm{\textbf{x}})$ be the $(k-1)^{\text{th}}$ interpolant of $\bar{\psi}^{h}(\bm{\textbf{x}})$ with $k\geq 2$ denoting the FE interpolating polynomial. As opposed to the approximation error $|E-E^{h}|$, we will work with $|E-\widetilde{E}^{h}|$ where $\widetilde{E}^{h}$ denotes the discrete Kohn-Sham ground-state energy obtained using the $(k-1)^{\text{th}}$ interpolants of the FE solution of the Kohn-Sham problem (equation (34)). To this end, following along similar lines as the derivation of equation  (36), the dominant term in $|E-\widetilde{E}^{h}|$ can be derived to be

[TABLE]

We recall, in the above expression, $|\,.\,|_{1,\Omega}$ denotes the $H^{1}$ semi-norm over $\Omega$. Using triangle inequality, we have

[TABLE]

We note that, asymptotically, as $h\rightarrow 0$, the first term, ${|\bar{\psi}_{i}-\bar{\psi}^{h}_{i}|}_{1,\Omega}$ , in the above inequality (41) is $\order{h^{k}}$  [34], while the second term, ${|\bar{\psi}^{h}_{i}-\widetilde{\psi}^{h}_{i}|}_{1,\Omega}$ , is of order $\order{h^{k-1}}$. Hence, the error is dominated by ${|\bar{\psi}^{h}_{i}-\widetilde{\psi}^{h}_{i}|}_{1,\Omega}$ , which in turn can be bounded in terms of the FE mesh size as [89]

[TABLE]

where $e$ denotes an element in the FE discretization, with mesh size $h_{e}$ covering the domain $\Omega_{e}$. Hence, the energy error estimate $|E-\widetilde{E}^{h}|$ in equation  (40) can be written as

[TABLE]

The energy error estimate  (43) motivates $r_{e}=(h_{e})^{2(k-1)}\sum_{i}{|D^{k}\bar{\psi}_{i}^{h}(\bm{\textbf{x}})|}^{2}$ as a useful local error indicator to be used for the automatic adaptive mesh refinement procedure. This error indicator is ideally suited for a-posteriori mesh adaption, i.e., use the solution from the current mesh to conduct mesh adaption. However, such a procedure will be impractical for large-scale calculations, where the Kohn-Sham problem has to be solved for every iterate of the refinement procedure. Thus, in DFT-FE, we adopt the strategy to use this error indicator to conduct mesh refinement a-priori. To this end, we use single-atom wavefunctions interpolated onto the FE mesh (denoted as $\psi_{i}^{h}(\bm{\textbf{x}})$) as good approximations to $\bar{\psi}_{i}^{h}(\bm{\textbf{x}})$ for the purpose of mesh-refinement, and use these in the computation of the local error indicator.

Thus, the automatic adaptive mesh refinement procedure in DFT-FE starts with an initial triangulation $\mathcal{T}^{h^{(0)}}$, and a sequence of nested triangulations $\mathcal{T}^{h^{(n)}}$ are generated using the following iterative procedure (Algorithm 1): Interpolate single-atom DFT wavefunctions $\rightarrow$ Estimate local error $\rightarrow$ Mark for refinement $\rightarrow$ Refine $\rightarrow$ Check convergence. In particular, the triangulation $\mathcal{T}^{h^{(n+1)}}$ is constructed from $\mathcal{T}^{h^{(n)}}$ by interpolating the single-atom DFT wavefunction data onto the current mesh $\mathcal{T}^{h^{(n)}}$. The local error indicator ${r_{e}}^{(n)}$ is computed for each FE cell, and the cells are ranked in the order of decreasing error. A predefined fraction ($\beta$) of cells are marked for refinement, by selecting those at the top of the ordered cells (corresponding to the highest local error), and the refinement procedure is carried out to generate the triangulation $\mathcal{T}^{h^{(n+1)}}$. As the AAMR is executed as an a-priori adaption scheme, a key aspect of AAMR is to devise a stopping criterion for the refinement algorithm. To this end, as the sequence of nested triangulations get generated in AAMR, we examine the convergence of kinetic energy term, i.e, $\sum_{i}\int|\nabla\psi_{i}^{h^{(n)}}(\bm{\textbf{x}})|^{2}\,\,d\bm{\textbf{x}}$. This is motivated from  (40), as the leading order error in the ground-state energy results from the error in the kinetic energy. The refinement algorithm is terminated when the kinetic energy is converged to within a prescribed tolerance.

In order to compare AAMR with UDAMR, we conduct a comparative study between these mesh adaption schemes. To this end, we choose representative benchmark examples involving non-periodic norm conserving (ONCV) pseudopotential, all-electron DFT calculations on SiF4 molecule, and a periodic all-electron DFT calculation on Si diamond unit-cell. The mesh attributes, degrees of freedom per atom, and the discretization errors in ground-state energies, atomic forces and periodic unit-cell hydrostatic stresses are tabulated for the above benchmark systems (cf. Tables 1–3). We note that the reference data for the ground-state energy per atom ($E_{0}$), force vector ($\bm{f}_{0}$) and hydrostatic stress ($\sigma_{0}$) to measure the discretization errors are obtained using highly refined FE calculations. These results indicate that AAMR procedure as described in Algorithm 1 provides FE meshes with significantly lesser degrees of freedom than the UDAMR procedure for a range of discretization errors. In the case of pseudopotential DFT calculations on SiF4 molecule, we observe that AAMR scheme resulted in a FE mesh with $3$ times lesser degrees of freedom than the UDAMR scheme for discretization errors of the $\order{10^{-5}}$ in ground-state energies and forces. While in the case of all-electron DFT calculations on the same benchmark system, AAMR scheme resulted in a FE mesh with $20$ times lesser degrees of freedom than the UDAMR scheme for discretization errors of the $\order{10^{-4}}$ in ground-state energies and forces. Fig. 2 shows the slices of a 3D FE meshes obtained using AAMR scheme at discretization errors of around $10^{-4}$ Ha/atom and $10^{-5}$ Ha/atom in the case of SiF4 all-electron DFT calculations. Further, in the case of all-electron DFT periodic calculations on Si diamond unit-cell, AAMR scheme resulted in a FE mesh with $3$ times lesser degrees of freedom for discretization errors of the $\order{10^{-4}}$ in ground-state energies.

We also demonstrate the systematic convergence of energies and forces obtained from the sequence of FE meshes generated using AAMR. To this end, we plot in Fig.  3(a) and Fig. 3(b) the absolute discretization errors in ground-state energy and force vector (on a specific atom) for the sequence of increasingly refined meshes obtained from AAMR. Fig. 3 demonstrates that absolute discretization errors as low as $\order{10^{-5}}$ Ha/atom in the ground-state energy and $\order{10^{-5}}$ Ha/Bohr in the forces can be attained using AAMR, which are significantly better than chemical accuracy. We note that the discretization errors stagnate for increasingly finer meshes in Fig. 3 and this behavior is a consequence of using a priori error indicator involving single atom wavefunctions in AAMR.

3.2 SCF Algorithm

The discrete nonlinear Hermitian eigenvalue problem is solved self-consistently along with Poisson equations (see equation  (34)) to compute the Kohn-Sham ground-state solution. Algorithm 2 lists all the steps in the SCF procedure followed in DFT-FE. We use adaptive higher order spectral finite-elements in conjunction with computationally efficient and scalable Chebyshev filtered subspace iteration technique (ChFSI) [57; 34] to evaluate the occupied eigenspace of the discrete Kohn-Sham Hamiltonian. We further employ Anderson and Broyden schemes [90; 91] for electron-density mixing, and the finite-temperature Fermi-Dirac smearing [19] to avoid the charge sloshing associated with metallic systems.

The ChFSI procedure in Algorithm 2 involves the Chebyshev filtering (CF), orthonormalization (CholGS), and the Rayleigh-Ritz procedure (RR). We note that CF scales quadratically with number of atoms, while CholGS and RR scale cubically with number of atoms. Thus, for small to medium scale system sizes CF is the dominant computational cost, while for larger system sizes the computational cost of CholGS and RR dominates. To this end, the numerical implementation in DFT-FE focuses on reducing the prefactor and improving scalability of the ChFSI procedure by exploiting efficient methods and cache-friendly data-structures like FE cell level matrix-matrix multiplications, mixed precision strategies and spectrum splitting approach, as will be discussed subsequently. Furthermore, the electrostatic potentials are computed by solving a Poisson problem, which employs a matrix-free framework of the deal.II finite-element library [92; 76] in conjunction with a Jacobi preconditioned conjugate gradient solver. We note that the above matrix-free framework computes the matrix-vector product of the FE operator on the fly without ever storing it as a sparse matrix. Such on the fly computations benefit from significantly lower memory access costs and have been demonstrated to outperform global sparse-matrix based methods on modern computing architectures [92].

3.3 Chebyshev filtering

In DFT-FE, Chebyshev polynomial filtering technique [59] is used to adaptively approximate the wanted eigenspace (the lowest $N$ occupied eigenfunctions) of the FE discretized Hamiltonian ${\bf\widetilde{H}}$ [34]. In practice, $N$ is typically chosen as $N_{e}/2+b$ to allow for finite-temperature Fermi-Dirac smearing, where $b$ is usually $\rm{(5-10)\%}$ of $N_{e}/2$. In a given SCF iteration step, a scaled Hamiltonian ${\bf\bar{H}}$ is obtained by scaling and shifting ${\bf\widetilde{H}}$ such that the unwanted spectrum of ${\bf\widetilde{H}}$ is mapped on to $[-1,1]$, and the wanted spectrum is mapped on to $\left(-\infty,-1\right)$ to exploit the fast growth property of Chebyshev polynomials in this region. Subsequently, the action of a degree $m$ Chebyshev polynomial filter, $T_{m}({\bf\bar{H}})$, on the input subspace, $\bm{\widetilde{\Psi}}$, is computed recursively as

[TABLE]

We use an adaptive filtering strategy in which multiple sweeps of ChFSI procedure are performed till the residual norm of the eigenpair closest to the Fermi energy reaches below a specified tolerance $\delta$, chosen to between $1\times 10^{-2}-5\times 10^{-2}$. Our numerical experiments in the case of pseudopotential electronic ground-state calculations show that while multiple calls to ChFSI are triggered in the first few SCF iterations, there is an overall reduction in the number of ChFSI calls (due to reduced number of SCF iterations) when employing the adaptive filtering strategy in comparison to employing a single sweep in all SCF iterations. We remark that despite using the adaptive filtering strategy, for atomic relaxations or molecular dynamics simulations, multiple Chebyshev filtering calls are typically not triggered as the wavefunctions from the previous electronic ground-state calculation are reused as a starting guess. We note that the choice of the Chebyshev polynomial degree $m$ in equation (44) is based on the upper bound of the spectrum of ${\bf\widetilde{H}}$, which is governed by the smallest mesh size employed in the finite element discretization. A Chebyshev polynomial degree between 20–50 is typically used in DFT-FE for pseudopotential calculations, whereas significantly higher Chebyshev polynomial degrees ($\sim 500-1000$) are required for all-electron calculations.

3.3.1 Practical implementation aspects of Chebyshev filtering

The computational complexity of Chebyshev filtering scales as $\mathcal{O}(MN)$, where $M$ is the size of the discretized Hamiltonian ${\bf\widetilde{H}}$ and $N$ is the number of occupied states. Since Chebyshev filtering is the dominant computational cost in DFT-FE for small to medium sized systems (up to 20,000 electrons), we optimize the core kernel in the Chebyshev filtering procedure, which involves the computation of ${\bf\bar{H}}{\bf X}$ in equation (44), with ${\bf X}$ denoting a trial subspace in the course of the Chebyshev recursive iteration. To this end, we first explicitly compute and store the FE Hamiltonian matrices (cell level Hamiltonian matrices), and subsequently extract the cell level wavefunction matrices from the global wavefunction vectors ${\bf X}$. We then employ BLAS Xgemm routines to compute the matrix-matrix products involving cell Hamiltonian and wavefunction matrices, and assemble them to get the global wavefunction vectors. We note that global FE sparse matrix approaches, particularly when dealing with large number of wavefunction vectors, are more memory-bandwidth limited555The cell level matrix approach is similar in spirit to matrix-free based approaches, which have been demonstrated to have lower memory access costs than global FE sparse-matrix based methods [92]. and incur a higher communication cost666The global FE sparse matrix framework in deal.II library currently does not take advantage of performing MPI communication of multiple vectors in a single communication call. than the cell level matrix approach employed above.

In case of large problems with many thousands of wavefunction vectors, the peak memory during Chebyshev filtering can be quite high if implemented naively by filtering all the wavefunction vectors simultaneously, as multiple temporary memories of size ${\bf X}$ are needed in the course of the Chebyshev recursive iteration. Hence, to reduce the peak memory, we use a blocked approach by filtering blocks of wavefunction vectors, ${\bf X}_{b}$ with block size denoted by $B^{f}$, based on the rationale that Chebyshev filtering can be performed on each wavefunction vector independently. Further, the blocked approach also allows us to take advantage of batched Xgemm 777Batched operations are efficient for performing many small matrix-matrix multiplications concurrently on multiple threads. Currently such routines are available in vendor optimized BLAS libraries such as Intel MKL. routines in ${\bf\widetilde{H}}{\bf X}_{b}$ to perform the aforementioned cell level matrix-matrix products concurrently on multiple threads, which we found to be faster than using multiple threads on standard Xgemm calls involving very skewed matrix dimensions when blocked approach is not used. Additionally, we use a single contiguous memory block to store the global wavefunction vectors as well as the block wavefunction vectors, where the data layout is such that for each degree of freedom the corresponding wavefunction values are stored contiguously. This leads to more cache-friendly data access while copying the data between the global wavefunction vectors and the cell wavefunction matrices. Furthermore, we exploit the fact that all wavefunction vectors have identical communication pattern to minimize the total number of MPI point-to-point communication calls in ${\bf\widetilde{H}}{\bf X}_{b}$, which reduces the network latency.

The optimal value of the Chebyshev filtering block size, $B^{f}$, depends on two competing factors—very small sizes lead to higher memory access overheads and communication latency, whereas very large sizes increase peak memory and reduce the efficiency of batched Xgemm routines. Based on numerical experiments, we find the optimal range of $B^{f}$ to be between 300–400, which is set as the default in DFT-FE.

3.4 Cholesky factorization based Gram-Schmidt orthonormalization

ChFSI involves orthonormalization procedure after the Chebyshev filtering step to prevent the ill-conditioning of the filtered vectors in the course of the subspace iteration procedure. This procedure scales cubically with number of electrons and becomes one of the dominant computational costs in large-scale problems (greater than 20,000 electrons). To this end, we employ Cholesky factorization based Gram-Schmidt (CholGS) orthonormalization technique in DFT-FE. This is shown to be more efficient and scalable [46; 93] than the commonly used classical Gram-Schmidt procedure. Algorithm 3 shows the steps involved in the CholGS procedure. The $\mathcal{O}(N^{2})$ dot products involved in classical Gram-Schmidt are replaced by more cache-friendly matrix-matrix multiplications in CholGS (steps 1 and 4). Furthermore, the single communication call involved in the computation of overlap matrix ${\bf S}$ in CholGS has a much lower communication latency in comparison to $\mathcal{O}(N^{2})$ communication calls in classical Gram-Schmidt.

3.4.1 Parallel implementation aspects of CholGS in Algorithm 3

Computation of overlap matrix

We first note that $\bm{\widetilde{\Psi}}$ is stored in parallel as a $M_{\rm loc}\times N$ matrix, where $M_{\rm loc}$ is the number of FE nodes owned locally by a given MPI task. Accordingly, a straightforward approach to compute the overlap matrix ${\bf S}$ in step 1 involves the evaluation of local contributions of $\bm{\widetilde{\Psi}^{\dagger}}\bm{\widetilde{\Psi}}$ (a $N\times N$ matrix) on each MPI task, and then accumulating the local contributions to ${\bf S}$ using the MPI_Allreduce collective routine. However, this approach requires memory corresponding to a $N\times N$ matrix on each MPI task, and hence is not practically applicable for large-scale problems $\sim(N>20,000)$. To avoid this large memory footprint in both storage of ${\bf S}$ as well as computation of the local contributions, we use the popular 2D cyclic block grid distribution of ScaLAPACK library [94] to distribute the memory of ${\bf S}$, and use a blocked approach to compute the local contributions of $\bm{\widetilde{\Psi}^{\dagger}}\bm{\widetilde{\Psi}}$ to ${\bf S}$. Further in the blocked approach, we also exploit the Hermiticity of ${\bf S}$, by computing only the lower triangular portion of ${\bf S}$. Fig. 4 shows the schematic of the blocked approach with block size $B^{v}$, where $[i,N]\times[i,i+B^{v}]$ sub-matrices of ${\bf S}$ are computed successively one after another. Computation of each sub-matrix first involves computation of the local contribution in each MPI task by performing matrix-matrix multiplication between $[i,N]\times M_{\rm loc}$ block of $\bm{\widetilde{\Psi}^{\dagger}}$ and $M_{\rm loc}\times[i,i+B^{v}]$ of $\bm{\widetilde{\Psi}}$ using BLAS Xgemm routine, followed by accumulation of the local contributions using the MPI_Allreduce collective. Subsequently, the corresponding sub-matrix entries of the ScaLAPACK parallelized ${\bf S}$ are filled. Overall, the above blocked approach combined with ScaLAPACK parallelization of ${\bf S}$ provides both memory optimization and efficiency improvements.

Computation of inverse of Cholesky factor

Cholesky factorization of ${\bf S}$ in step 2 and inversion of the Cholesky factor ${\bf L}$ in step 3 are performed using ScaLAPACK routines pXpotrf and pXtrtri, respectively. Based on the numerical experiments conducted on a large benchmark systems, we find that the steps 2 and 3 are a minor cost compared to other steps in CholGS. For instance, the cost of steps 2 and 3 combined contributed to about 7% of the total wall time for CholGS for a system containing 61,502 electrons (see Fig. 6).

Construction of orthonormal vectors

Similar to step 1, computation of the orthonormalized basis ${\bm{{\widetilde{\Psi}}}^{{\rm o}}}$ in step 4 also has a large memory footprint when performed simply as a matrix-matrix multiplication between the local portion ($M_{\rm loc}\times N$ matrix) of the parallel distributed ${\bm{\widetilde{\Psi}}}$ and the full ${\bf{L^{-1}}^{\dagger}}$ ($N\times N$ matrix) on every MPI task. For large-scale problems this leads to a high peak memory due to storage of the full ${\bf{L^{-1}}^{\dagger}}$ on every MPI task, and also to store the computed ${\bm{{\widetilde{\Psi}}}^{{\rm o}}}$, which requires the same memory size as ${\bm{\widetilde{\Psi}}}$. Hence, we compute ${\bm{{\widetilde{\Psi}}}^{{\rm o}}}$ using two blocked levels to address both of these memory issues, as shown schematically in Fig. 5. First, we employ an outer blocked level over $M_{\rm loc}$ with block size $B^{d}$, which allows reuse of the memory of $\bm{{\widetilde{\Psi}}}$ to store ${\bm{{\widetilde{\Psi}}}^{{\rm o}}}$. In particular, we compute $[i,i+B^{d}]\times[1,N]$ sub-matrices of ${\bm{{\widetilde{\Psi}}}^{{\rm o}}}$ one after the other and copy the orthonormalized sub-matrices back on to $\bm{{\widetilde{\Psi}}}$, thereby requiring only an additional $B^{d}\times N$ memory. Secondly, we employ an inner blocked level where each $[i,i+B^{d}]\times[1,N]$ sub-matrix of ${\bm{{\widetilde{\Psi}}}^{{\rm o}}}$ is further divided into $[i,i+B^{d}]\times[j,j+B^{v}]$ sub-matrices and successively computed. Similar to the blocked approach used in step 1, this inner blocked level removes the requirement to store the full ${\bf{L^{-1}}^{\dagger}}$ while also exploiting the triangular matrix property of ${\bf{L^{-1}}^{\dagger}}$. Each $[i,i+B^{d}]\times[j,j+B^{v}]$ sub-matrix in the inner blocked level is computed by performing a matrix-matrix multiplication between a $[i,i+B^{d}]\times[1,j+B^{v}]$ sub-matrix of ${\bm{\widetilde{\Psi}}}$ and a $[1,j+B^{v}]\times[j,j+B^{v}]$ sub-matrix of ${\bf{L^{-1}}^{\dagger}}$. We note that ${\bf{L^{-1}}^{\dagger}}$ is stored in a ScaLAPACK parallel format after the end of step 3. Thus to obtain the $[1,j+B^{v}]\times[j,j+B^{v}]$ sub-matrix of ${\bf{L^{-1}}^{\dagger}}$ in each MPI task, we first use the local portion of the parallel ${\bf{L^{-1}}^{\dagger}}$ to fill the corresponding entries in the sub-matrix and the rest as zeros, and subsequently use the MPI_Allreduce collective to gather and communicate the filled sub-matrix to all MPI tasks.

Remarks on block sizes

We now discuss few considerations regarding the choice of optimal values for the block sizes $B^{v}$ (used above in steps 1 and 4) and $B^{d}$ (used above in step 4). Too small values of $B^{v}$ will lead to computational overheads in the Xgemm calls due to the highly skewed matrix dimensions which are not cache-friendly, and, further, the total number of MPI collective communication calls will increase leading to higher communication latency. On the other hand too large values of $B^{v}$ will deprecate the efficiency benefit of exploiting the Hermiticity of ${\bf S}$ in step 1 and the triangular matrix nature of ${\bf{L^{-1}}^{\dagger}}$ in step 4. Based on numerical experiments, we find that value of $B^{v}$ between 350–500 is optimal. Similarly, the choice of $B^{d}$ is based on two competing factors— too small values of $B^{d}$ incur higher computational and communication overheads due to repeated access of ${\bf{L^{-1}}^{\dagger}}$ for every outer level block computation, whereas larger values increase the peak memory required in step 4. We find that $B^{d}$ values between 2000–3000 have very negligible overhead costs while still providing memory efficiency when $M_{\rm loc}$ is much larger than $B^{d}$.

3.4.2 Mixed precision approaches in CholGS

To further reduce the prefactor of the CholGS algorithm, we make use of mixed precision arithmetic in steps 1 and 4 of Algorithm 3, which are the dominant costs in the CholGS algorithm. Mixed precision approaches for orthonormalization in the context of electronic-structure calculations have been explored previously by [95]. We first develop a mixed precision approach for step 1, where the computation of the overlap matrix, S can be split into computation of the diagonal and the off-diagonal parts:

[TABLE]

where ${\bf S_{d}}$ is a matrix containing the diagonal entries of ${\bf S}$. We take advantage of the fact that ${\bf S_{od}}\to\bf{0}$ as the SCF approaches convergence and hence compute ${\bf S_{od}}$ using single precision BLAS Xgemm routines, while the computation of diagonal entries of ${\bf S_{d}}$ is performed using double precision BLAS routines at negligible computational cost. Similarly, step 4 can be split into

[TABLE]

where ${\bf{L^{-1}_{d}}^{\dagger}}$ is a matrix containing the diagonal entries of ${\bf{L^{-1}}^{\dagger}}$. Taking advantage of the fact that ${\bf{L^{-1}_{od}}^{\dagger}}\to\bf{0}$ as the SCF approaches convergence, we compute ${\bm{\widetilde{\Psi}}}{\bf{L^{-1}_{od}}^{\dagger}}$ using single precision BLAS Xgemm routines, while the computation of ${\bm{\widetilde{\Psi}}}{\bf{L^{-1}_{d}}^{\dagger}}$ is performed as a double precision scaling operation at negligible computational cost. We remark that, in addition to the reduction of computational costs, the use of mixed precision also reduces the communication costs in steps 1 and 4 as the MPI_Allreduce collectives employed in these steps communicate the relevant single precision data with half the MPI message size (bytes), in comparison to their double precision counterparts.

The computational cost reduction in steps 1 and 4 of the mixed precision approach is demonstrated in Fig. 6 for large-scale benchmark problems involving 39,900 and 61,502 electrons. We find this approach to be around 2 times faster in comparison to double precision approach. Furthermore, we also examine the accuracy and robustness of the mixed precision algorithm in the overall SCF convergence in Section 3.5.2, and is discussed in detail subsequently.

3.5 Rayleigh-Ritz procedure and electron-density computation

Rayleigh-Ritz (RR) procedure in ChFSI involves the following steps: i) computation of the projected Hamiltonian, ${\bf\hat{H}}=\bm{{\widetilde{\Psi}}^{{\rm o}^{\dagger}}}{\bf\widetilde{H}}\bm{{\widetilde{\Psi}}^{{\rm o}}}$ into the space spanned by the orthonormalized wavefunctions $\bm{{\widetilde{\Psi}}^{{\rm o}}}$, ii) diagonalization of ${\bf\hat{H}}$: ${\bf\hat{H}}\bm{{\rm Q}}=\bm{{\rm Q}}\bm{{\rm D}}$, where $\bm{{\rm D}}$ contains all the eigenvalues of ${\bf\hat{H}}$ in ascending order and $\bm{{\rm Q}}$ contains the corresponding eigenvectors, iii) subspace rotation of $\bm{{\widetilde{\Psi}}^{{\rm o}}}$: $\bm{{\widetilde{\Psi}}}^{\bf R}=\bm{{\widetilde{\Psi}}^{{\rm o}}}\bm{{\rm Q}}$. Subsequently, the output electron-density at a point ${\bf x}$ belonging to a FE cell $e$ is computed as

[TABLE]

where ${\left\{N^{e}_{1}({\bf x}),\,N^{e}_{2}({\bf x}),\,\cdots\,N^{e}_{M_{e}}({\bf x})\right\}}$ denote the FE basis functions associated with the given cell ($M_{e}$ denoting the number of nodes in the cell), and ${\left\{{\psi}_{i}^{e,1},\,{\psi}_{i}^{e,2},\,\cdots\,{\psi}_{i}^{e,M_{e}}\right\}}$ denote the corresponding nodal values of the $i^{\rm th}$ wavefunction, $\bm{{\psi}}^{h}_{i}({\bf x})$ in FE cell $e$. Using the subspace rotated wavefunctions $\bm{{\widetilde{\Psi}}}^{\bf R}$, Equation 47 can be re-written as

[TABLE]

where

[TABLE]

and the matrix $\bm{{\hat{\Psi}}^{\bf R}_{e}}$ contains the FE cell column vectors extracted from $\bm{\hat{\Psi}}^{\bf R}$ which is given by

[TABLE]

In the above, the computational complexity of steps i), ii) and iii) of the Rayleigh-Ritz procedure scales as $\mathcal{O}(MN^{2})$, $\mathcal{O}(N^{3})$, and $\mathcal{O}(MN^{2})$, respectively, while the electron-density computation scales as $\mathcal{O}(MN)$. Rayleigh-Ritz procedure is, thus, one of the significant bottlenecks for large-scale problems. To this end, we employ two strategies in DFT-FE: spectrum-splitting and mixed precision, to reduce the prefactor of Rayleigh-Ritz procedure, as discussed below.

3.5.1 Spectrum-splitting in RR

The key idea behind spectrum-splitting is that the eigenvalues and eigenvectors of the projected Hamiltonian ${\bf\hat{H}}$ with orbital occupancy function $f_{i}=1$ are not explicitly necessary for the computation of the electron-density in equation (48). This can be exploited to achieve significant computational savings when most of the Kohn-Sham states are fully occupied as is the case for typically used Fermi-Dirac smearing temperatures $\sim 500\,{\rm K}$. Such methods have been developed in previous works in the context of both pseudopotential [41; 96] and all-electron DFT [97] calculations, which we have adapted in DFT-FE. Furthermore, we additionally take advantage of spectrum-splitting to develop a mixed precision technique to reduce the computational cost of the projected Hamiltonian computation. We discuss below the implementation of the spectrum-splitting algorithm in DFT-FE.

Let $N_{\rm oc}$ denote the number of Kohn-Sham eigenstates with full occupancies ($f_{i}=1$), and $N_{\rm fr}=N-N_{\rm oc}$ denote the number of remaining states with partial occupancies. We consider the following split in the diagonalization of ${\bf\hat{H}}$:

[TABLE]

where ${\bf Q_{oc}}$ contains the eigenvectors corresponding to $N_{\rm oc}$ eigenvalues of ${\bf\hat{H}}$, which are stored as the diagonal entries of ${\bf D_{oc}}$. On the other hand, ${\bf Q_{fr}}$ contains the eigenvectors corresponding to remaining $N_{\rm fr}$ eigenvalues of ${\bf\hat{H}}$, which are stored as the diagonal entries of ${\bf D_{fr}}$. Similarly $f\left({\bf D},\mu\right)$ can be split as

[TABLE]

Using the above equation (52) along with the scaling step in equation (50) and subspace rotation: $\bm{{\widetilde{\Psi}}}^{\bf R}=\bm{{\widetilde{\Psi}}^{{\rm o}}}\bm{{\rm Q}}$, equation (48) can be written as

[TABLE]

where $\bm{{\hat{\Psi}}^{{\rm o}}_{\bm{e}}}$ denotes the FE cell level vectors of ${\bm{{\hat{\Psi}}^{\rm o}}}={\bf M}^{-1/2}\bm{{\widetilde{\Psi}}^{{\rm o}}}$. We note that $f\left({\bf D_{oc}},\mu\right)={\bf I_{oc}}$, an $N_{\rm oc}\times N_{\rm oc}$ identity matrix and hence equation (53) can be recast in the following way:

[TABLE]

where $\bm{{\hat{\Psi}}_{{\rm fr},e}}^{\bf R}$ denotes the FE cell level vectors of $\bm{{\hat{\Psi}}_{\rm fr}}^{\bf R}={\bf M}^{-1/2}\bm{{\widetilde{\Psi}}^{{\rm o}}}\bm{{\rm Q}_{\rm fr}}$.

In the above, it is evident that the electron-density computation requires only the $N_{\rm fr}$ largest eigenstates of ${\bf\hat{H}}$. Accordingly, the spectrum-splitting based algorithms for the Rayleigh-Ritz procedure and electron-density computation in DFT-FE are given in Algorithm 4 and Algorithm 5, respectively. Even with finite-temperature Fermi-Dirac smearing, $N_{\rm fr}$ is usually a small fraction of $N$. From our numerical experiments, we find that $N_{\rm fr}$ is 10–15% of $N$ for metallic systems, and much smaller percentage ($<5\%$) for insulating and semi-conducting systems. This translates to significant cost savings in the subspace rotation step as shown in Fig. 7. This is because the usual full subspace rotation: $\bm{{\widetilde{\Psi}}}^{\bf R}=\bm{{\widetilde{\Psi}}^{{\rm o}}}\bm{{\rm Q}}$, which scales as $\mathcal{O}(MN^{2})$ is now replaced by a significantly cheaper partial subspace rotation step: $\bm{{\widetilde{\Psi}}_{\rm fr}}^{\bf R}=\bm{{\widetilde{\Psi}}^{{\rm o}}}\bm{{\rm Q}_{\rm fr}}$ (step 3 of Algorithm 4), which scales as $\mathcal{O}(MNN_{\rm fr})$. Furthermore, step 2, which now amounts to a partial diagonalization of ${\bf\hat{H}}$ to compute the $N_{\rm fr}$ largest eigenstates, can be exploited to reduce diagonalization cost. In the literature, iterative approaches like LOBPCG [41], and inner Chebyshev filtering [96] are shown to be better than ScaLAPACK’s direct eigensolver for partial diagonalization. However, iterative approaches may not be robust for metallic systems in the limit of vanishing band gaps. Hence in DFT-FE, we perform partial diagonalization using the ELPA library’s [98; 99; 100] direct eigensolver, which is more scalable than ScaLAPACK’s eigensolver and competes with the aforementioned iterative approaches with respect to minimum solution time. Fig. 7 shows the direct diagonalization times888The ELPA diagonalization times quoted here are run on NERSC Cori Intel KNL nodes which have 1.4 GHz clock frequency. On a higher clock frequency machine (eg: IBM Power and Intel Skylake architectures), these diagonalization timings are faster by a factor of 2–3 [100]. (step 2) for very large system sizes with 39,990 electrons and 61,502 electrons.

3.5.2 Mixed precision in RR

We observe that the computation of the projected Hamiltonian, ${\bf\hat{H}}=\bm{{{\widetilde{\Psi}}^{{\rm o}}}^{\dagger}}{\bf\widetilde{H}}\bm{{\widetilde{\Psi}}^{{\rm o}}}$ is the most dominant cost in the Rayleigh Ritz procedure using the spectrum-splitting technique (see Fig. 7). To this end, we develop a mixed precision algorithm to reduce the prefactor of the computation of ${\bf\hat{H}}$ and illustrate the procedure here. We first consider the following split of the orthonormalized wavefunctions $\bm{{\widetilde{\Psi}}^{{\rm o}}}$

[TABLE]

where the columns $\bm{{\widetilde{\Psi}}_{\rm oc}^{{\rm o}}}$ and $\bm{{\widetilde{\Psi}}_{\rm fr}^{{\rm o}}}$ contain the first $N_{\rm oc}$ and the remaining $N_{\rm fr}$ wavefunctions, respectively. We next rewrite the partial eigendecomposition of ${\bf\hat{H}}:{\bf\hat{H}}\bm{{\rm Q}_{\rm fr}}=\bm{{\rm Q}_{\rm fr}}\bm{{\rm D}_{\rm fr}}$ (see step 2 of Algorithm 4) as

[TABLE]

where ${\bf\hat{H}}_{\bf oc-oc}=\bm{{\widetilde{\Psi}}_{\rm oc}^{{\rm o}^{\dagger}}}{\bf\widetilde{H}}\bm{{\widetilde{\Psi}}_{\rm oc}^{{\rm o}}}$, ${\bf\hat{H}}_{\bf fr-fr}=\bm{{\widetilde{\Psi}}_{\rm fr}^{{\rm o}^{\dagger}}}{\bf\widetilde{H}}\bm{{\widetilde{\Psi}}_{\rm fr}^{{\rm o}}}$, ${\bf\hat{H}}_{\bf fr-oc}=\bm{{\widetilde{\Psi}}_{\rm fr}^{{\rm o}^{\dagger}}}{\bf\widetilde{H}}\bm{{\widetilde{\Psi}}_{\rm oc}^{{\rm o}}}$, and ${\bf\hat{H}}_{\bf oc-fr}=\bm{{\widetilde{\Psi}}_{\rm oc}^{{\rm o}^{\dagger}}}{\bf\widetilde{H}}\bm{{\widetilde{\Psi}}_{\rm fr}^{{\rm o}}}$. As the SCF approaches convergence, $\bm{{\widetilde{\Psi}}^{{\rm o}}}$ tends to the eigenfunctions of ${\bf\widetilde{H}}$, and hence the limiting behaviour of equation (62) can be written as

[TABLE]

Using equation (63) the limiting behaviour of the partial eigendecomposition of ${\bf\hat{H}}:{\bf\hat{H}}\bm{{\rm Q}_{\rm fr}}=\bm{{\rm Q}_{\rm fr}}\bm{{\rm D}_{\rm fr}}$ is written as

[TABLE]

Equation (64) provides the rationale to design a mixed precision algorithm to compute ${\bf\hat{H}}$ by employing double precision BLAS Xgemm routine to compute the ${\bf\hat{H}}_{\bf fr-fr}$ sub-matrix, while all the other sub-matrices: ${\bf\hat{H}}_{\bf oc-oc}$, ${\bf\hat{H}}_{\bf fr-oc}$ and ${\bf\hat{H}}_{\bf oc-fr}$ are computed using single precision BLAS Xgemm routine. Since $N_{\rm fr}$ is typically less than 15% of $N$, the computation of ${\bf\hat{H}}_{\bf fr-fr}$ using double precision is a very small computational cost in this approach. This leads to overall computational savings by a factor of around 2 in computation of ${\bf\hat{H}}$ as shown in Fig. 8. Additionally, in Table 4, we examine the accuracy and robustness in employing mixed precision algorithms for both Rayleigh-Ritz and orthonormalization (see Section 3.4) steps on various benchmark systems in DFT-FE. These benchmark system are chosen such that the FE discretization errors are $\sim 10^{-4}$ Ha/atom in ground-state energy and $\sim 10^{-4}$ Ha/Bohr in ionic forces. The results in Table 4 show that number of SCFs do not change between mixed precision and double precision approaches, and further the mixed precision algorithms incur negligible errors in both energies and forces in comparison to the double precision calculations. Notably, these errors are about two orders of magnitude lower than the discretization errors.

In addition to using mixed precision in computation of ${\bf\hat{H}}$, we also use a blocked approach for memory and computational efficiency improvements. The computational efficiency improvement in using the blocked approach arises from exploiting the Hermiticity of ${\bf\hat{H}}$ as shown in Fig. 4. The implementation of the blocked approach used here is similar to the implementation of the blocked approach in the overlap matrix computation (see Section 3.4). Finally, we remark that the use of spectrum splitting technique in conjunction with the mixed precision algorithm in the Rayleigh-Ritz procedure provides efficiency gains by a factor of around 3 for the large benchmark systems considered in Fig. 7 and  8.

3.6 Mixing schemes

The SCF iteration procedure for solving the Kohn-Sham eigenvalue problem can be viewed as a fixed point iteration in electron density or effective potential. In terms of electron density, this fixed point problem can be written as $\rho=F[\rho]$, where $F[\rho]$ involves computing the occupied eigenspace for a given $\rho$. This fixed point iteration can be accelerated by mixing the electron density using an appropriate mixing scheme [90; 91; 101; 102; 103; 104; 105]. In the present DFT-FE software release, we implement two kinds of mixing schemes: n-stage Anderson mixing scheme [90] and Broyden mixing scheme [91]. We have used Anderson mixing scheme in all the simulations conducted in the present work. In a future release, we plan to implement more advanced mixing strategies [101; 104], which provide improved SCF convergence rate independent of system size.

3.7 Parallelization

The primary level of parallelization in the DFT–FE code is based on domain decomposition of the adaptive FE mesh into partitions and distributing them to different MPI tasks. This is accomplished through the deal.II finite element library with p4est [106]. We note that the FE basis is localized with a compact support on the cells shared by a FE node. Hence only the FE nodes on the processor boundaries need to be communicated, which has a significantly smaller communication cost in comparison to the all-to-all communication required in global basis sets like plane-waves. This allows for excellent parallel scalability of DFT–FE, which we demonstrate subsequently in Section 4.2.

To further improve scalability, we implement a second level of parallelization over wavefunctions (band parallelization) in each of the key computational steps of the SCF iteration: Chebyshev filtering, CholGS (steps 1 and 4 of Algorithm 3), Rayleigh–Ritz procedure (steps 1 and 3 of Algorithm 4), and electron density computation. In particular, computations over the total number of wavefunctions ($\bm{\widetilde{\Psi}}$) are divided into groups of wavefunctions (band groups) and distributed among a group of MPI sub–communicators, with each sub–communicator doing computations on a single band group of size $N_{b_{i}}$, denoted by $\bm{\widetilde{\Psi}_{b_{i}}}$, where $i$ denotes the band group index. The use of band parallelization in Chebyshev filtering and electron density computation does not involve any inter band group communication of $\bm{\widetilde{\Psi}_{b_{i}}}$’s. However, the computation of the electron density requires accumulation of the electron density contribution from each band group incurring a very small communication cost. Further, we exploit band parallelization in the computation of overlap matrix and orthonormal basis construction in CholGS as ${\bf S}_{\bm{b_{i}}}=\bm{\widetilde{\Psi}^{\dagger}}\bm{\widetilde{\Psi}_{b_{i}}}$ and ${\bm{{\widetilde{\Psi}}_{b_{i}}^{\rm o}}}={\bm{\widetilde{\Psi}}}{{{\bf L^{-1}_{\bm{b_{i}}}}}^{\bf\dagger}}$, respectively. We note that ${{{\bf L^{-1}_{\bm{b_{i}}}}}^{\bf\dagger}}$ denotes $N\times N_{b_{i}}$ sub-matrix of ${\bf{L^{-1}}^{\dagger}}$. Similarly, band parallelized computation of projected Hamiltonian in Rayleigh–Ritz procedure is performed as ${\bf\hat{H}}_{\bm{b_{i}}}=\bm{{\widetilde{\Psi}}^{{\rm o}^{\dagger}}{\bf\widetilde{H}}}\bm{{\widetilde{\Psi}}_{b_{i}}^{{\rm o}}}$. We note that all-to-all communications of $\bm{\widetilde{\Psi}_{b_{i}}}$’s and $\bm{\widetilde{\Psi}^{\rm o}_{b_{i}}}$’s across band groups are performed before beginning the orthonormalization and Rayleigh-Ritz procedure, respectively. Furthermore, all-to-all communications of ${\bf S}_{\bm{b_{i}}}$’s and ${\bf\hat{H}}_{\bm{b_{i}}}$’s across band groups are also performed to compute the matrices ${\bf S}$ and ${\bf\hat{H}}$ ($N\times N$ dimensions). The above all-to-all communications involve large memory sizes, and hence the communication cost can increase significantly with increase in band parallelization groups, thus affecting parallel scaling efficiency. However, a modest amount of band parallelization can be combined with domain decomposition parallelization to extend the parallel scalability in DFT-FE, as discussed below. Additionally, we have also implemented parallelization over k points for problems involving multiple k-point sampling over the Brillouin zone in periodic calculations.

We now compare the scalability of three different parallelization approaches in DFT-FE: (i) only domain decomposition parallelization (P1), (ii) primarily band parallelization with just enough domain decomposition parallelization to fit the memory (P2), and (iii) domain decomposition parallelization till parallel scaling efficiency of $\sim 70\%$ followed by moderate band parallelization (P3). We conduct comparative studies on a large benchmark system containing 3999 atoms (39,990 electrons) using the above three approaches, and the results are shown in Fig. 9. We observe that only domain decomposition parallelization (P1) provides better parallel scalability than the primarily band parallelization approach (P2)—$73\%$ efficiency vs. $52\%$ efficiency at 51,200 MPI. This is attributed to the significant increase in MPI collective communication cost of wavefunctions as the number of band parallelization groups increase. However, the use of band parallelization is beneficial for medium–large system sizes when appreciable scaling from domain decomposition parallelization has already been extracted, as is evident from Fig. 9, where the best parallel scaling efficiency is obtained for the combined parallelization approach (P3). In particular, we use domain decomposition parallelization till 51,200 MPI tasks ($73\%$ efficiency), and band parallelization using two band parallelization groups to achieve $49\%$ efficiency at 102,400 MPI tasks. We remark that only $41\%$ efficiency is achieved for the same 102,400 tasks by solely using domain decomposition parallelization approach (P1). Based on the above comparison, in the remainder of this work, we primarily use the combined parallelization approach (P3) to scale DFT-FE calculations, particularly for large system sizes.

3.8 Geometry optimization algorithm

Geometry optimization involves relaxation of ionic forces and periodic unit-cell stresses, which are computed using the configurational forces approach as discussed in Section 2.3. In DFT-FE, geometry optimization is performed using Polak–Ribiere-Polyak non-linear conjugate gradient algorithm (PRP-CG) [107] with a secant line search. Furthermore, in the case of atomic force relaxation, we regenerate the finite-element mesh corresponding to the new atomic positions for every geometry update step. We note that the parallel adaptive mesh generation in DFT-FE, which is performed via the deal.II package is a very minor cost in comparison to the electronic ground-state solve. Additionally, we have also implemented an interface to checkpoint and restart the non-linear conjugate gradient solver for the geometry relaxation.

In the future, we plan to improve the secant line search in our non-linear conjugate gradient algorithm implementation with the more robust Brent’s method [108]. We also plan to implement alternate geometry optimization algorithms like Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) [109] and Fast Inertial Relaxation Engine (FIRE) [110], which demonstrate faster convergence than PRP-CG in many cases.

4 Results and discussion

In this section, we demonstrate the accuracy, parallel scaling performance and computational efficiency of the developed DFT-FE code on various benchmark systems involving both pseudopotential and all-electron DFT calculations. GGA [63] exchange correlation of the PBE form [66] is employed in all the calculations, and additionally ONCV [69] pseudopotentials from the SG15 database [111] are employed in all the pseudopotential DFT calculations. Further, we use Fermi-Dirac smearing with temperature $T=500$ $K$, and the $n$-stage Anderson mixing scheme [90] for mixing the electron-density in the SCF iteration procedure.

We first validate the accuracy of DFT-FE with widely used DFT basis sets on benchmark materials systems involving periodic and non-periodic pseudopotential and all-electron DFT calculations. Second, we demonstrate the parallel scalability of DFT-FE on pseudopotential benchmark systems with sizes ranging from $2,550$ electrons to $39,990$ electrons. Third, we assess the computational efficiency of the DFT-FE code by comparing to popular plane-wave based codes—QUANTUM ESPRESSO (QE) [17; 61], and ABINIT [18]—on periodic and non-periodic pseudopotential benchmark systems with sizes ranging from $2,550$ to $20,470$ electrons. Finally, we also conduct large-scale DFT calculations on sizes ranging from $27,986$ to $61,502$ electrons using DFT-FE that are computationally prohibitive using plane-wave codes. The discretization parameters in the above computational efficiency studies are chosen to be commensurate with chemical accuracy (discretization errors of $\sim 10^{-4}$ Ha and $\sim 10^{-4}$ Ha/Bohr in energy per atom and ionic forces respectively), based on the validation studies on the same benchmark systems at smaller sizes.

All the numerical simulations with computational times reported in this work were executed on the Cori supercomputer at the National Energy Research Scientific Computing (NERSC) center. In particular, we used Cori’s Phase @slowromancapii@ partition containing 9,688 compute nodes based on Intel Xeon Phi processors. Each compute node has the following specifications: single-socket Intel Xeon Phi $7250$ (“Knights Landing”) processor with $68$ physical cores per node @ $1.4$ GHz and $96$ GB of memory per node. Cori uses a Cray Aries with Dragonfly topology for inter-node communication with $45.0$ TB/s global peak bisection bandwidth.

The numerical simulations for validation of DFT-FE with respect to QE were performed on the Theta supercomputer at Argonne Leadership Computing Facility (ALCF). Theta comprises 4,392 compute nodes, with each compute node having the following specifications: single 1.3 GHz Intel Xeon Phi 7230 SKU chip (“Knights Landing”) with 64 cores, 16 GB Multi-Channel DRAM (MCDRAM) and 192 GB DDR4 memory. The interconnect topology is a dual place Dragonfly with ten groups. Each group consists of two cabinets or racks. The total bi-section bandwidth is 7.2 TB/sec.

The DFT-FE simulations reported in this section were run using $32$ MPI tasks per node and $2$ OpenMP threads (for BLAS operations), except for a few small system sizes (less than $3,000$ electrons) where using $64$ MPI tasks per node and a single OpenMP thread was found to be more efficient. Similarly, QE and ABINIT simulations were run using optimal MPI tasks-OpenMP threads combinations based on the problem size— $32$ MPI tasks per node and $2$ OpenMP threads for smaller problem sizes (less than $3,000$ electrons), and $16$ MPI tasks per node and $4$ OpenMP threads for larger problem sizes where more memory per MPI task is required. Furthermore, MPI task to core binding was also appropriately set for all the above combinations. Our numerical experiments showed that using more than $4$ OpenMP threads provided negligible performance gains in DFT-FE as well as in QE and ABINIT.

4.1 Validation

We now validate the accuracy of DFT-FE on periodic and non-periodic pseudopotential and all-electron benchmark systems. The pseudopotential benchmark systems are validated using plane-wave code QE, while all-electron benchmark systems are validated with exciting [20], which uses linearized augmented plane-wave (LAPW) method, and NWChem [25] using Gaussian basis set.

4.1.1 Pseudopotential DFT calculations

We consider three different benchmark systems: (i) hexagonal close packed (hcp) Mg periodic supercells with a mono-vacancy, (ii) body centered cubic (bcc) Mo periodic supercells with a mono-vacancy, and (iii) non-periodic Icosahedron Cu nano-particles. In each of the benchmark systems, we take increasingly refined basis sets in DFT-FE and QE and compare the ground-state energy per atom, ionic forces, and cell stresses at two different accuracy levels: a) medium accuracy—discretization errors of $\sim 10^{-4}$ Ha/atom in ground-state energy, $\sim 10^{-4}$ Ha/Bohr in ionic forces, and $\sim 5\times 10^{-6}$ Ha/${\rm Bohr}^{3}$ in cell stress (in periodic benchmark systems), which we consider as chemical accuracy; b) high accuracy—using a more refined basis set in both DFT-FE and QE to demonstrate much closer agreement between them. Further, we note that the validation studies for the aforementioned periodic benchmark systems are conducted using a Gamma point.

In the benchmark system involving Mg, we consider periodic supercells constructed from orthogonal unit cells (containing $4$ atoms) of hcp Mg with lattice constants: $a=5.882$ Bohr and $c=9.586$ Bohr. We consider a mono-vacancy in two supercells: $2\times 2\times 2$ denoted by $\textrm{Mg}_{\textrm{2x2x2}}$ containing $31$ atoms ($310$ electrons) and $4\times 4\times 4$ denoted by $\textrm{Mg}_{\textrm{4x4x4}}$ containing $255$ atoms ($2,550$ electrons). The relevant mesh parameters for DFT-FE and cut-off energies for QE are shown in Table 5. Table 5 also provides the comparison between DFT-FE and QE at medium and high accuracy levels. At the medium accuracy level, the agreement in ground-state energy is $\mathcal{O}(10^{-5})$ Ha/atom, ionic forces is $\mathcal{O}(10^{-4})$ Ha/Bohr, and hydrostatic stress is $\mathcal{O}(10^{-6})$ Hartree/${\rm Bohr}^{3}$. Similarly, at high accuracy level the agreement in ground-state energy is $\mathcal{O}(10^{-5})$ Ha/atom, ionic forces is $\mathcal{O}(10^{-5})$ Ha/Bohr, and hydrostatic stress is $\mathcal{O}(10^{-7})$ Hartree/${\rm Bohr}^{3}$. We additionally remark that Mg has a hard ONCV pseudopotential, which is reflected in basis set parameters.

Next, in the benchmark system involving Mo, we consider periodic supercells constructed from bcc Mo unit cells with lattice constant of $5.95$ Bohr. We consider a mono-vacancy in two supercell sizes—$2\times 2\times 2$ denoted by $\textrm{Mo}_{\textrm{2x2x2}}$ containing $15$ atoms ($210$ electrons) and $4\times 4\times 4$ denoted by $\textrm{Mo}_{\textrm{4x4x4}}$ containing 127 atoms (1,778 electrons). Table 5 shows the comparison between DFT-FE and QE, which demonstrates a similar excellent agreement as in the case of the Mg benchmark system. We note that the lower plane-wave cut-off or larger $h_{\rm min}$ in DFT-FE in comparison to the Mg benchmark system is attributed to Mo having a softer ONCV pseudopotential than Mg.

Finally, in benchmark system involving Cu, we consider three-dimensional non-periodic Icosahedron Cu nano-particles [112]. The Icosahedron nano-particles are constructed with nearest neighbour bond length of $6.8$ Bohr and varying the number of shells. We consider two Cu nano-particle sizes: $\textrm{Cu}_{\textrm{3-Shell}}$ containing $147$ atoms ($2,793$ electrons), $\textrm{Cu}_{\textrm{4-shell}}$ containing $309$ atoms ($5,871$ electrons). For the DFT-FE simulations, we choose a non-periodic domain containing the Cu nano-particle and impose homogeneous Dirichlet boundary conditions on the boundary of the domain. On the other hand, for the QE simulations, we choose an artificial periodic domain containing the Cu nano-particle. The energy and forces are converged with respect to the domain size in both DFT-FE and QE simulations. Table 5 shows the comparison between DFT-FE and QE, which demonstrates excellent agreement between the two codes in a non-periodic setting.

Overall, from Table 5, we show excellent agreement between DFT-FE and QE for both medium and high accuracy calculations, with the difference between the codes systematically reducing with increasing discretization.

4.1.2 All electron DFT calculations

We consider three different all-electron benchmark systems: (i) diamond cubic Si periodic unit cell, (ii) diamond cubic Si periodic supercell with a mono-vacancy, and (iii) ${\rm C}_{60}$ (Buckminsterfullerene) molecule. The first two benchmark systems, which are periodic, are validated with the exciting code using LAPW basis, while the third non-periodic benchmark system is validated with the NWChem code using Gaussian basis.

Firstly, we consider a periodic diamond cubic Si unit cell containing $8$ Si atoms ($114$ electrons) and with lattice constant $10.0065225$ Bohr. We conduct the following DFT simulations for the validation study: Gamma point denoted by $\textrm{Si}_{\rm uc}^{\rm Gamma}$, and multiple k-points using a $5\times 5\times 5$ Monkhorst-Pack Grid denoted by $\textrm{Si}_{\rm uc}^{\rm Multi-kpt}$. The mesh parameters for DFT-FE and the relevant parameters for exciting are shown in Table 6(a). In the case of exciting, rgkmax in Table 6(a) represents the product of the muffin-tin radius and the plane-wave cut-off in the interstitial region, and is the key parameter determining the accuracy of the calculation in the LAPW basis set. Further, we use the default muffin-tin radius for Si. Table 6(a) shows the comparison between DFT-FE and exciting for both $\textrm{Si}_{\rm uc}^{\rm Gamma}$ and $\textrm{Si}_{\rm uc}^{\rm Multi-kpt}$, where the agreement in ground-state energy is $\sim 1\times 10^{-4}$ Ha/atom and hydrostatic stress is $\sim 5\times 10^{-6}$ Hartree/${\rm Bohr}^{3}$.

Next, we consider a mono-vacancy in a $2\times 2\times 2$ supercell constructed from diamond cubic Si unit cell with lattice constant $10.0065225$ Bohr. This benchmark system is denoted as $\textrm{SiVac}_{2\times 2\times 2}$ and contains $63$ atoms and $882$ electrons. Table 6(a) shows the comparison between DFT-FE and exciting where the agreement in ground-state energy is $\sim 1\times 10^{-4}$ Ha/atom, ionic forces is $\sim 5\times 10^{-4}$ Ha/Bohr, and hydrostatic stress is $\sim 5\times 10^{-6}$ Hartree/${\rm Bohr}^{3}$.

Finally, we consider the non-periodic ${\rm C}_{60}$ Buckminsterfullerene molecule with hexagon-hexagon ring bond length of $2.6258$ Bohr and pentagon-hexagon ring bond length of $2.6197$ Bohr. In both NWChem and DFT-FE we use a non-periodic domain with homogeneous Dirichlet boundary conditions on the boundary of the domain, and the domain size is taken to be large enough for the boundary effects to be negligible. The appropriate mesh parameters for DFT-FE and choice of Gaussian basis for NWChem are shown in Table 6(b). The results demonstrate an excellent agreement between DFT-FE and NWChem in the ground-state energy.

Overall, from Table 6, we show good agreement between DFT-FE and exciting and also between DFT-FE and NWChem, underlining a significant capability of DFT-FE to perform highly accurate all-electron calculations with arbitrary boundary conditions.

4.2 Parallel scaling performance

Here we demonstrate the parallel scalability of DFT-FE on various system sizes. We consider hexagonal close packed (hcp) Mg periodic super cells with a mono-vacancy, and study the strong scaling behavior on three system sizes: (a) small—$\textrm{Mg}_{\textrm{4x4x4}}$ (255 atoms, 2,550 electrons), (b) medium—$\textrm{Mg}_{\textrm{8x8x8}}$ (2,047 atoms, 20,470 electrons), and (c) large—$\textrm{Mg}_{\textrm{10x10x10}}$ (3,999 atoms, 39,990 electrons). In particular, we study the strong scaling behavior by measuring the relative speedup with increasing number of MPI tasks while keeping the discretization fixed for all the three systems. The speedup is measured relative to the wall time taken on 512 MPI tasks, 3,200 MPI tasks, and 12,800 MPI tasks for $\textrm{Mg}_{\textrm{4x4x4}}$, $\textrm{Mg}_{\textrm{8x8x8}}$ and $\textrm{Mg}_{\textrm{10x10x10}}$, respectively. We note that lower number of MPI tasks were not possible due to memory constraints. This is primarily because of low memory per core ($\sim$1.4 GB) of the many-core KNL architecture in the Cori supercomputer. Further, the FE mesh in the above studies is chosen such that the discretization errors in energy and forces are $\sim 10^{-4}$ Ha per atom and $\sim 10^{-4}$ Ha/Bohr, respectively.

First, we consider the parallel scalability of the smaller system size $\textrm{Mg}_{\textrm{4x4x4}}$ as demonstrated in Fig. 10(a). Here we use only domain decomposition parallelization to scale up to 4,096 MPI tasks at 75$\%$ efficiency, with an average of 1,629 dofs per MPI task. The corresponding wall time for a single SCF iteration step is 28 seconds. Such excellent parallel scalability is possible due to the aforementioned low communication cost in the FE discretized Hamiltonian matrix and wavefunction vector products involved in Chebyshev filtering, which is the dominant computational cost for small to medium system sizes.

Next, we consider the parallel scalability of the medium to large system sizes: $\textrm{Mg}_{\textrm{8x8x8}}$ and $\textrm{Mg}_{\textrm{10x10x10}}$, which are shown in Fig. 10(b) and  10(c), respectively. Based on the comparison of three different parallelization strategies in Section 3.7 for achieving maximum parallel scalability for large system sizes, we use the parallelization strategy of combined domain decomposition and band parallelization. In particular, in the case of $\textrm{Mg}_{\textrm{8x8x8}}$, we use domain decomposition parallelization from 3,200 to 32,000 MPI tasks and then use band parallelization with two band parallelization groups to further scale to 64,000 MPI tasks at 43$\%$ efficiency. At 64,000 MPI tasks, we use an average of 1,436 dofs per MPI task, and obtain a wall time of 91 seconds for a single SCF iteration step. Similarly, in the case of $\textrm{Mg}_{\textrm{10x10x10}}$, we use domain decomposition parallelization till 51,200 MPI tasks, and then use band parallelization to further scale to 102,400 MPI tasks at 49$\%$ efficiency. At 102,400 MPI tasks, we use an average of 1,835 dofs per MPI task, and obtain a wall time of 237 seconds for a single SCF iteration step. The parallel scaling of the above medium to large system sizes is dependent on the scalability of the major computational steps: CF, CholGS and RR (section 3). We remark that in spite of the computational complexity of CF scaling quadratically in comparison to cubic scaling of CholGS and RR, CF’s excellent parallel scalability afforded by FE discretization continues to be crucial for parallel scalability at medium to large system sizes. This is evident from Fig. 11 showing the strong scaling of the various computational steps in DFT-FE for $\textrm{Mg}_{\textrm{8x8x8}}$ and $\textrm{Mg}_{\textrm{10x10x10}}$, where we note that CF is still a significant portion of the total wall-times, and further CF also demonstrates excellent parallel scalability. Fig. 11 also demonstrates good parallel scaling of CholGS and RR steps, where the use of mixed precision arithmetic based algorithms play a key role in reducing communication costs. Overall, DFT-FE’s massive parallel scalability, as demonstrated here, is a result of the locality of the FE basis as well as an effective parallel implementation of the various algorithms in DFT-FE, as discussed in section 3, that reduce communication costs and latency.

4.3 Computational efficiency and wall time comparison with plane-wave codes

We now consider three different benchmark systems with sizes ranging from $2,550$ to $61,502$ electrons to compare the computational efficiency (CPU-time and minimum wall-time) of the DFT-FE code against the plane-wave codes—QUANTUM ESPRESSO (QE) v6.3 [17; 61], and ABINIT v8.8.4 [18]. In particular, we consider: (i) hexagonal close packed (hcp) Mg periodic super cells with a mono-vacancy, (ii) body centered cubic (bcc) Mo periodic super cells with a mono-vacancy, and (iii) non-periodic Icosahedron Cu nano-particles. The details on these benchmark systems are discussed previously in Section 4.1, with a wider range of system sizes considered here. We note that the FE mesh parameters and plane-wave cut-off energies in all the benchmark calculations are chosen to be commensurate with chemical accuracy, based on the validation studies on smaller system sizes for the same benchmark system types. We do not explicitly measure the discretization errors here as highly refined calculations required to do so will use significant computational resources given the many benchmark systems along with large sizes (up to $61,502$ electrons) considered here. Further, we note that all QE and ABINIT timings reported below are for an optimal combination of FFT grid parallelization and band parallelization.

The DFT-FE simulations for the above benchmark systems use the following values of Chebyshev polynomial degree $m$ (see Section 3.3): $m=45$ for benchmark systems (i) and (ii), and $m=50$ for benchmark system (iii). Further, $N_{\rm fr}$, which is used in the RR step (see Section 3.5), is chosen to be 15 % of $N$. Additionally, in all simulations (DFT-FE, QE and ABINIT) $N$ is chosen as $N_{e}/2+b$, with $b\sim 10\%$ of $N_{e}/2$ for benchmark systems (i) and (ii), and $b\sim 5\%$ of $N_{e}/2$ for benchmark system (iii).

CPU-time comparisons

In Table 7, we compare the average computational CPU-times per SCF iteration step999Measured by taking the average of a few SCF iteration steps after the first $2-3$ SCF iteration steps, which are excluded as their timings can be variable depending on the starting wavefunctions guess to the SCF procedure. Furthermore, for few system sizes ($<10,000$ electrons) in each benchmark system, we verify that the choice of the DFT-FE parameters are adequate to achieve convergence in similar number of SCF iteration steps as taken by QE and ABINIT. We do not use advanced mixing strategies like Kerker preconditioning in QE and ABINIT simulations as such strategies are currently not implemented in DFT-FE. in the above benchmark systems between DFT-FE and the plane-wave codes QE and ABINIT. The CPU-times are reported in Node-Hrs, which is obtained by multiplying the total number of compute nodes used with the average wall-time per SCF iteration measured in hours. We also compare the number of basis functions used to achieve the desired chemical accuracy in energy and forces. We note that all the simulations for Table 7 are run using the minimum number of compute nodes required to fit the peak memory of the simulation, and additionally remark that the dashes in the table corresponding to QE and ABINIT benchmark simulations that are not performed as they are computationally prohibitive. First, we consider the periodic benchmark problems in Tables 7(a) and  7(b), where we observe that DFT-FE is more computationally efficient than both QE and ABINIT beyond system sizes of $\sim 3000$ electrons ($300$ atoms) for hcp Mg supercells, and $\sim 6,000$ electrons ($428$ atoms) for bcc Mo supercells. There are a couple of reasons for DFT-FE’s efficiency gains over QE in spite of the number of basis functions advantage of plane-wave basis for periodic problems. Firstly, simulations for medium to large system sizes require more compute nodes to fit the peak memory on many-core architectures like the Cori KNL nodes. This increases the CPU-time of plane-wave codes, relative to DFT-FE, due to the better parallel scaling of DFT-FE. Secondly, the efficient and scalable numerical implementation CF, CholGS and RR in DFT-FE (see Section 3) is also a key factor. Next, we consider the non-periodic benchmark problem: Icosahedron nano-particles of varying sizes, in Table 7(c). Here we observe that that DFT-FE is more computationally efficient than QE beyond a very small system size of $147$ atoms ($2,793$ electrons). We note that the spatial adaptivity of DFT-FE provides a key advantage in non-periodic systems where the FE mesh can be coarse-grained into the vacuum as opposed to a uniform spatial resolution of the plane-wave basis. Furthermore, the spatial adaptivity of the FE basis is also an advantage in systems having hard pseudopotentials such as Cu. Overall, from Table 7, we observe that DFT-FE’s efficiency gains over QE increases with increasing system size, achieving efficiency gains of $5.7\times$, $12.4\times$, and $11.9\times$ for $\textrm{Mg}_{\textrm{8x8x8}}$, $\textrm{Mo}_{\textrm{10x10x10}}$, and $\textrm{Cu}_{\textrm{5-shell}}$, respectively, which are the largest benchmark systems considered for CPU-time comparison. We note that ABINIT is slower than QE for all benchmark systems considered in Table 7. Finally, another key observation is that due to the efficient numerical implementation of cubic-scaling CholGS and RR steps in DFT-FE, the range of close to quadratic scaling in computational complexity with respect to number of electrons ($N_{e}$) is extended to much larger system sizes—$\mathcal{O}(N_{e}^{2.12})$ up to $N_{e}=39,990$ for hcp Mg super cells, $\mathcal{O}(N_{e}^{2.32})$ up to $N_{e}=27,986$ for bcc Mo super cells, and $\mathcal{O}(N_{e}^{2.04})$ up to $N_{e}=17,537$ for Cu nano-particles.

Wall-time comparisons

Next, in Table 8 we compare the average minimum wall-times per SCF iteration step in the above benchmark systems between DFT-FE and QE, with the restriction that the parallel scaling efficiency is above $40\%$. We observe that DFT-FE wall-times are smaller than QE wall-times for all the benchmark systems considered. Furthermore, the speedups in DFT-FE over QE increases with system size, with substantial speedups of $9\times$, $16.1\times$ and $6.9\times$ for $\textrm{Mg}_{\textrm{8x8x8}}$, $\textrm{Mo}_{\textrm{10x10x10}}$, and $\textrm{Cu}_{\textrm{5-shell}}$, respectively. Even at the smallest system sizes DFT-FE is still significantly faster than QE, with speedups of $1.5\times$, $7.5\times$ and $3.3\times$ for $\textrm{Mg}_{\textrm{4x4x4}}$, $\textrm{Mo}_{\textrm{6x6x6}}$, and $\textrm{Cu}_{\textrm{3-shell}}$, respectively. Additionally, in Table 8, we also report some very large scale simulations conducted using DFT-FE: $\textrm{Mg}_{\textrm{10x10x10}}$, $\textrm{Mo}_{\textrm{13x13x13}}$, and $\textrm{Cu}_{\textrm{6-shell}}$, obtaining very modest minimum wall-times of $203$, $277$ and $70$ seconds, respectively (with parallel scaling efficiencies above $40\%$). We note that such large system sizes are computationally prohibitive using QE, and thereby QE simulations for these systems are not performed. Finally, in Table 9 we show the breakdown of DFT-FE wall-times into key computational steps for the large system sizes in the above benchmark problems. Overall, we have demonstrated that DFT-FE is faster than plane-wave codes QE and ABINIT for system sizes beyond $2,000$ electrons with significant speedups at larger system sizes, and large-scale DFT simulations on generic material systems are practically feasible using DFT-FE for system sizes ranging up to $50,000$–$100,000$ electrons.

5 Demonstration of DFT-FE’s capabilities

5.1 Geometry optimization

We demonstrate here the capability of DFT-FE to conduct ionic relaxation using the methodology described in Section  3.8, and, further, we validate the ionic relaxation by comparing with QE. As our benchmark system, we consider tris (bipyridine) ruthenium (TBR), which belongs to a class of transition metal complexes with distinctive optical properties [113]. TBR contains $61$ atoms ($290$ electrons) in total, comprising of $30$ Carbon atoms, $24$ Hydrogen atoms, $6$ Nitrogen atoms and 1 Ruthenium atom. A schematic of the starting structure is shown in Figure 12. In both the DFT-FE and QE simulations, the discretization errors in ground-state energy and forces are converged to $\sim 10^{-4}$ Ha per atom and $\sim 10^{-4}$ Ha/Bohr, respectively. Furthermore, GGA exchange correlation and ONCV norm conserving pseudopotentials are employed along with use of Fermi-Dirac smearing with temperature $T=500$ K. In Table 10(a), we first consider the un-relaxed TBR structure, and validate that the ground-state energy and forces obtained from DFT-FE are in very good agreement with QE values. Next we conduct ionic relaxation of the TBR structure in both DFT-FE and QE until all ionic force component magnitudes are below $1\times 10^{-3}$ Ha/Bohr. The ground-state energy change of the relaxed TBR structure compared to the un-relaxed TBR structure obtained from DFT-FE and QE are in very good agreement as shown in Table 10(b). The maximum difference between relaxed and the initial un-relaxed TBR coordinates obtained from DFT-FE after removing the rigid body modes is $\Delta_{\rm max}{R^{\rm relax}}=0.276$ Bohr. Finally, using a similar metric we compare the relaxed TBR coordinates obtained from DFT-FE and QE, where we obtain $\Delta_{\rm max}{R^{\rm DFT-FE/QE}}=0.023$ Bohr. We note that only the core region of the TBR: $6$ Nitrogen atoms and 1 Ruthenium atom is considered when measuring $\Delta_{\rm max}{R^{\rm DFT-FE/QE}}$, as the outer ligand structure consisting of Carbon and Hydrogen atoms have close to zero energy modes. This coupled with the use of different relaxation solvers (BFGS in QE and CG in DFT-FE) can result in DFT-FE and QE giving significantly different coordinates for the outer ligand structure.

5.2 Molecular dynamics

We demonstrate the capability of DFT-FE to conduct Born-Oppenheimer molecular dynamics (MD) calculations [114]. In particular, we choose fcc Al $3\times 3\times 3$ supercell with fcc lattice constant $7.571944$ Bohr and containing 108 atoms and conduct a NVE molecular dynamics simulation. We use an initial ionic temperature of $T=1,500$ K and employ a time step of $0.5$ fs using a velocity Verlet time integration algorithm. This simulation is conducted using norm-conserving Troullier-Martins pseudopotential and LDA exchange correlation. We assign initial velocities using a Maxwell-Boltzmann distribution corresponding to the initial temperature of $1,500$ K. Further, we choose the electronic temperature to be equal to the ionic temperature at every step and the MD simulation is conducted until $800$ fs.

Fig. 13 shows the variation of the total energy of the system up to $800$ fs. The mean and the standard deviation of the total energy is computed to be $-2.0948$ Ha/atom and $1.482\times 10^{-5}$ Ha/atom, respectively. We also compute the drift in the total energy by evaluating the slope of linear fit, which is found to be $1.8045\times 10^{-9}$ Ha/atom-fs.

5.3 Band-structure calculations

We demonstrate the bandstructure calculation capability in DFT-FE with a diamond unit cell of silicon (lattice constant 10.26 Bohr). The Brillouin zone (BZ) is shown in Fig. 14 (a) along with two high-symmetry lines. ONCV pseudopotential is used along with PBE [66] functionals to carry out the DFT calculation. The BZ is sampled by a $8\times 8\times 8$ Monkhorst-Pack grid. The KS eigenvalues along the two high-symmetry lines are extracted from the converged Hamiltonian and are shown in Fig. 14 (b). For the sake of comparison, a similar calculation is done using QE, and the corresponding bandstructure is overlaid in Fig. 14 (b).

Defects in crystalline solids often give rise to midgap states that act as recombination centers and have drastic consequences on many applications in electronics and photonics. DFT-FE can be used to explore novel defects and their energy locations within the optical gap. The excellent scaling capability of DFT-FE with increasing system size helps to probe the physics of a system with very low defect concentrations, for example unintentional doping in semiconductors. We demonstrate a simple example of a $4\times 4\times 4$ supercell of diamond-Si with a silicon vacancy (VSi) as illustrated in Fig. 14 (c). This gives rise to spin-polarized midgap states shown on Fig. 14 (d). While the majority spin channel has a relatively deep state the minority spin has a very shallow state almost at the conduction band minima. The occupied majority spin states in the midgap lie at energy of $0.54700$ eV ($0.54714$ eV) and the corresponding empty minority spin states lie at $0.79010$ eV ($0.79016$ eV), where the top of the valence band is taken as the energy reference in both cases. The numbers within the parenthesis represent data obtained from a similar calculation in QE.

6 Conclusions

In this work, we have developed DFT-FE (Density Functional Theory with Finite-Elements), an accurate, computationally efficient and scalable finite-element (FE) based code for large-scale first-principles based materials modeling using Kohn-Sham DFT. The DFT-FE code can conduct both pseudopotential and all-electron calculations on non-periodic, semi-periodic and periodic systems, which is a unique feature that has been possible due to the real-space formulation employed in this work in conjunction with the versatility of the FE basis. Besides the systematic convergence afforded by the FE basis, the spatial adaptivity of the FE basis (realized through ‘p4est’ library in deal.II package) and the higher-order spectral finite-elements employed in DFT-FE play an important role in the computational efficiency of the code. DFT-FE offers spectral FE basis up to ${12}^{\text{th}}$ order, and the spatial adaptivity can be realized either via user defined mesh parameters, or by employing the automatic adaptive mesh refinement algorithm implemented in DFT-FE that is based on estimates of local error.

The computational efficiency and scalability of the DFT-FE code can largely be attributed to the locality of the FE basis, the algorithms employed in the solution of the discrete Kohn-Sham problem, and a careful numerical implementation of the algorithms—minimizing floating point operations, communication costs and latency—some of which leverage the attributes of the FE basis. In particular, the solution to the Kohn-Sham problem is efficiently computed by: (i) employing Chebyshev filtered subspace iteration technique to compute the eigensubspace of interest; (ii) using Cholesky factorization based Gram-Schmidt orthonormalization (CholGS) to compute an orthonormal basis spanning the subspace; (iii) employing Rayleigh-Ritz (RR) procedure to diagonalize the Hamiltonian in the projected subspace and compute the electron density to continue the self-consistent field (SCF) iteration. We developed and implemented mixed precision arithmetic based algorithms for the various steps in the solution procedure, and also employed spectrum splitting technique, that significantly reduced the computational prefactors for CholGS and RR procedure by factors of around $2$ and $3$, respectively, for large-scale systems. Since the computational complexity of CholGS and RR procedure scale cubically with system size, these efficiency gains delay the onset of cubic computational complexity in DFT-FE to very large system sizes. Our numerical studies on benchmark problems demonstrate that DFT-FE exhibits close to quadratic-scaling in system size even for those as large as $40,000$ electrons. The scalability of DFT-FE has been tested on up to $192,000$ MPI tasks, with a parallel efficiency of $42\%$ realized on a $\sim 60,000$ electron system. The locality of the FE basis is an important factor in the excellent parallel scalability of the DFT-FE code. However, a careful implementation of the various aspects of the algorithms—such as the use of elemental level matrix-matrix products and blocked approach in Chebyshev filtering, the use mixed precision arithmetic based algorithms in the CholGS and RR procedure—have been instrumental in reducing communication costs and latency, thus resulting in better scalability and reduced wall-times.

In order to assess the accuracy of DFT-FE with respect to state-of-the-art DFT codes, we have conducted a comprehensive comparison study. For validation of pseudopotential calculations we compared with QUANTUM ESPRESSO (QE), a widely used pseudopotential plane-wave code. The benchmark systems considered included periodic metallic systems and non-periodic nano-particles of various sizes. Importantly, the agreement between DFT-FE and QE on ground-state energies, ionic forces and stresses was excellent, with the differences being significantly lower than the discretization errors of the codes. In order to validate all-electron calculations, we compared with exciting for periodic calculations and NWChem for non-periodic calculations. The ground-state energies, forces and stresses obtained from DFT-FE on benchmark all-electron systems were again in excellent agreement with these codes.

In order to compare the computational efficiency and scalability afforded by DFT-FE with respect to widely used DFT codes, we conducted an extensive comparison study with QE using ONCV pseudopotentials. We used the CPU-time per SCF iteration on a wide range of benchmark problems as a metric to assess the computational efficiency, and the minimum wall-time (with at least $40-45\%$ parallel efficiency) as a metric to also assess the effectiveness of the scalability. These benchmark studies suggest that DFT-FE is more efficient than QE for periodic system sizes beyond 3,000 electrons. Furthermore, for larger systems ($10,659-20,470$ electrons), we find DFT-FE to be substantially more efficient than QE by $4.5-12\times$ in terms of CPU-times. We used the same benchmark systems to also compare the minimum wall-time per SCF iteration. Importantly, DFT-FE was significantly faster than QE for all the systems considered, with $1.5-7.5\times$ speedups for smaller system sizes ($2,550-6,034$ electrons) to $7-16\times$ speedups for larger system sizes ($14,322-20,470$ electrons). We also considered three very large periodic metallic systems with $27,986$, $39,990$ and $61,502$ electrons, where we achieve very modest average minimum wall-times per SCF iteration of $1.9$, $3.4$ and $4.6$ minutes, respectively.

DFT-FE code development has adopted best software development practices like unit and regression tests, code modularity, code documentation (using Doxygen and comments), and strict procedures for code review and merging of new changes. Such steps are very critical for sustainable and inclusive software development of DFT-FE and maintaining reproducibilty of the accuracy and performance benchmarks. We have compiled and tested DFT-FE on the following supercomputers: Flux at University of Michigan, XSEDE Comet, TACC Stampede2, NERSC Cori (both Intel KNL and Haswell nodes), ALCF Theta, OLCF Summit (IBM Power 9 CPUs), and BSC MareNostrum (Intel Xeon Platinum).

Overall, DFT-FE provides a practical capability to perform accurate massively parallel large-scale pseudopotential DFT calculations (reaching 100,000 electrons) in generic material systems with arbitrary boundary conditions and complex geometries. An attractive feature of DFT-FE is also the ability to perform all-electron calculations in the same framework, which can aid transferability studies on pseudopotentials. Further, the framework, in principle, can support mixed pseudopotential and all-electron calculations, where some atoms are treated at the all-electron level and others are treated at the pseudopotential level. This can be useful in a wide range of applications from using all-electron calculations for certain atoms with unreliable pseudopotentials to the computation of spin Hamiltonian parameters that require an all-electron treatment around the defect states (e.g. NV center in diamond) [115]. Implementation of enriched finite-element basis [38] in DFT-FE, which can enable large-scale efficient all-electron DFT calculations, is currently being pursued. Further, implementation of advanced exchange-correlation functionals (hybrid and dispersion corrected), advanced mixing schemes, spin-orbit coupling, and implementation of polarizability and dielectric calculations in DFT-FE are other efforts that are being pursued.

Acknowledgements

We thank Bikash Kanungo, Nelson Rufus and Shukan Parekh for independently testing the DFT-FE code, providing useful feedback at many stages of the code development, and providing the reference data for validating the all-electron calculations reported in this manuscript. We gratefully acknowledge the support from the Department of Energy, Office of Basic Energy Sciences, under Award numbers DE-SC0008637 and DE-SC0017380 for supporting the pseudopotential and all-electron development, respectively. We are also grateful for the support from Toyota Research Institute that funded the development efforts in the later stages of this effort. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant number ACI-1053575. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This research also used resources of Argonne Leadership Computing Facility (ALCF) resources, a DOE Office of Science User Facility. V.G. also gratefully acknowledges the support of the Army Research Office through the DURIP grant W911NF1810242, which also provided computational resources for this work. D.D. acknowledges the financial support of the German Research Foundation (DFG) grant DA 1664/2-1, that made this collaboration possible.