DFT-FE -- A massively parallel adaptive finite-element code for large-scale density functional theory calculations
Phani Motamarri Sambit Das Shiva Rudraraju Krishnendu Ghosh Denis Davydov Vikram Gavini
































a]Phani Motamarrifn1
a]Sambit Dasfn1
b]Shiva Rudraraju
a]Krishnendu Ghosh
c]Denis Davydov
a,d]Vikram Gaviniauthor
[fn1]Phani Motamarri and Sambit Das contributed equally to this work.
[author] Corresponding author.
E-mail address: vikramg@umich.edu
[a]Department of Mechanical Engineering, University of Michigan, Ann Arbor, MI 48109, USA
[b]Department of Mechanical Engineering, University of Wisconsin, Madison, WI 53706, USA
[c]Chair of Applied Mechanics, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen 91058, Germany
[d]Department of Materials Science & Engineering, University of Michigan, Ann Arbor, MI 48109, USA


We present an accurate, efficient and massively parallel finite-element code, DFT-FE, for large-scale ab-initio calculations (reaching ∼ 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.

 



Electronic structure, real-space, spectral finite-elements, ionic forces, mixed-precision arithmetic, pseudopotential, all-electron, ab-initio molecular dynamics, band structure.



 








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.

  
   






                               
                               








§ 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 <cit.> 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 <cit.> 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 <cit.> in metals and semiconductors, understanding ion conduction mechanisms and computing diffusivities in solid-state electrolytes <cit.>, and large-scale bio-molecular simulations <cit.>, 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 <cit.> or atomic-orbital type basis sets <cit.> 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 <cit.>, finite-difference <cit.>, wavelets <cit.>, psinc functions <cit.>, and other reduced order basis techniques  <cit.>.

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 <cit.> 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 <cit.>, 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 <cit.> 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  <cit.> 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  <cit.>, 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 <cit.>, 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 <cit.> 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 . 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  <cit.>. 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() 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) <cit.> 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 ∼50% on 102,400 MPI tasks for a system containing ∼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) <cit.> and ABINIT <cit.>—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  <cit.> and NWChem <cit.> 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 (∼10^-4 Ha/atom and ∼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×. 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× for larger systems. Finally, we also demonstrate the capability of DFT-FE to simulate very large systems reaching up to ∼ 60,000 electrons, which are computationally prohibitive using plane-wave codes. For the largest benchmark system containing ∼ 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 <ref> 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 <ref> 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 <ref> 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 <ref> highlights other capabilities of the DFT-FE code, and we conclude with an outlook in Section <ref>.


















     
      
    
    







§ REAL SPACE KOHN-SHAM DFT FORMULATION




 §.§ Governing equations in DFT

We consider a materials system with N_e electrons and N_a atoms whose position vectors are denoted by = {_1, _2, ⋯_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 <cit.>:

    (-1/2^2 + V_eff(ρ,)) ψ_i = ϵ_iψ_i,      i = 1,2,⋯ N   with   N > N_e/2 ,

where
ϵ_i and ψ_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 <cit.> is straightforward, and incorporated in DFT-FE. The electron density ρ in equation  (<ref>) can be expressed in terms of the orbital occupancy function  f(ϵ,μ)  and the canonical wavefunctions as

    ρ() = 2∑_i=1^N f(ϵ_i,μ)|ψ_i()|^2   .

The range of f(ϵ_i,μ) lies in the interval [0,1], and μ 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 <cit.> given by

    f(ϵ,μ) = 1/1 + exp(ϵ - μ/σ) .

In the above, σ = k_B T denotes the regularization parameter with k_B denoting the Boltzmann constant and T representing the finite temperature. We note that as σ→ 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 μ, and is given by

    ∫ρ()  = 2∑_i f(ϵ_i,μ) = N_e .

We note that f(ϵ_i,μ) is denoted as f_i in the remainder of the manuscript. The effective single-electron potential, V_eff(ρ,), in the Hamiltonian in equation (<ref>) is given by

    V_eff(ρ,) = V_xc(ρ) + V_el(ρ,) = δ E_xc/δρ + δ E_el/δρ

In the above, V_xc(ρ) denotes the exchange-correlation potential that accounts for quantum-mechanical interactions between electrons <cit.>, and is given by the first variational derivative of the exchange-correlation energy E_xc. We adopt the generalized gradient approximation (GGA) <cit.> 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

    E_xc(ρ) =∫ϵ_xc (ρ,ρ) ρ() .

Numerous forms for ϵ_xc (ρ,ρ) have been proposed, and the three widely used forms are Becke (B88)  <cit.>, Perdew and Wang (PW91)  <cit.> and Perdew, Burke and Enzerhof (PBE)  <cit.>. 

The term V_el(ρ), in the effective single-electron potential (equation (<ref>)), accounts for the electrostatic interactions. In particular, it is the variational derivative of the classical electrostatic interaction energy between electrons and nuclei, E_el, which can further be decomposed as

    E_el(ρ,) = E_H(ρ) + E_ext(ρ,) + E_zz() .

In the above, E_H, E_ext and E_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

    E_H = 1/2∫∫ρ()ρ()/| - |   ,   E_ext = - ∑_J∫ρ()  Z_J/|-_J|,  E_zz= 1/2∑_I,J ≠ IZ_I Z_J/|_I-_J|,

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 (<ref>) are over , 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  in equation (<ref>) are over the periodic domain (supercell), whereas the integrals involving  are over . 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(,) = - ∑_IZ_Iδ̃(| - _I|) with 
δ̃( - _I) denoting a regularized Dirac distribution centered at _I (and similarly b(,) = - ∑_JZ_Jδ̃(| - _J|)) to reformulate the repulsive energy  E_zz() as

    E_zz   = 1/2∫∫b(,) b(,)/| - |  - E_self
       = 1/2∫∫( b(,) b(,)/| - | - ∑_IZ_I^2δ̃(|-_I|)δ̃(|-_I|)/|-|)  ,

where E_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 𝒱_PS = 𝒱_loc + 𝒱_nl = ∑_J (𝒱_loc^J + 𝒱_nl^J), where 𝒱_loc^J and 𝒱_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, 𝒱_nl^J can be constructed as a separable pseudopotential operator  <cit.> of the form ∑_lpm|χ_lpm⟩h_lp⟨χ_lpm| , with |χ_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 𝒱_PS in the  basis, E_ext in pseudopotential Kohn-Sham DFT can be expressed as

    E_ext= 2 ∑_i=1^N∫∫ f_i  ψ^*_i()  V_PS(,,) ψ_i()  .

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

    V_nl ψ_i := ∫ V_nl(,,)ψ_i()   = ∑_J∑_lp∑_m C^J,i_lpm  h^J_lp  χ_lpm^J(,_J) ,


    with       C^J,i_lpm = ∫χ_lpm^J(,_J)ψ_i(),     1/h^J_lp = ξ_lm^Jχ_lpm^J .

In the above, |ξ_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  (<ref>) are normalized to unity. We remark that equation  (<ref>)  reduces to Troullier-Martins (TM) pseudopotential <cit.> in the Kleinman-Bylander form <cit.> for one projector component, i.e. p=1 for every l, while in the case of optimized norm conserving Vanderbilt pseudopotential (ONCV)  <cit.> 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 <cit.>.

We note that the various components of the electrostatic interaction energy in  (<ref>) and  (<ref>) are non-local in real-space, and these extended interactions are reformulated as local variational problems as discussed in   <cit.>. To this end, we define the electrostatic potential corresponding to the I^th nuclear charge Z_Iδ̃(| - _I|) to be I() and the electrostatic potential corresponding to the total charge distribution (ρ + b) to be φ(,), and these potentials are given by:

    I() =∫-Z_Iδ̃(|-_I|)/|-|,    φ(,) = ∫ρ() + b(,)/| - | .

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  (<ref>) and the expressions for different components of electrostatic energy in  (<ref>)-(<ref>), we can rewrite the electrostatic energy E_el = E_H + E_ext + E_zz as <cit.>

    E_el   = 1/2∫ (ρ()  +  b(,)) φ(,) - 1/2∑_I∫ -Z_Iδ̃(|-_I|) I()
       +∑_J∫(V^J_loc(|-_J|) -  J(|-_J|))ρ() + 2 ∑_i=1^N∫∫ f_i  ψ^*_i()  V_nl(,) ψ_i()  .

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


    (-1/2^2 +  V_xc + φ + ∑_J(V^J_loc -  J) + V_nl)ψ_i = ϵ_i ψ_i,
    
    -1/4 π^2 φ(,) = ρ() +b(,) ,    -1/4 π^2   I(,_I) = -Z_Iδ̃(|-_I|) ,
    
    2∑_if(ϵ_i,μ) = N_e ,    ρ() = 2∑_if(ϵ_i,μ)|ψ_i()|^2.


Though, the above equations  (<ref>) and  (<ref>) represent a pseudopotential treatment, we note that an all-electron treatment can be realized by setting V^J_loc =  J and V_nl = 0. We further remark that the equations  (<ref>) and  (<ref>) 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  <cit.>, 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 ψ_n,() can be expressed as

    ψ_n,() = e^i u_n,()  ,

where i = √(-1) and  u_n,() is a function that is periodic on the unit-cell while  denotes a reciprocal space point in the Brillouin zone. Using  (<ref>), the computation of electron-density in equation  (<ref>) is given by 

    ρ() = 2∑_n=1^N_BZf_n, |u_n,()|^2  ,

where _BZ denotes the volume average of the integral over the Brillouin zone corresponding to the periodic unit-cell Ω_p and f_n, denotes the orbital occupancy function corresponding to u_n,. The contribution to the electrostatic interaction energy arising from the non-local pseudopotential using the Bloch theorem can be expressed as

    E_ext^nl = 2 ∑_n=1^N_BZ ∫_ ∫_ f_n,  u^*_n,() e^-i  V_nl(,,) e^i u_n,()  .

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

    V_nl  u_n, := e^-i∫_ V_nl(,,) e^i u_n,k()  = ∑_a∑_lpm∑_r e^-i.( - _r) C_lpm^a,n h^a_lp  χ_lpm^a(,_a+_r) ,

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 in the above equation is given by

    C_lpm^a,n =  ∫_∑_re^i.( - _r) χ^a_lpm(,_a+_r)  u_n,()  .

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:


    (-1/2^2 - i  .  + 1/2 ||^2 +  V_xc + φ + ∑_J(V^J_loc -  J) + V_nl)u_n, = ϵ_n, u_n,  on  Ω_p, 
    
    -1/4 π^2 φ(,) = ρ() +b(,)  on  Ω_p ,    -1/4 π^2   I(,_I) = -Z_Iδ̃(|-_I|)  on   , 
    
    2∑_n=1^N_BZf_n, = N_e ,    ρ() = 2∑_n=1^N_BZf_n, |u_n,()|^2  .


In the above, periodic boundary conditions are imposed on Ω_p for the fields u_n(,) and φ(). Furthermore, we remark that the equations  (<ref>) 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   points in the first BZ, and this sampling of the BZ is done using the Monkhorst-pack (MP) grid  <cit.> 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  <cit.>. Finally, the set of equations  (<ref>) is solved self-consistently, with the eigenvalue problem  (<ref>) being solved for N lowest bands for every -point in IBZ. We remark that, as noted previously, periodic all-electron calculations can be realized by setting V^J_loc =  J and V_nl = 0.



 §.§ 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 <cit.>:

    ℱ_0() = min_∈ [0,1]^Nmin_∈ (𝒴)^N max_φ∈𝒴 ℒ(,,φ;)  such that  ∫ψ^*_iψ_j  = δ_ij,  2∑_if_i = N_e ,




    where   ℒ(,,φ;) = ℒ(,) 
    + min_𝒱∈ (H^1())^N_aℒ_el(,,φ,𝒱;) ,




    with  ℒ(,) = T_s(,)+ E_xc(ρ) + E_ent().

We note that = {ψ_1(),ψ_2(),ψ_3(), ⋯, ψ_N()}, and = {f_1,f_2,f_3⋯ f_N} denotes the vector of orbital occupancy factors, while
𝒱 = {V^1,V^2,⋯,V^N_a} denotes the vector containing the trial electrostatic potentials corresponding to all nuclear charges in the simulation domain. Here, T_s(,) denotes the kinetic energy of non-interacting electrons and E_ent() denotes the electronic entropy contribution, and the corresponding expressions are given by

    T_s(,) = 2 ∑_i=1^N∫ f_i ψ_i^*() (-1/2^2) ψ_i()  ,
    
     E_ent = -2 σ ∑_i=1^N[f_i ln f_i + (1- f_i)ln(1 - f_i)] .

The energy functional corresponding to electrostatic energy, ℒ_el, can be expressed in the local form as <cit.>

    ℒ_el(,,φ,𝒱,) = ∫[-1/8π |φ()|^2 + (ρ() + b(,))φ()]
       + ∑_I∫[1/8π | V^I()|^2  +  Z_I δ̃(| - _I|) V^I() ]  + ∑_J∫(V^J_loc(|-_J|) -  J(|-_J|))ρ()
        + 2 ∑_i=1^N∫∫ f_i  ψ_i^*()V_nl(,,) ψ_i() ,

where J denotes the electrostatic potential corresponding to the J^th nuclear charge (see equation (<ref>)), or analogously 𝒱̅_δ = { 1, 2,⋯,N_a}=arg  min_𝒱∈ (H^1())^N_aℒ_el(,,φ,𝒱;). Further, we note that, 𝒴 in equation  (<ref>) 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-cell[while the variational problem in equation (<ref>) 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  <cit.>.]. Denoting such an appropriate bounded domain by Ω subsequently, 𝒴=H^1_0(Ω) in the case of non-periodic problems, and 𝒴=H^1_per(Ω) in the case of periodic problems.



 §.§ Configurational Forces

We employ configurational forces approach <cit.> 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  (<ref>) with respect to the position of a material point . 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  (<ref>) with respect to . We refer to <cit.> for the derivation of configurational forces, and a comprehensive discussion on this topic. 

Let : → represent the infinitesimal perturbation of the underlying space, mapping a material point  to a new point ' such that ^0 = I. Further, let the generator of this mapping be = such that  is constrained to rigid body deformations in the compact support of the regularized nuclear charge distribution () in order to preserve the integral constraint ∫δ̃( - _I)  = 1. Denoting ℱ_0() to be the ground-state free energy in the perturbed space, the configurational force is evaluated by computing the Gâteaux derivative of  ℱ_0() given by

    dℱ_0()/dε = ∫_Ω: ()  + ∑_I∫_'^I: ()  + F^PSP ,


where  and  '^I denote the Eshelby tensors whose expressions in terms of the solutions of the saddle point problem  (<ref>) on the original space are provided below.  If = {_1,_2,_3 ⋯_N } denote the Kohn-Sham eigenfunctions corresponding to the lowest N eigenvalues ϵ̅={ϵ̅_1,ϵ̅_2,⋯,ϵ̅_N} with occupancies ={_1,_2 ⋯_N}, and  denotes the electrostatic potential (all solutions of the saddle point problem  (<ref>)), the expressions for the Eshelby tensor  and F^PSP in equation  (<ref>) are given by

    = (∑_i=1^N(_i_i^*()·_i() - 2  _i ϵ̅_i ^*_i() _i())  + ε_xc(,)()  -1/8π |()|^2 + () ()  
       + ∑_J (V_loc^J -  J)()  +  E^nl + E^nl^*) - ∑_i=1^N_i[_i^*() ⊗_i()
     + _i() ⊗_i^*()]  
        - ∂/∂ρ(ε_xc(,)()) ⊗  + 1/4π() ⊗() ,
       '^I = 1/8π |_^I()|^2  - 1/4π_^I() ⊗_^I () ,

where



    () = 2∑_i_i ^*_i() _i() ,    _i = 1/1 + exp(ϵ̅_i -μ/k_B T) ,
    E^nl = 2 ∑_i=1^N∑_J∑_lpm_i  h_lp^J _i^*() χ^J_lpm(,_J) ∫_Ωχ^J_lpm(,_J) _i()  .

Further,

    F^PSP = ∑_J∫_Ω()((V^J_loc(|-_J|) -  J(|-_J|)))·(() - (_J) ) + F^nl + F^nl^* ,

where

    F_nl = 2 ∑_i=1^N∑_J∑_lpm_i h_lp^J[∫_Ω_i^*() χ^J_lpm(,_J)(() - (_J))]
    [∫_Ωχ^J_lpm(,_J)_i() ]  .

Though equation (<ref>) represents the pseudopotential case, the expression for an all-electron case is realized by setting E^nl = 0, F^PSP = 0 and  V^J_loc(|-_J|) =  J(|-_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 (). In particular, the force on any given atom is computed via equation  (<ref>) by choosing the compact support of () to contain the atom of interest, and the periodic unit-cell stress tensor is evaluated by choosing _i = C_ijx_j in   (<ref>) as explained in  <cit.>.






 §.§ 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 <cit.>. 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 <cit.> 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  (<ref>) 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 𝒯^h and 𝒯^h_el with the characteristic mesh-sizes denoted by h and h_el, respectively. We consider 𝒯^h_el to be a uniform subdivision of 𝒯^h. Denoting the subspaces spanned by the FE basis corresponding to triangulations 𝒯^h and 𝒯^h_el to be 𝕍_h^M (with dimension M) and 𝕍_h_el^M_el (with dimension M_el>M), we note that 𝕍_h^M⊂𝕍_h_el^M_el. Finally, the representation of the various fields in the Kohn-Sham problem (<ref>)—the wavefunctions and the electrostatic potentials—in the FE basis is given by

    ψ^h_i() = ∑_j=1^M N^h_j() ψ^j_i  ,    φ^h_el() = ∑_j=1^M_el N^h_el_j() φ^j ,    J^h_el() = ∑_j=1^M_el N^h_el_j()J^j ,

where N^h_j:1≤ j ≤ M denotes the FE basis spanning 𝕍_h^M and N^h_el_j:1≤ j ≤ M_el denotes the FE basis spanning 𝕍_h_el^M_el . We note that ψ^h_i, φ^h_el and J^h_el denote the FE discretized fields, with ψ^j_i, φ^j and J^j 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 (<ref>) results in a generalized eigenvalue problem given by

_i = ϵ^h_i_i

where   denotes the discrete Hamiltonian matrix with matrix elements H_jk,  denotes the overlap matrix (or commonly referred to as the mass matrix in finite element literature) with matrix elements M_jk, and ϵ^h_i denotes the i^th eigenvalue corresponding to the discrete eigenvector _i. The expression for the discrete Hamiltonian matrix, H_jk = H_jk^loc + H_jk^nl, is given in terms of  

    H_jk^loc = 1/2∫_Ω N^h_j () . N^h_k()  +∫_Ω V^h_eff,loc(,)  N^h_j() N^h_k()  .





In the above, V_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_xc^h, total electrostatic potential φ^h_el() and the local pseudopotential term as follows:

    V_eff,loc^h(,) = V_xc^h() + φ^h_el() + ∑_J (V^J_loc(|-_J|) - J^h_el(|-_J|)) .


In the case of all-electron calculations, V_eff,loc^h(,) = V_xc^h() + φ^h_el() and H_jk^nl is zero. In the case of pseudopotential calculations, H_jk^nl is given by

    H_jk^nl = ∑_J = 1^N_a∑_lpm C^J_lpm,jh^J_lpC^J_lpm,k ,  where   C^J_lpm,j =  ∫_Ωχ^J_lpm(,_J) N^h_j()  .

Finally, the matrix elements of the overlap matrix  are given by M_jk=∫_Ω N^h_j() N^h_k().
We note that the matrices ^loc and  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 ^nl are also sparse since the projectors χ^J_lpm(,_J) have a compact support, thus rendering a sparse structure to the discrete Hamiltonian .

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

    _i = ϵ^h_i_i  ⇒    _i = ϵ^h_i^1/2^1/2_i  ⇒    _i = ϵ^h_i_i ,


    where   _i = ^1/2_i ,     = ^-1/2^-1/2 .

We note that  is a Hermitian matrix, and (<ref>) represents a standard Hermitian eigenvalue problem. The actual eigenvectors are recovered by the transformation _i = ^-1/2_i. 
Furthermore, we note that the matrix ^-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 <cit.>. This renders the matrix   the same sparsity structure as the matrix .

Finally, for the given positions of nuclei, the discrete Kohn-Sham eigenvalue problem along with the discretized Poisson equations for the electrostatic potentials (φ^h_el and J^h_el) are to be solved self-consistently, and are given by:


    ^-1/2^-1/2_i = ϵ^h_i_i ,
    ∑_j=1^M_el[1/4π∫_Ω∇ N_i^h_el(). ∇ N_j^h_el() ]φ^j
    =∫_Ω(ρ^h() + b^h_el(,))N_i^h_el( x)   d x ,   
    ∑_j=1^M_el[1/4π∫_Ω_J∇ N_i^h_el(). ∇ N_j^h_el()  d x]V̅^J^j_δ̃
    =∫_Ω_J( b^h_el_J(| - _J|) )N_i^h_el( x)   ,   ∀ J  , 
    
    
    2∑_if(ϵ^h_i,μ) = N_e ,    ρ^h() = 2∑_if(ϵ^h_i,μ)|ψ^h_i()|^2  .


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(,) in equation (<ref>) is given by b^h_el(,) = ∑_I b_I^h_el(| - _I|) with b_I^h_el = - Z_Iδ(| - _I|) where δ(| - _I|) denotes the Dirac-delta distribution centered at the position of the atom _I.  The boundary conditions used for the computation of the discrete potential field φ^h_el() in equation  (<ref>) 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 V̅^J^h_el_δ̃ associated with individual nuclear charge J  is solved using the discrete Poisson equation  (<ref>) subject to Dirichlet boundary conditions  with prescribed Coulomb potential applied on a domain Ω_J enclosing the atom J.  After obtaining the electronic ground-state from the solution of the discrete Kohn-Sham problem (equations (<ref>)), we compute the discrete total ground-state energy  E^h  in terms of the discrete solution fields (ϵ̅^h_i, ^h, ^h_el,I^h_el) as follows:

    E^h = E^h_band - E^h_pot +  E^h_xc(^h,ρ^h) + E^h_el_el ,

where, 

    E^h_band = 2 ∑_i f(^h_i,μ) ^h_i ,    
     E^h_pot  = ∫_Ω^h() (V^h_xc() + ^h_el())  ,
    
    E^h_el_el   =  ∫_Ω[-1/8π |^h_el()|^2 + (^h_el() + b^h_el(,))^h_el()]
       + ∑_I∫_[1/8π |I^h_el()|^2  -  b^h_el_I(| - _I|) I^h_el() ] .






§ NUMERICAL IMPLEMENTATION


DFT-FE is built over the deal.II open-source finite-element library <cit.>, 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. 


 §.§ 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  <cit.> as well as the development of coarse-graining techniques that seamlessly bridge electronic structure calculations with continuum <cit.>. 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  <cit.> via deal.II open-source finite-element library <cit.>. Spatial discretization via this octree refinement produces a non-conforming mesh (cf. Fig. <ref>), 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.



  §.§.§ 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_atom and coarsens away.  This is accomplished by discretizing the given simulation domain using an initial coarse uniform mesh with mesh size h_base.  Subsequently, the FE cells whose centroids are within a distance of r_atom units from each of the atomic positions are marked for refinement. These marked cells are refined until a target mesh size h_atom is achieved. Hence, the mesh size parameters h_base,  h_atom  and the parameter r_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_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 SiF_4 molecule for all-electron DFT calculation is shown in Fig. <ref> .

We note that the mesh parameters h_base, h_atom, h_fine for a given materials system can be estimated from the optimal mesh size distribution h() 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  <cit.> to determine this mesh-size distribution. Let E be the ground-state energy for the continuous problem  (<ref>) with = {_1 , _2  ⋯ _N} representing the electronic wavefunction solutions of the continuous problem  (<ref>), and, further, let E_h denote the discretized ground-state energy with  ^h = {_1^h , _2^h ⋯ _N^h} representing the electronic wavefunction solutions of the discrete problem (<ref>). As discussed in the previous section, we note that the triangulation 𝒯^h_e employed in the solution of the electrostatics problems  (<ref>) and  (<ref>) is more resolved than the triangulation 𝒯^h used for solving Kohn-Sham eigenvalue problem  (<ref>), and is usually chosen to be the uniform subdivision of the triangulation 𝒯^h. Thus, the dominant discretization error will correspond to the error in the wavefunctions, and as derived in equation (47) of  <cit.>, the discretization error |E - E_h| (retaining only the leading order terms) is given by 

    |E - E^h| ≤ C  (∑_i|_i - ^h_i|^2_1,Ω) ≤𝒞∑_eh_e^2k[∑_i|_i|_k+1,Ω_e^2] ,

where e denotes an element in the FE discretization, with mesh size h_e covering the domain Ω_e, and | . |_k+1,Ω_e denotes the (k+1) semi-norm over Ω_e. Using the definition of semi-norms, equation (<ref>) can be written as

    |E_h - E| ≤𝒞∑_e[ h_e^2k∫_Ω_e[∑_i|D^k+1_i()|^2] ] ≤𝒞' h^2k()[∑_i|D^k+1_i()|^2]  ,

where h() denotes the element size distribution function defining the target element size at point  in the simulation domain. The mesh distribution h() can be estimated by minimizing the approximation error in energy in  (<ref>) subject to fixed number of elements, which is given by

    min_h{ h^2k(  )[∑_i|D^k+1_i()|^2] }  /h^3()=N_E .

The solution to this variational problem is given by

    h() = A (∑_i|D^k+1_i()|^2)^-1/(2k+3) ,

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() in equation  (<ref>) involves the knowledge of _i(), 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_base, h_atom, h_fine can be estimated from h(), 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_atom and h_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. 







  §.§.§ Automatic adaptive mesh refinement (AAMR)

A number of recent works <cit.> 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  <cit.>. Further, error indicators based on H^1 semi-norms of the wavefunctions <cit.> and those based on coarsening mesh approaches <cit.> 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 [computational 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 (<ref>) that involves the continuous wavefunctions. We first derive this energy error estimate following the mesh adaption ideas in  <cit.>, and subsequently present the algorithm implemented in DFT-FE for automatic adaptive mesh refinement.





















Let ^h()  be the (k-1)^th interpolant of ^h() with k≥ 2 denoting the FE interpolating polynomial. As opposed to the approximation error |E-E^h|, we will work with |E-E^h| where E^h denotes the discrete Kohn-Sham ground-state energy obtained using the (k-1)^th interpolants of the FE solution of the Kohn-Sham problem (equation (<ref>)). To this end, following along similar lines as the derivation of equation  (<ref>), the dominant term in |E-E^h| can be derived to be

    |E - E^h| ≤C(∑_i|_i - ^h_i|^2_1,Ω) .

We recall, in the above expression, | . |_1,Ω denotes the H^1 semi-norm over Ω. Using triangle inequality, we have

    |_i - ^h_i|_1,Ω = |(_i - ^h_i) + (^h_i - ^h_i)|_1,Ω≤|_i - ^h_i |_1,Ω + |^h_i - ^h_i|_1,Ω

We note that, asymptotically, as h → 0, the first term,  |_i - ^h_i|_1,Ω , in the above inequality (<ref>) is h^k  <cit.>, while the second term, |^h_i - ^h_i|_1,Ω , is of order h^k-1. Hence, the error is dominated by |^h_i - ^h_i|_1,Ω , which in turn can be bounded in terms of the FE mesh size as <cit.>

    |_i - ^h_i|_1,Ω≤Ĉ|^h_i - ^h_i|_1,Ω≤Ĉ^'∑_e (h_e)^k-1|^h_i|_k,Ω_e ,

where e denotes an element in the FE discretization, with mesh size h_e covering the domain Ω_e. Hence, the energy error estimate |E - E^h| in equation  (<ref>) can be written as

    |E - E^h| ≤ C  ∑_e  (h_e)^2(k-1)(∑_i |^h_i|^2_k,Ω_e) .


The energy error estimate  (<ref>) motivates r_e= (h_e)^2(k-1)∑_i |D^k_i^h()|^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 ψ_i^h()) as good approximations to _i^h() 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 𝒯^h^(0), and a sequence of nested triangulations 𝒯^h^(n) are generated using the following iterative procedure (Algorithm <ref>): Interpolate single-atom DFT wavefunctions → Estimate local error → Mark for refinement → Refine → Check convergence.
In particular, the triangulation 𝒯^h^(n+1) is constructed from 𝒯^h^(n) by interpolating the single-atom DFT wavefunction data onto the current mesh 𝒯^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 (β) 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 𝒯^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, ∑_i ∫ |ψ_i^h^(n)()|^2.  This is motivated from  (<ref>), 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 SiF_4 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 <ref>–<ref>). We note that the reference data for the ground-state energy per atom (E_0), force vector (_0) and hydrostatic stress (σ_0) to measure the discretization errors are obtained using highly refined FE calculations. These results indicate that AAMR procedure as described in Algorithm <ref> 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 SiF_4 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 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 10^-4 in ground-state energies and forces. Fig. <ref> 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 SiF_4 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 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.  <ref> and  Fig. <ref> 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. <ref> demonstrates that absolute discretization errors as low as 10^-5 Ha/atom in the ground-state energy and 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. <ref> and this behavior is a consequence of using a priori error indicator involving single atom wavefunctions in AAMR.
































 §.§ SCF Algorithm


The discrete nonlinear Hermitian eigenvalue problem is solved self-consistently along with Poisson equations (see equation  (<ref>)) to compute the Kohn-Sham ground-state solution. Algorithm <ref> 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) <cit.> to evaluate the occupied eigenspace of the discrete Kohn-Sham Hamiltonian. We further employ Anderson and Broyden schemes <cit.> for electron-density mixing, and the finite-temperature Fermi-Dirac smearing <cit.> to avoid the charge sloshing associated with metallic systems. 

The ChFSI procedure in Algorithm <ref> 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 <cit.> 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 <cit.>.





 §.§ Chebyshev filtering

In DFT-FE, Chebyshev polynomial filtering technique <cit.> is used to adaptively approximate the wanted eigenspace (the lowest N occupied eigenfunctions) of the FE discretized Hamiltonian H <cit.>. In practice, N is typically chosen as N_e/2+ b to allow for finite-temperature Fermi-Dirac smearing, where b is usually (5-10) % of N_e/2. In a given SCF iteration step, a scaled Hamiltonian H̅ is obtained by scaling and shifting H  such that the unwanted spectrum of H is mapped on to [-1, 1], and the wanted spectrum is mapped on to (-∞,-1) to exploit the fast growth property of Chebyshev polynomials in this region. Subsequently, the action of a degree m Chebyshev polynomial filter, T_m(H̅), on the input subspace, Ψ, is computed recursively as

    T_m(H̅)  Ψ=[2H̅ T_m-1(H̅)-T_m-2(H̅)]Ψ .

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 δ, chosen to between 1 × 10^-2 - 5 × 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 (<ref>) is based on the upper bound of the spectrum of 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 (∼ 500-1000) are required for all-electron calculations.   




  §.§.§ Practical implementation aspects of Chebyshev filtering

The computational complexity of Chebyshev filtering scales as 𝒪(MN), where M is the size of the discretized Hamiltonian 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 H̅ X in equation (<ref>), with 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 X. We then employ BLAS  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 limited[The 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 <cit.>.] and incur a higher communication cost[The 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 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, 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  [Batched 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 H 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  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  H 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  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. 



 §.§ 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 <cit.> than the commonly used classical Gram-Schmidt procedure. Algorithm <ref> shows the steps involved in the CholGS procedure. The 𝒪(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 S in CholGS has a much lower communication latency in comparison to 𝒪(N^2) communication calls in classical Gram-Schmidt.













  §.§.§ Parallel implementation aspects of CholGS in Algorithm <ref>



  
Computation of overlap matrix
We first note that Ψ is stored in parallel as a M_ loc× N matrix, where M_ loc is the number of FE nodes owned locally by a given MPI task. Accordingly, a straightforward approach to compute the overlap matrix  S in step 1 involves the evaluation of local contributions of Ψ^†Ψ (a N × N matrix) on each MPI task, and then accumulating the local contributions to  S using the  collective routine. However, this approach requires memory corresponding to a N × N matrix on each MPI task, and hence is not practically applicable for large-scale problems ∼(N>20,000). To avoid this large memory footprint in both storage of S as well as computation of the local contributions, we use the popular 2D cyclic block grid distribution of ScaLAPACK library <cit.> to distribute the memory of S, and use a blocked approach to compute the local contributions of Ψ^†Ψ to S. Further in the blocked approach, we also exploit the Hermiticity of S, by computing only the lower triangular portion of S. Fig. <ref> shows the schematic of the blocked approach with block size B^v, where [i,N] × [i,i+B^v] sub-matrices of 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] × M_ loc block of Ψ^† and M_ loc× [i,i+B^v] of Ψ using BLAS  routine, followed by accumulation of the local contributions using the  collective. Subsequently, the corresponding sub-matrix entries of the ScaLAPACK parallelized S are filled. Overall, the above blocked approach combined with ScaLAPACK parallelization of S provides both memory optimization and efficiency improvements. 




  
Computation of inverse of Cholesky factor
Cholesky factorization of S in step 2 and inversion of the Cholesky factor L in step 3 are performed using ScaLAPACK routines  and , 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. <ref>). 

 


  
Construction of orthonormal vectors
Similar to step 1, computation of the orthonormalized basis Ψ^ o in step 4 also has a large memory footprint when performed simply as a matrix-matrix multiplication between the local portion (M_ loc× N matrix) of the parallel distributed Ψ and the full L^-1^† (N × N matrix) on every MPI task. For large-scale problems this leads to a high peak memory due to storage of the full L^-1^† on every MPI task, and also to store the computed Ψ^ o, which requires the same memory size as Ψ. Hence, we compute Ψ^ o using two blocked levels to address both of these memory issues, as shown schematically in Fig. <ref>. First, we employ an outer blocked level over M_ loc with block size B^d, which allows reuse of the memory of Ψ to store Ψ^ o. In particular, we compute [i,i+B^d] × [1,N] sub-matrices of Ψ^ o one after the other and copy the orthonormalized sub-matrices back on to Ψ, thereby requiring only an additional B^d × N memory. Secondly, we employ an inner blocked level where each [i,i+B^d] × [1,N] sub-matrix of Ψ^ o is further divided into [i,i+B^d] × [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 L^-1^† while also exploiting the triangular matrix property of L^-1^†. Each [i,i+B^d] × [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] × [1,j+B^v] sub-matrix of Ψ and a [1,j+B^v] × [j,j+B^v] sub-matrix of L^-1^†. We note that L^-1^† is stored in a ScaLAPACK parallel format after the end of step 3. Thus to obtain the [1,j+B^v] × [j,j+B^v] sub-matrix of L^-1^† in each MPI task, we first use the local portion of the parallel L^-1^† to fill the corresponding entries in the sub-matrix and the rest as zeros, and subsequently use the  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  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 S in step 1 and the triangular matrix nature of L^-1^† 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 L^-1^† 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_ loc is much larger than B^d.







  §.§.§ 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 <ref>, 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 <cit.>. 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:

    S= S_d +  S_od,

where S_d is a matrix containing the diagonal entries of S. We take advantage of the fact that S_od→0 as the SCF approaches convergence and hence compute S_od using single precision BLAS  routines, while the computation of diagonal entries of S_d is performed using double precision BLAS routines at negligible computational cost.  Similarly, step 4 can be split into

    Ψ^ o=ΨL^-1_d^†+ΨL^-1_od^†,

where L^-1_d^† is a matrix containing the diagonal entries of L^-1^†.  Taking advantage of the fact that L^-1_od^†→0 as the SCF approaches convergence, we compute ΨL^-1_od^†
using single precision BLAS  routines, while the computation of ΨL^-1_d^† 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  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. <ref> 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 <ref>, and is discussed in detail subsequently.







 §.§ Rayleigh-Ritz procedure and electron-density computation

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

    ρ^h_ out( x)= 2∑_i=1^N f(ϵ^h_i,μ)|ψ^h_i()|^2
    =2∑_i=1^N[f(ϵ^h_i,μ) (∑_j=1^M_eψ_i^e,j N^e_j( x) ) ( ∑_k=1^M_eψ^e,k^*_i N^e_k( x) ) ],

where {N^e_1( x) , N^e_2( x),  ⋯  N^e_M_e( x)} denote the FE basis functions associated with the given cell (M_e denoting the number of nodes in the cell),  and  {ψ_i^e,1,  ψ_i^e,2,  ⋯ ψ_i^e,M_e} denote the corresponding nodal values of the i^ th wavefunction, ψ^h_i( x) in FE cell e. 





Using the subspace rotated wavefunctions Ψ^ R, Equation <ref> can be re-written as

    ρ^h_ out( x)=2  n^e^T( x)
     Ψ̂^ R_e f (  D, μ) Ψ̂_e^ R^†n^e( x),

where  
    n^e( x)=[N^e_1( x)  N^e_2( x)  ⋯  N^e_M_e( x)]^T,

and the matrix Ψ̂^ R_e contains the FE cell column vectors extracted from Ψ̂^ R which is given by

    Ψ̂^ R=  M^-1/2Ψ^ R.















In the above, the computational complexity of steps i), ii) and iii) of the Rayleigh-Ritz procedure scales as 𝒪(MN^2), 𝒪(N^3), and 𝒪(MN^2), respectively, while the electron-density computation scales as 𝒪(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.



  §.§.§ Spectrum-splitting in RR



The key idea behind spectrum-splitting is that the eigenvalues and eigenvectors of the projected Hamiltonian Ĥ with orbital occupancy function f_i=1 are not explicitly necessary for the computation of the electron-density in equation (<ref>). 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 ∼ 500   K.  Such methods have been developed in previous works in the context of both pseudopotential <cit.> and all-electron DFT <cit.> 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_ oc denote the number of Kohn-Sham eigenstates with full occupancies (f_i=1), and N_ fr = N - N_ oc denote the number of remaining states with partial occupancies. We consider the following split in the diagonalization of Ĥ:

    Ĥ
    =
    [
    [ Q_oc Q_fr ]] 
    [
    [ D_oc    0;    0 D_fr ]]
    [
    [ Q_oc^†; Q_fr^† ]],

where Q_oc contains the eigenvectors corresponding to N_ oc eigenvalues of Ĥ, which are stored as the diagonal entries of D_oc. On the other hand, Q_fr contains the eigenvectors corresponding to remaining N_ fr eigenvalues of Ĥ, which are stored as the diagonal entries of D_fr. Similarly f (  D, μ) can be split as

    f (  D, μ)
     =
      [
     [ f ( D_oc,μ)           0;           0 f ( D_fr,μ) ]].

Using the above equation (<ref>) along with the scaling step in equation (<ref>) and subspace rotation: Ψ^ R=Ψ^ o Q, equation (<ref>) can be written as

    ρ^h_ out( x)=2 n^e^T( x)
       Ψ̂_ e^ o[
    [ Q_oc Q_fr ]] 
      [
     [ f ( D_oc,μ)           0;           0 f ( D_fr,μ) ]]
    [
    [ Q_oc^†; Q_fr^† ]]
    Ψ̂_e^ o^†n^e( x),

where Ψ̂^ o_ e denotes the FE cell level vectors of Ψ̂^ o= M^-1/2Ψ^ o.
We note that f ( D_oc,μ)= I_oc, an N_ oc× N_ oc identity matrix and hence  equation (<ref>) can be recast in the following way:

    ρ^h_ out( x)=     2  n^e^T( x)
     Ψ̂_ e^ o[
    [ Q_oc Q_fr ]] 
      [
       I+
      (
     [                 0                 0;                 0 f ( D_fr,μ)- I_fr ])
    ]
    [
    [ Q_oc^†; Q_fr^† ]]
    Ψ̂_e^ o^†n^e( x)
    
    =    2  n^e^T( x)[  Ψ̂_ e^ oΨ̂_ e^ o^†
    + Ψ̂^ o_ e Q_fr(f( D_ fr,μ) - I_fr) Q_fr^†Ψ̂_ e^ o^†]n^e( x)
    
    =    2  n^e^T( x)[  Ψ̂_ e^ oΨ̂_e^ o^†
    + Ψ̂_ fr,e^ R(f( D_ fr,μ) - I_fr)Ψ̂_ fr,e^ R^†]n^e( x),

where Ψ̂_ fr,e^ R denotes the FE cell level vectors of Ψ̂_ fr^ R=  M^-1/2Ψ^ o Q_ fr.


In the above, it is evident that the electron-density computation requires only the N_ fr largest eigenstates of Ĥ. Accordingly, the spectrum-splitting based algorithms for the Rayleigh-Ritz procedure and electron-density computation in DFT-FE are given in Algorithm <ref> and Algorithm <ref>, respectively. Even with finite-temperature Fermi-Dirac smearing, N_ fr is usually a small fraction of N. From our numerical experiments, we find that N_ 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. <ref>. This is because the usual full subspace rotation: Ψ^ R=  Ψ^ o Q, which scales as 𝒪(MN^2) is now replaced by a significantly cheaper partial subspace rotation step: Ψ_ fr^ R= Ψ^ o Q_ fr (step 3 of Algorithm <ref>), which scales as 𝒪(MNN_ fr). Furthermore, step 2, which now amounts to a partial diagonalization of Ĥ to compute the N_ fr largest eigenstates, can be exploited to reduce diagonalization cost. In the literature, iterative approaches like LOBPCG <cit.>, and inner Chebyshev filtering <cit.> 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 <cit.> direct eigensolver, which is more scalable than ScaLAPACK's eigensolver and competes with the aforementioned iterative approaches with respect to minimum solution time. Fig. <ref> shows the direct diagonalization times[The 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 <cit.>.] (step 2) for very large system sizes with 39,990 electrons and 61,502 electrons. 












  §.§.§ Mixed precision in RR


We observe that the computation of the projected Hamiltonian, Ĥ=Ψ^ o^†HΨ^ o is the most dominant cost in the Rayleigh Ritz procedure using the spectrum-splitting technique (see Fig. <ref>). To this end, we develop a mixed precision  algorithm to reduce the prefactor of the computation of Ĥ and illustrate the procedure here. 

We first consider the following split of the orthonormalized wavefunctions  Ψ^ o

    Ψ^ o=[Ψ_ oc^ o Ψ_ fr^ o],

where the columns Ψ_ oc^ o and Ψ_ fr^ o contain the first N_ oc and the remaining N_ fr wavefunctions, respectively. We next rewrite the partial eigendecomposition of  Ĥ:Ĥ Q_ fr= Q_ fr D_ fr (see step 2 of Algorithm <ref>) as

    [
    [ Ĥ_ oc-oc Ĥ_ oc-fr; Ĥ_ fr-oc Ĥ_ fr-fr ]]
    [
    [ Q_ fr^a; Q_ fr^b ]]
    =
    [
    [ Q_ fr^a; Q_ fr^b ]]
     D_ fr,

where Ĥ_ oc-oc=Ψ_ oc^ o^†HΨ_ oc^ o,  Ĥ_ fr-fr=Ψ_ fr^ o^†HΨ_ fr^ o,  Ĥ_ fr-oc=Ψ_ fr^ o^†HΨ_ oc^ o, and Ĥ_ oc-fr=Ψ_ oc^ o^†HΨ_ fr^ o. As the SCF approaches convergence, Ψ^ o tends to the eigenfunctions of H, and hence the limiting behaviour of equation (<ref>) can be written as

    [
    [ Ĥ_ oc-oc→ D_oc    Ĥ_ oc-fr→ 0;    Ĥ_ fr-oc→ 0 Ĥ_ fr-fr→ D_fr ]]
    [
    [    Q_ fr^a→ 0; Q_ fr^b→ I_fr ]]
    = 
    [
    [    Q_ fr^a→ 0; Q_ fr^b→ I_fr ]]
     D_ fr.

Using equation (<ref>) the limiting behaviour of the partial eigendecomposition of  Ĥ:Ĥ Q_ fr= Q_ fr D_ fr is written as

    Ĥ_ fr-fr Q_ fr^b =  Q_ fr^b D_ fr.

Equation (<ref>) provides the rationale to design a mixed precision algorithm to compute Ĥ by employing double precision BLAS  routine to compute the Ĥ_ fr-fr sub-matrix, while all the other sub-matrices: Ĥ_ oc-oc, Ĥ_ fr-oc and Ĥ_ oc-fr  are computed using single precision BLAS  routine. Since N_ fr is typically less than 15% of N, the computation of Ĥ_ 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 Ĥ as shown in Fig. <ref>. Additionally, in Table <ref>, we examine the accuracy and robustness in employing mixed precision algorithms for both Rayleigh-Ritz and orthonormalization (see Section <ref>) steps on various benchmark systems in DFT-FE. These benchmark system are chosen such that the FE discretization errors are ∼ 10^-4 Ha/atom in ground-state energy and ∼ 10^-4 Ha/Bohr in ionic forces. The results in Table <ref> 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 Ĥ, 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 Ĥ as shown in Fig. <ref>. The implementation of the blocked approach used here is similar to the implementation of the blocked approach in the overlap matrix computation (see Section <ref>). 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. <ref> and  <ref>.





 §.§ 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 ρ=F[ρ],
where F[ρ] involves computing the occupied eigenspace for a given ρ. This fixed point iteration can be accelerated by mixing the electron density using an appropriate mixing scheme <cit.>. In the present DFT-FE software release, we implement two kinds of mixing schemes: n-stage Anderson mixing scheme <cit.> and Broyden mixing scheme <cit.>. 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 <cit.>, which provide improved SCF convergence rate independent of system size. 


 §.§ 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 <cit.>. 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 <ref>. 

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 <ref>), Rayleigh–Ritz procedure (steps 1 and 3 of Algorithm <ref>), and electron density computation. In particular, computations over the total number of wavefunctions (Ψ)  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 Ψ_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 Ψ_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 S_b_i=Ψ^†Ψ_b_i and Ψ_b_i^ o=Ψ L^-1_b_i^†, respectively. We note that L^-1_b_i^† denotes N × N_b_i sub-matrix of L^-1^†. Similarly, band parallelized computation of projected Hamiltonian in Rayleigh–Ritz procedure is performed as Ĥ_b_i=Ψ^ o^†HΨ_b_i^ o. We note that all-to-all communications of Ψ_b_i's and Ψ^ o_b_i's across band groups are performed before beginning the orthonormalization and Rayleigh-Ritz procedure, respectively. Furthermore, all-to-all communications of S_b_i's and Ĥ_b_i's across band groups are also performed to compute the matrices S and Ĥ (N × 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 ∼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. <ref>. 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. <ref>, 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.
 
 




 §.§ 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 <ref>. In DFT-FE, geometry optimization is performed using Polak–Ribiere-Polyak non-linear conjugate gradient algorithm (PRP-CG) <cit.> 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 <cit.>. We also plan to implement alternate geometry optimization algorithms like Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) <cit.> and Fast Inertial Relaxation Engine (FIRE) <cit.>, which demonstrate faster convergence than PRP-CG in many cases.


§ 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 <cit.> exchange correlation of the PBE form <cit.> is employed in all the calculations, and additionally ONCV <cit.> pseudopotentials from the SG15 database <cit.> 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 <cit.> 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) <cit.>, and ABINIT <cit.>—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 ∼ 10^-4 Ha and ∼ 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 2 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  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.


 §.§ 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  <cit.>, which uses linearized augmented plane-wave (LAPW) method, and NWChem <cit.> using Gaussian basis set.



  §.§.§ 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 ∼ 10^-4 Ha/atom in ground-state energy, ∼ 10^-4 Ha/Bohr in ionic forces, and ∼ 5 × 10^-6 Ha/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×2×2 denoted by Mg_2x2x2 containing 31 atoms (310 electrons) and 4×4×4 denoted by Mg_4x4x4 containing 255 atoms (2,550 electrons). The relevant mesh parameters for DFT-FE and cut-off energies for QE are shown in Table <ref>. Table <ref> 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 𝒪(10^-5) Ha/atom, ionic forces is 𝒪(10^-4) Ha/Bohr, and hydrostatic stress is 𝒪(10^-6) Hartree/Bohr^3. Similarly, at high accuracy level the agreement in ground-state energy is 𝒪(10^-5) Ha/atom, ionic forces is 𝒪(10^-5) Ha/Bohr, and hydrostatic stress is 𝒪(10^-7) Hartree/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×2×2 denoted by Mo_2x2x2 containing 15 atoms (210 electrons) and 4×4×4 denoted by Mo_4x4x4 containing 127 atoms (1,778 electrons). Table <ref> 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_ 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 <cit.>. 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: Cu_3-Shell containing 147 atoms (2,793 electrons), Cu_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 <ref> shows the comparison between DFT-FE and QE, which demonstrates excellent agreement between the two codes in a non-periodic setting. 



Overall, from Table <ref>, 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. 




  §.§.§ 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) C_60 (Buckminsterfullerene) molecule. The first two benchmark systems, which are periodic, are validated with the  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 Si_ uc^ Gamma, and multiple k-points using a 5 × 5 × 5 Monkhorst-Pack Grid denoted by Si_ uc^ Multi-kpt. The mesh parameters for DFT-FE and the relevant parameters for  are shown in Table <ref>. In the case of ,   in Table <ref> 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 <ref> shows the comparison between DFT-FE and  for both Si_ uc^ Gamma and  Si_ uc^ Multi-kpt, where the agreement in ground-state energy is ∼ 1 × 10^-4 Ha/atom and hydrostatic stress is ∼ 5 × 10^-6 Hartree/Bohr^3. 



Next, we consider a mono-vacancy in a 2 × 2 × 2 supercell constructed from diamond cubic Si unit cell with lattice constant 10.0065225 Bohr. This benchmark system is denoted as SiVac_2×2×2 and contains 63 atoms and 882 electrons. Table <ref> shows the comparison between DFT-FE and  where the agreement in ground-state energy is ∼ 1 × 10^-4 Ha/atom, ionic forces is ∼ 5 × 10^-4 Ha/Bohr, and hydrostatic stress is ∼ 5 × 10^-6 Hartree/Bohr^3. 

Finally, we consider the non-periodic 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 <ref>. The results demonstrate an excellent agreement between DFT-FE and NWChem in the ground-state energy.





Overall, from Table <ref>, we show good agreement between DFT-FE and  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. 





 §.§ 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—Mg_4x4x4 (255 atoms, 2,550 electrons), (b) medium—Mg_8x8x8 (2,047 atoms, 20,470 electrons), and (c) large—Mg_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 Mg_4x4x4, Mg_8x8x8 and Mg_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 (∼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 ∼ 10^-4 Ha per atom and ∼ 10^-4 Ha/Bohr, respectively.

First, we consider the parallel scalability of the smaller system size Mg_4x4x4 as demonstrated in Fig. <ref>. 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: Mg_8x8x8 and Mg_10x10x10, which are shown in Fig. <ref> and  <ref>, respectively.
Based on the comparison of three different parallelization strategies in Section <ref> 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 Mg_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 Mg_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 <ref>). 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. <ref> showing the strong scaling of the various computational steps in DFT-FE for Mg_8x8x8 and Mg_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. <ref> 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 <ref>, that reduce communication costs and latency.






 §.§ 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 <cit.>, and ABINIT v8.8.4 <cit.>. 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 <ref>, 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 <ref>): m=45 for benchmark systems (i) and (ii), and m=50 for benchmark system (iii). Further, N_ fr, which is used in the RR step (see Section <ref>), 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∼10% of N_e/2 for benchmark systems (i) and (ii), and b∼5% of N_e/2 for benchmark system (iii).



  
CPU-time comparisons
In Table <ref>, we compare the average computational CPU-times per SCF iteration step[Measured 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 <ref> 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 <ref> and  <ref>, where we observe that DFT-FE is more computationally efficient than both QE and ABINIT beyond system sizes of ∼ 3000 electrons (300 atoms)  for hcp Mg supercells, and ∼ 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 <ref>) is also a key factor. Next, we consider the non-periodic benchmark problem: Icosahedron nano-particles of varying sizes, in Table <ref>. 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 <ref>, we observe that DFT-FE's efficiency gains over QE increases with increasing system size, achieving efficiency gains of 5.7×, 12.4×, and 11.9× for Mg_8x8x8, Mo_10x10x10, and Cu_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 <ref>. 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—𝒪(N_e^2.12) up to N_e=39,990 for hcp Mg super cells, 𝒪(N_e^2.32) up to N_e=27,986 for bcc Mo super cells, and 𝒪(N_e^2.04) up to N_e=17,537 for Cu nano-particles.





  
Wall-time comparisons
Next, in Table <ref> 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×, 16.1× and 6.9× for Mg_8x8x8, Mo_10x10x10, and Cu_5-shell, respectively. Even at the smallest system sizes DFT-FE is still significantly faster than QE, with speedups of 1.5×, 7.5× and 3.3× for Mg_4x4x4, Mo_6x6x6, and Cu_3-shell, respectively. Additionally, in Table <ref>, we also report some very large scale simulations conducted using DFT-FE: Mg_10x10x10, Mo_13x13x13, and Cu_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 <ref> 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.










§ DEMONSTRATION OF DFT-FE'S CAPABILITIES




 §.§ Geometry optimization



We demonstrate here the capability of DFT-FE to conduct ionic relaxation using the methodology described in Section  <ref>, 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 <cit.>. 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 <ref>. In both the DFT-FE and QE simulations, the discretization errors in ground-state energy and forces are converged to ∼ 10^-4 Ha per atom and ∼ 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 <ref>, 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 × 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 <ref>. The maximum difference between relaxed and the initial un-relaxed TBR coordinates obtained from DFT-FE after removing the rigid body modes is Δ_ maxR^ relax=0.276 Bohr. Finally, using a similar metric we compare the relaxed TBR coordinates obtained from DFT-FE and QE, where we obtain Δ_ maxR^ 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 Δ_ maxR^ 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.






 §.§ Molecular dynamics

We demonstrate the capability of DFT-FE to conduct Born-Oppenheimer molecular dynamics (MD) calculations <cit.>. In particular, we choose fcc Al 3 × 3 × 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. <ref> 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 × 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 × 10^-9 Ha/atom-fs.





 §.§ 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. <ref> (a) along with two high-symmetry lines. ONCV pseudopotential is used along with PBE <cit.> functionals to carry out the DFT calculation. The BZ is sampled by a 8× 8× 8 Monkhorst-Pack grid. The KS eigenvalues along the two high-symmetry lines are extracted from the converged Hamiltonian and are shown in Fig. <ref> (b). For the sake of comparison, a similar calculation is done using QE, and the corresponding bandstructure is overlaid in Fig. <ref> (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× 4× 4 supercell of diamond-Si with a silicon vacancy (VSi) as illustrated in Fig. <ref> (c). This gives rise to spin-polarized midgap states shown on Fig. <ref> (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.   








§ 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^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 ∼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  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× 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× speedups for smaller system sizes (2,550-6,034 electrons) to 7-16× 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) <cit.>. Implementation of enriched finite-element basis <cit.> 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. 

















*
elsarticle-num

















100
url<#>1urlprefixURL href#1#2#2 #1#1

kohn65
W. Kohn, L. J. Sham, Self-consistent equations including exchange and
  correlation effects, Phys. Rev. 140 (4A) (1965) A1133.

kohn96
W. Kohn, Density functional and density matrix method scaling linearly with the
  number of atoms, Phys. Rev. Lett. 76 (1996) 3168–3171.

XCReview2005
G. E. Scuseria, V. N. Staroverov, Progress in the development of
  exchange-correlation functionals (2005).

rodney2017
D. Rodney, L. Ventelon, E. Clouet, L. Pizzagalli, F. Willaime, Ab initio
  modeling of dislocation core properties in metals and semiconductors, Acta
  Mater. 124 (2017) 633–659.

Arias2000
S. Ismail-Beigi, T. A. Arias, Ab initio study of screw dislocations in mo and
  ta: A new picture of plasticity in bcc transition metals, Phys. Rev. Lett. 84
  (2000) 1499–1502.

Trinkle2005
D. R. Trinkle, C. Woodward, The chemistry of deformation: How solutes soften
  pure metals, Science 310 (5754) (2005) 1665–1667.

Trinkle2008
C. Woodward, D. R. Trinkle, L. G. Hector, D. L. Olmsted, Prediction of
  dislocation cores in aluminum from density functional theory, Phys. Rev.
  Lett. 100 (2008) 045507.

Clouet2009
E. Clouet, L. Ventelon, F. Willaime, Dislocation core energies and core fields
  from first principles, Phys. Rev. Lett. 102 (2009) 055502.

shin2011
I. Shin, E. A. Carter, Orbital-free density functional theory simulations of
  dislocations in magnesium, Modelling and Simulation in Materials Science and
  Engineering 20 (1) (2011) 015006.

shin2013
I. Shin, E. A. Carter, Possible origin of the discrepancy in peierls stresses
  of fcc metals: First-principles simulations of dislocation mobility in
  aluminum, Phys. Rev. B 88 (2013) 064106.

iyer2015
M. Iyer, B. Radhakrishnan, V. Gavini, Electronic-structure study of an edge
  dislocation in Aluminum and the role of macroscopic deformations on its
  energetics, J. Mech. Phys. Solids 76 (2015) 260–275.

radhakrishnan2016
B. Radhakrishnan, V. Gavini, Orbital-free density functional theory study of
  the energetics of vacancy clustering and prismatic dislocation loop
  nucleation in Aluminium, Philos. Mag. 96 (2016) 2468–2487.

das2017
S. Das, V. Gavini, Electronic structure study of screw dislocation core
  energetics in Aluminum and core energetics informed forces in a dislocation
  aggregate, J. Mech. Phys. Solids 104 (2017) 115–143.

dawson2018
J. A. Dawson, H. Chen, M. S. Islam, Composition screening of Lithium- and
  Sodium-rich anti-perovskites for fast-conducting solid electrolytes, J.
  Phys. Chem. C 122 (42) (2018) 23978–23984.

dive2018
A. Dive, C. Benmore, M. Wilding, S. W. Martin, S. Beckman, S. Banerjee,
  Molecular dynamics modeling of the structure and Na^+-ion
  transport in Na_2S + SiS_2
  glassy electrolytes, J. Phys. Chem. B 122 (30) (2018) 7597–7608.

cole2016
D. J. Cole, N. D. M. Hine, Applications of large-scale density functional
  theory in biology, J. Phys.: Condens. Matter 28 (39) (2016) 393001.

qe2009
P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni,
  D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso,
  S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann,
  C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri,
  R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia,
  S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, R. M.
  Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project
  for quantum simulations of materials, J. Phys. Condens. Matter 21 (39) (2009)
  395502.

gonze2002
X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese,
  L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami,
  P. Ghosez, J.-Y. Raty, D. Allan, First-principles computation of material
  properties: the {ABINIT} software project, Comput. Mater. Sci. 25 (3)
  (2002) 478–492.

VASP
G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio
  total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (16)
  (1996) 11169–11186.

exciting2014
A. Gulans, S. Kontur, C. Meisenbichler, D. Nabok, P. Pavone, S. Rigamonti,
  S. Sagmeister, U. Werner, C. Draxl, exciting: a full-potential all-electron
  package implementing density-functional theory and many-body perturbation
  theory, J. Phys.: Condens. Matter 26 (36) (2014) 363202.

Pople
W. J. Hehre, R. F. Stewart, J. A. Pople, Self-consistent molecular-orbital
  methods. i. Use of Gaussian expansions of Slater-type atomic orbitals,
  J. Chem. Phys. 51 (6) (1969) 2657–2664.

jensen2002
F. Jensen, Polarization consistent basis sets. ii. estimating the
  Kohn–Sham basis set limit, J. Chem. Phys. 116 (17) (2002) 7372–7379.

cp2k2014
J. Hutter, M. Iannuzzi, F. Schiffmann, J. VandeVondele, cp2k: atomistic
  simulations of condensed matter systems, Wiley Interdiscip. Rev.: Comput.
  Mol. Sci. 4.

blum2009
V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter,
  M. Scheffler, Ab initio molecular simulations with numeric atom-centered
  orbitals, Comput. Phys. Commun. 180 (11) (2009) 2175–2196.

nwchem
M. Valiev, E. Bylaska, N. Govind, K. Kowalski, T. Straatsma, H. V. Dam,
  D. Wang, J. Nieplocha, E. Apra, T. Windus, W. de Jong, Nwchem: A
  comprehensive and scalable open-source solution for large scale molecular
  simulations, Comput. Phys. Commun. 181 (9) (2010) 1477–1489.

tsuchida1995
E. Tsuchida, M. Tsukada, Electronic-structure calculations based on the
  finite-element method, Phys. Rev. B 52 (1995) 5573–5578.

tsuchida1996
E. Tsuchida, M. Tsukada, Adaptive finite-element method for
  electronic-structure calculations, Phys. Rev. B 54 (1996) 7602–7605.

tsuchida1998
E. Tsuchida, M. Tsukada, Large-scale electronic-structure calculations based on
  the adaptive finite-element method, Journal of the Physical Society of Japan
  67 (11) (1998) 3844–3858.

pask1999
J. E. Pask, B. M. Klein, C. Y. Fong, P. A. Sterne, Real-space local polynomial
  basis for solid-state electronic-structure calculations: A finite-element
  approach, Phys. Rev. B 59 (1999) 12352–12358.

pask2005
J. E. Pask, P. A. Sterne, Finite element methods in ab initio electronic
  structure calculations, Modell. Simul. Mater. Sci. Eng. 13 (3) (2005) R71.

sukumar2009
N. Sukumar, J. E. Pask, Classical and enriched finite element formulations for
  bloch-periodic boundary conditions, Int. J. Numer. Methods Eng. 77 (8) (2009)
  1121–1138.

suryanarayana2010
P. Suryanarayana, V. Gavini, T. Blesgen, K. Bhattacharya, M. Ortiz,
  Non-periodic finite-element formulation of Kohn–Sham density functional
  theory, J. Mech. Phys. Solids 58 (2010) 256–280.

zhou2011
H. Chen, L. He, A. Zhou, Finite element approximations of nonlinear eigenvalue
  problems in quantum physics, Comput. Methods in Appl. Mech. Eng. 200 (21)
  (2011) 1846–1865.

motamarri2013
P. Motamarri, M. Nowak, K. Leiter, J. Knap, V. Gavini, Higher-order adaptive
  finite-element methods for Kohn-Sham density functional theory, J.
  Comput. Phys. 253 (2013) 308–343.

SCHAUER2013644
V. Schauer, C. Linder, All-electron Kohn–Sham density functional theory on
  hierarchic finite element spaces, J. Comput. Phys. 250 (2013) 644–664.

zhou2014
H. Chen, X. Dai, X. Gong, L. He, A. Zhou, Adaptive finite element
  approximations for Kohn–Sham models, Multiscale Model. Simul. 12 (4)
  (2014) 1828–1869.

denis2016
D. Davydov, T. D. Young, P. Steinmann, On the adaptive finite element analysis
  of the Kohn–Sham equations: methods, algorithms, and implementation, Int.
  J. Numer. Methods Eng. 106 (11) (2016) 863–888.

kanungo2017
B. Kanungo, V. Gavini, Large-scale all-electron density functional theory
  calculations using an enriched finite-element basis, Phys. Rev. B 95 (2017)
  035112.

Davydov2018
D. Davydov, T. Heister, M. Kronbichler, P. Steinmann, Matrix-free locally
  adaptive finite element solution of density-functional theory with
  nonorthogonal orbitals and multigrid preconditioning, Physica Status Solidi
  B: Basic Solid State Physics 255 (9).
http://dx.doi.org/10.1002/pssb.201800069
  doi:10.1002/pssb.201800069.

parsec2006
L. Kronik, A. Makmal, M. L. Tiago, M. M. G. Alemany, M. Jain, X. Huang,
  Y. Saad, J. R. Chelikowsky, PARSEC — the pseudopotential algorithm for
  real-space electronic structure calculations: recent advances and novel
  applications to nano-structures, Phys. Status Solidi B 243 (5) (2006)
  1063–1079.

rescu2016
V. Michaud-Rioux, L. Zhang, H. Guo, RESCU: A real space electronic structure
  method, J. Comput. Phys. 307 (2016) 593–613.

sparc2017a
S. Ghosh, P. Suryanarayana, SPARC: Accurate and efficient finite-difference
  formulation and parallel implementation of density functional theory:
  Isolated clusters, Comput. Phys. Commun. 212 (2017) 189–204.

sparc2017b
S. Ghosh, P. Suryanarayana, SPARC: Accurate and efficient finite-difference
  formulation and parallel implementation of density functional theory:
  Extended systems, Comput. Phys. Commun. 216 (2017) 109–125.

octopus2015
X. Andrade, D. Strubbe, U. De Giovannini, A. H. Larsen, M. J. T. Oliveira,
  J. Alberdi-Rodriguez, A. Varas, I. Theophilou, N. Helbig, M. J. Verstraete,
  L. Stella, F. Nogueira, A. Aspuru-Guzik, A. Castro, M. A. L. Marques,
  A. Rubio, Real-space grids and the Octopus code as tools for the
  development of new simulation approaches for electronic systems, Phys. Chem.
  Chem. Phys. 17 (2015) 31371–31396.

gpaw
J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen, M. Dułak, L. Ferrighi,
  J. Gavnholt, C. Glinsvad, V. Haikola, H. A. Hansen, H. H. Kristoffersen,
  M. Kuisma, A. H. Larsen, L. Lehtovaara, M. Ljungberg, O. Lopez-Acevedo, P. G.
  Moses, J. Ojanen, T. Olsen, V. Petzold, N. A. Romero, J. Stausholm-Møller,
  M. Strange, G. A. Tritsaris, M. Vanin, M. Walter, B. Hammer, H. HÃ¤kkinen,
  G. K. H. Madsen, R. M. Nieminen, J. K. Nørskov, M. Puska, T. T. Rantala,
  J. Schiøtz, K. S. Thygesen, K. W. Jacobsen, Electronic structure
  calculations with GPAW: a real-space implementation of the projector
  augmented-wave method, Journal of Physics: Condensed Matter 22 (25) (2010)
  253202.

luigi2008
L. Genovese, A. Neelov, S. Goedecker, T. Deutsch, S. A. Ghasemi, A. Willand,
  D. Caliste, O. Zilberberg, M. Rayson, A. Bergman, R. Schneider, Daubechies
  wavelets as a basis set for density functional pseudopotential calculations,
  J. Chem. Phys. 129 (1) (2008) 014109.

skylaris2005
C.-K. Skylaris, P. D. Haynes, A. A. Mostofi, M. C. Payne, Introducing ONETEP:
  Linear-scaling density functional simulations on parallel computers, J. Chem.
  Phys. 122 (8) (2005) 084119.

dgdft2015
W. Hu, L. Lin, C. Yang, DGDFT: A massively parallel method for large scale
  density functional theory calculations, J. Chem. Phys. 143 (12) (2015)
  124110.

motamarri2015
P. Motamarri, V. Gavini, Tucker-tensor algorithm for large-scale Kohn-Sham
  density functional theory calculations, Phys. Rev. B 93 (2017) 035111.

godecker99
S. Goedecker, Linear scaling electronic structure methods, Rev. Mod. Phys. 71
  (1999) 1085–1123.

bowler2012
D. R. Bowler, T. Miyazaki, 𝒪(𝒩) methods in electronic structure
  calculations, Rep. Prog. Phys. 75 (3) (2012) 036503.

fattebert2006
J.-L. Fattebert, F. Gygi, Linear-scaling first-principles molecular dynamics
  with plane-waves accuracy, Phys. Rev. B 73 (2006) 115124.

ls3df2008
L.-W. Wang, Z. Zhao, J. Meza, Linear-scaling three-dimensional fragment method
  for large-scale electronic structure calculations, Phys. Rev. B 77 (2008)
  165113.

motamarri2014
P. Motamarri, V. Gavini, Subquadratic-scaling subspace projection method for
  large-scale Kohn-Sham density functional theory calculations using
  spectral finite-element discretization, Phys. Rev. B 90 (2014) 115127.

skylaris2018
J. Aarons, C.-K. Skylaris, Electronic annealing Fermi operator expansion for
  DFT calculations on metallic systems, J. Chem. Phys. 148 (7) (2018) 074107.

mohr2018
S. Mohr, M. Eixarch, M. Amsler, M. J. Mantsinen, L. Genovese, Linear scaling
  DFT calculations for large tungsten systems using an optimized local basis,
  Nucl. Mater. Energy 15 (2018) 64 – 70.

saad2006
Y. Zhou, Y. Saad, M. L. Tiago, J. R. Chelikowsky, Parallel
  self-consistent-field calculations via Chebyshev-filtered subspace
  acceleration, Phys. Rev. E 74 (6) (2006) 066704.

saad2012
G. Schofield, J. R. Chelikowsky, Y. Saad, A spectrum slicing method for the
  Kohn–Sham problem, Comp. Phys. Comm 183 (3) (2012) 497–505.

zhou2006
Y. Zhou, Y. Saad, M. L. Tiago, J. R. Chelikowsky, Self-consistent-field
  calculations using Chebyshev-filtered subspace iteration, J. Comput. Phys.
  219 (1) (2006) 172–184.

motamarri2018
P. Motamarri, V. Gavini, Configurational forces in electronic structure
  calculations using Kohn-Sham density functional theory, Phys. Rev. B. 97
  (2018) 165132.

qe2017
P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra,
  R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo,
  A. D. Corso, S. de Gironcoli, P. Delugas, R. A. D. Jr, A. Ferretti,
  A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino,
  T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. KÃ¼Ã§Ã¼kbenli,
  M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen,
  A. O. de-la Roza, L. Paulatto, S. PoncÃ©, D. Rocca, R. Sabatini, B. Santra,
  M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari,
  N. Vast, X. Wu, S. Baroni, Advanced capabilities for materials modelling with
  QUANTUM ESPRESSO, J. Phys. Condens. Matter 29 (46) (2017) 465901.

rmartin
R. M. Martin, Electronic structure: basic theory and practical methods,
  Cambridge university press, Cambridge, UK, 2004.

gga1
D. C. Langreth, M. J. Mehl, Beyond the local-density approximation in
  calculations of ground-state electronic properties, Phys. Rev. B 28 (1983)
  1809–1834.

becke
A. D. Becke, Density-functional exchange-energy approximation with correct
  asymptotic behavior, Phys. Rev. A. 38 (1988) 3098–3100.

pw
J. P. Perdew, Y. Wang, Accurate and simple analytic representation of the
  electron-gas correlation energy, Phys. Rev. B. 45 (1992) 13244–13249.

pbe
J. P. Perdew, K. Burke, M. Ernzerhof, Generalized Gradient Approximation made
  simple, Phys. Rev. Lett. 77 (1996) 3865–3868.

kb82
L. Kleinman, D. M. Bylander, Efficacious form for model pseudopotentials, Phy.
  Rev. Lett. 48 (20) (1982) 1425.

tm91
N. Troullier, J. L. Martins, Efficient pseudopotentials for plane-wave
  calculations, Phys. Rev. B 43 (1991) 1993–2006.

oncv2013
D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials, Phys.
  Rev. B 88 (1995) 239906.

science2016
K. Lejaeghere, et al., Reproducibility in density functional theory
  calculations of solids, Sci. 351.

gavini2007
V. Gavini, J. Knap, K. Bhattacharya, M. Ortiz, Non-periodic finite-element
  formulation of orbital free density functional theory, J. Mech. Phys. Solids
  55 (2007) 669–696.

das2015
S. Das, M. Iyer, V. Gavini, Real-space formulation of orbital-free density
  functional theory using finite-element discretization: The case for Al,
  Mg, and Al-Mg intermetallics, Phys. Rev. B 92 (2015) 014104.

mermin
N. W. Ashcroft, N. D. Mermin, Solid State Physics, Hartcourt College
  Publishers, San Diego, US, 1976.

mpgrid
H. J. Monkhorst, J. D. Pack, Special points for brillouin-zone integrations,
  Phys. Rev. B 13 (1976) 5188–5192.

brenner2002
S. C. Brenner, L. R. Scott, The Mathematical Theory of Finite-element Methods,
  Springer, New York, 2002.

dealII90
G. Alzetta, D. Arndt, W. Bangerth, V. Boddu, B. Brands, D. Davydov,
  R. Gassmöller, T. Heister, L. Heltai, K. Kormann, M. Kronbichler,
  M. Maier, J.-P. Pelteret, B. Turcksin, D. Wells, The 
  Library, Version 9.0, Journal of Numerical Mathematics 26 (4) (2018)
  173–183.

bylaska
E. J. Bylaska, M. Holst, J. H. Weare, Adaptive finite element method for
  solving the exact Kohn-Sham equation of density functional theory, J.
  Chem. Theory. Comput. 5 (4) (2009) 937–948.

lehtovaara
L. Lehtovaara, V. Havu, M. Puska, All-electron density functional theory and
  time-dependent density functional theory with high-order finite elements, J.
  Chem. Physics 131 (5) (2009) 054103.

qcofdft
V. Gavini, K. Bhattacharya, M. Ortiz, Quasi-continuum orbital-free
  density-functional theory: A route to multi-million atom non-periodic DFT
  calculation, J. Mech. Phys. Solids 55 (4) (2007) 697–718.

choly2005
N. Choly, G. Lu, W. E, E. Kaxiras, Multiscale simulations in simple metals: A
  density-functional-based methodology, Phys. Rev. B 71 (2005) 094101.

gang2006
G. Lu, E. B. Tadmor, E. Kaxiras, From electrons to finite elements: A
  concurrent multiscale approach for metals, Phys. Rev. B 73 (2006) 024108.

p4est2011
C. Burstedde, L. C. Wilcox, O. Ghattas, : Scalable algorithms
  for parallel adaptive mesh refinement on forests of octrees, SIAM J. Sci.
  Comput 33 (3) (2011) 1103–1133.

zhou2008
X. Dai, J. Xu, A. Zhou, Convergence and optimal complexity of adaptive finite
  element eigenvalue computations, Numerische Mathematik 110 (3) (2008)
  313–355.

chen2011
H. Chen, L. He, A. Zhou, Finite element approximations of nonlinear eigenvalue
  problems in quantum physics, Computer Methods in Applied Mechanics and
  Engineering 200 (21) (2011) 1846 – 1865.

bao2012
G. Bao, G. Hu, D. Liu, An h-adaptive finite element solver for the calculations
  of the electronic structures, Journal of Computational Physics 231 (14)
  (2012) 4967 – 4979.

chen2014
H. Chen, X. Dai, X. Gong, L. He, A. Zhou, Adaptive finite element
  approximations for kohn-sham models, Multiscale Modeling and Simulation
  12 (4) (2014) 1828–1869.

shen2018
Y. Shen, Y. Kuang, G. Hu, An asymptotics-based adaptive finite element method
  for kohn–sham equation, Journal of Scientific Computing.

radio
R. Radovitzky, M. Ortiz, Error estimation and adaptive meshing in strongly
  nonlinear dynamic problems, Comput. Methods Appl. Mech. Engrg. 172 (1) (1999)
  203–240.

ciarlet2002
P. G. Ciarlet, The Finite Element Method for Elliptic Problems, SIAM,
  Philadelphia, 2002.

anderson1965
D. G. Anderson, Iterative procedures for nonlinear integral equations, J.
  Assoc. Comput. Mach. 12 (4) (1965) 547–560.

broyden1965
C. G. Broyden, A class of methods for solving nonlinear simultaneous equations,
  Math. Comput. (1965) 577–593.

kronbichler2012
M. Kronbichler, K. Kormann, A generic interface for parallel cell-based finite
  element operator application, Comput. Fluids 63 (2012) 135–147.

bekas2010
C. Bekas, A. Curioni, Very large scale wavefunction orthogonalization in
  density functional theory electronic structure calculations, Comput. Phys.
  Commun. 181 (6) (2010) 1057–1068.

scalapack1997
L. S. Blackford, J. Choi, A. Cleary, E. D'Azevedo, J. Demmel, I. Dhillon,
  J. Dongarra, S. Hammarling, G. Henry, A. Petitet, K. Stanley, D. Walker,
  R. C. Whaley, ScaLAPACK Users' Guide, SIAM, Philadelphia, PA, 1997.

tsuchida2012
E. Tsuchida, Y.-K. Choe, Iterative diagonalization of symmetric matrices in
  mixed precision and its application to electronic structure calculations,
  Comput. Phys. Commun. 183 (4) (2012) 980–985.

amartya2018
A. S. Banerjee, L. Lin, P. Suryanarayana, C. Yang, J. E. Pask, Two-level
  chebyshev filter based complementary subspace method: Pushing the envelope of
  large-scale electronic structure calculations, J. Chem. Theory Comput. 14 (6)
  (2018) 2930–2946.

motamarri2017
P. Motamarri, V. Gavini, K. Bhattacharya, M. Ortiz, Spectrum-splitting approach
  for Fermi-operator expansion in all-electron Kohn-Sham DFT
  calculations, Phys. Rev. B 95 (2017) 035111.

elpa2014
A. Marek, V. Blum, R. Johanni, V. Havu, B. Lang, T. Auckenthaler, A. Heinecke,
  H.-J. Bungartz, H. Lederer, The ELPA library: scalable parallel eigenvalue
  solutions for electronic structure theory and computational science, J.
  Phys.: Condens. Matter 26 (21) (2014) 213201.

elpaCray
B. Cook, T. Kurth, J. Deslippe, P. Carrier, N. Hill, N. Wichmann, Eigensolver
  performance comparison on Cray XC systems, Concurr. Comp.-Pract. E. 0 (0)
  (2018) e4997.

elpaOpt
P. Kus, A. Marek, S. S. Koecher, H. Kowalski, C. Carbogno, C. Scheurer,
  K. Reuter, M. Scheffler, H. Lederer, Optimizations of the eigensolvers in the
  ELPA library, arXiv:1811.01277.

Kerker1981
G. P. Kerker, Efficient iteration scheme for self-consistent pseudopotential
  calculations, Phys. Rev. B 23 (1981) 3082–3084.

eyert1996
V. Eyert, A comparative study on methods for convergence acceleration of
  iterative vector sequences, J. Comp. Phys 124 (2) (1996) 271–285.

kudin2002
K. N. Kudin, G. E. Scuseria, E. Cancès, A black-box self-consistent
  field convergence algorithm: One step closer, J. Chem. Phys. 116 (19) (2002)
  8255–8261.

linlin2013b
L. Lin, C. Yang, Elliptic preconditioner for accelerating the self-consistent
  field iteration in Kohn–Sham density functional theory, SIAM J. Sci.
  Comput 35 (5) (2013) S277–S298.

zhou2018
Y. Zhou, H. Wang, Y. Liu, X. Gao, H. Song, Applicability of Kerker
  preconditioning scheme to the self-consistent density functional theory
  calculations of inhomogeneous systems, Phys. Rev. E 97 (2018) 033305.

BangerthBursteddeHeisterEtAl11
W. Bangerth, C. Burstedde, T. Heister, M. Kronbichler, Algorithms and data
  structures for massively parallel generic adaptive finite element codes, ACM
  Trans. Math. Software 38 (2) (2011) 14:1–14:28.

shewchuk1994
J. R. Shewchuk, An introduction to the conjugate gradient method without the
  agonizing pain, Tech. rep., Carnegie Mellon University, Pittsburgh, PA, USA
  (1994).

brent1971
R. P. Brent, An algorithm with guaranteed convergence for finding a zero of a
  function, Comput. J. 14 (4) (1971) 422–425.

liu1989
D. C. Liu, J. Nocedal, On the limited memory BFGS method for large scale
  optimization, Math. Program. 45 (1-3) (1989) 503–528.

bitzek2006
E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, P. Gumbsch, Structural
  relaxation made simple, Phys. Rev. Lett. 97 (2006) 170201.

oncv2015
M. Schlipf, F. Gygi, Optimization algorithm for the generation of ONCV
  pseudopotentials, Comput. Phys. Commun. 196 (2015) 36–44.

cuNano
J. M. Rahm, P. Erhart, Beyond magic numbers: Atomic scale equilibrium
  nanoparticle shapes for any size, Nano Lett. 17 (9) (2017) 5775–5781.

campagna2007
S. Campagna, F. Puntoriero, F. Nastasi, G. Bergamini, V. Balzani,
  Photochemistry and photophysics of coordination compounds: Ruthenium, Topics
  in current chemistry (2007) 117–214.

marx_hutter_2009
D. Marx, J. Hutter, Ab Initio Molecular Dynamics: Basic Theory and Advanced
  Methods, Cambridge University Press, 2009.

Ghosh2019
K. Ghosh, H. Ma, V. Gavini, G. Galli, All-electron density functional
  calculations for electron and nuclear spin interactions in molecules and
  solids, arXiv:1902.07377.

vb90
D. Vanderbilt, Soft self-consistent pseudopotentials in a generalized
  eigenvalue formalism, Phys. Rev. B 41 (1990) 7892–7895.



