SQDFT: Spectral Quadrature method for large-scale parallel $\mathcal{O}(N)$ Kohn-Sham calculations at high temperature
Phanish Suryanarayana Phanisri P. Pratapa Abhiraj Sharma John E. Pask






gatech]Phanish Suryanarayanacor
gatech]Phanisri P. Pratapa
gatech]Abhiraj Sharma
llnl]John E. Pask
[gatech]College of Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA
[llnl]Physics Division, Lawrence Livermore National Laboratory, Livermore, CA 94550, USA
[cor]Corresponding Author (phanish.suryanarayana@ce.gatech.edu) 



We present SQDFT: a large-scale parallel implementation of the Spectral Quadrature (SQ) method for 𝒪(N) Kohn-Sham Density Functional Theory (DFT) calculations at high temperature.

Specifically, we develop an efficient and scalable finite-difference implementation of the infinite-cell Clenshaw-Curtis SQ approach, in which results for the infinite crystal are obtained by expressing quantities of interest as bilinear forms or sums of bilinear forms, that are then approximated by spatially localized Clenshaw-Curtis quadrature rules. We demonstrate the accuracy of SQDFT by showing systematic convergence of energies and atomic forces with respect to SQ parameters to reference diagonalization results, and convergence with discretization to established planewave results, for both metallic and insulating systems. We further demonstrate that SQDFT achieves excellent strong and weak parallel scaling on computer systems consisting of tens of thousands of processors, with near perfect 𝒪(N) scaling with system size and wall times as low as a few seconds per self-consistent field iteration. Finally, we verify the accuracy of SQDFT in large-scale quantum molecular dynamics simulations of aluminum at high temperature.



Electronic structure, Linear scaling, Metallic systems, High temperature, Quantum molecular dynamics, High performance computing, Parallel computing








§ INTRODUCTION
  
Kohn-Sham Density Functional Theory (DFT) <cit.> is a powerful tool for predicting and understanding a wide range of materials properties, from the first principles of quantum mechanics, with no empirical or adjustable parameters. The tremendous popularity of DFT is a consequence of its high accuracy to cost ratio relative to other such ab initio theories. However, the solution of the Schrödinger type eigenproblem for the Kohn-Sham orbitals remains a challenging task. In particular, since the orbitals need to be orthogonal and increase in number linearly with the number of atoms N, the overall computational complexity of DFT calculations scales as 𝒪(N^3) and the memory requirement scales as 𝒪(N^2) (see, e.g., <cit.>).

The orthogonality constraint on the orbitals also results in global communications between processors in parallel computing, which limits parallel scalability. The need for high performance parallel computing is especially crucial for quantum molecular dynamics (QMD) calculations <cit.>, wherein tens or hundreds of thousands of Kohn-Sham solutions can be required to complete a single simulation. 

In order to overcome the critical 𝒪(N^3) scaling bottleneck, much research in the past two decades has been devoted to the development of 𝒪(N) solution strategies (see, e.g., <cit.> and references therein). Rather than calculate the orthonormal Kohn-Sham orbitals, these techniques directly determine the electron density, energy, and atomic forces in 𝒪(N) operations by exploiting the decay of the density matrix <cit.>.[The real-space density matrix has exponential decay for insulating systems as well as metallic systems at finite temperature <cit.>.] These efforts have yielded significant advances, culminating in mature implementations of a number of approaches <cit.>. However, significant challenges remain. In particular, the accuracy and stability of 𝒪(N) methods remain ongoing concerns due to the need for additional computational parameters, subtleties in determining sufficient numbers and/or centers of localized orbitals, limitations of underlying basis sets, and calculation of accurate atomic forces, as required for structural relaxation and molecular dynamics simulations <cit.>. In addition, efficient large-scale parallelization poses a significant challenge due to complex communications patterns and load balancing issues. Finally, and perhaps most importantly, the assumption of a band gap in the electronic structure makes existing methods inapplicable to metallic systems <cit.>.[This is also the case for insulating systems at sufficiently high temperature which results in conduction bands becoming partially occupied.] 






High-temperature DFT calculations present additional challenges <cit.>. 
Such calculations have a number of applications, including the study of warm dense matter and dense plasmas, as occur in laser experiments, and the interiors of giant planets and stars <cit.>. 



Particular challenges include the need for a significantly larger number of orbitals to be computed, as the number of partially occupied states increases, and need for more diffuse orbitals, as higher-energy states become less localized. Consequently, 𝒪(N^3) methods as well as local-orbital based 𝒪(N) methods have very large prefactors, which makes high-temperature QMD calculations for even small systems intractable. Recent work to address these challenges includes orbital-free molecular dynamics (OFMD) <cit.> wherein the standard Kohn-Sham kinetic energy is replaced by an approximation in terms of the density, extended first principles molecular dynamics (ext-FPMD) <cit.> wherein higher-energy states are approximated as planewaves rather than computed explicitly, and finite-temperature potential functional theory (PFT) <cit.> wherein an orbital-free free energy approximation is constructed through a coupling-constant formalism. While OFMD can miss electronic shell structure effects <cit.>, ext-FPMD and PFT have been shown to capture such effects in initial applications.

The recently developed Spectral Quadrature (SQ) method for 𝒪(N) Kohn-Sham calculations <cit.> addresses both scaling with number of atoms and scaling with temperature, while retaining systematic convergence to standard 𝒪(N^3) results for metals and insulators alike. 
In this approach, all quantities of interest are expressed as bilinear forms or sums of bilinear forms, which are then approximated by Clenshaw-Curtis quadrature rules that remain spatially localized by exploiting the locality of electronic interactions in real-space <cit.>, i.e, the exponential decay of the density matrix in real-space for insulators as well as metals at finite temperature. In conjunction with local reformulation of the electrostatics, this technique enables the 𝒪(N) evaluation of the electronic density, energy, and atomic forces. The computational cost of SQ decreases rapidly with increasing temperature due to the enhanced locality of the electronic interactions and the increased smoothness of the Fermi-Dirac function. Further, it is well suited to scalable high-performance parallel computing since a majority of the communication is localized to nearby processors, whose pattern remains fixed throughout the simulation. The SQ approach also permits infinite-crystal calculations without recourse to Brillouin zone integration or large supercells, a technique referred to as the infinite-cell method <cit.>.

In this paper, we present SQDFT: a parallel implementation of the SQ method for 𝒪(N) Kohn-Sham DFT calculations at high temperature.[Though we focus on high-temperature calculations in this work, SQDFT is also capable of performing 𝒪(N) DFT calculations at ambient temperature, with a larger prefactor (Appendix <ref>).] Specifically, we develop a finite-difference implementation of the  infinite-cell variant of the Clenshaw-Curtis SQ method that can efficiently scale on large-scale parallel computers. 

We verify the accuracy of SQDFT by showing systematic convergence of energies and atomic forces to reference diagonalization results, and convergence with discretization to established planewave results, for both metallic and insulating systems. We further show that SQDFT achieves excellent strong and weak parallel scaling on computer systems consisting of tens of thousands of cores, with near perfect 𝒪(N) scaling with system size and wall times as low as a few seconds per self-consistent field (SCF) iteration. Finally, we verify the accuracy of SQDFT in large-scale quantum molecular dynamics simulations of aluminum at high temperature.

The remainder of this paper is organized as follows. In Section <ref>, we review the 𝒪(N) density matrix formulation of DFT. Next, we discuss the formulation and implementation of SQDFT in Section <ref> and study its accuracy, efficiency, and scaling in Section <ref>. Finally, we provide concluding remarks in Section <ref>.




§ 𝒪(N) DENSITY FUNCTIONAL THEORY
 
Consider a cuboidal domain Ω containing N atoms, the unit cell of an infinite crystal. Let the nuclei be positioned at = {_1, _2, …, _N } and let there be a total of N_e valence electrons. Neglecting spin and Brillouin zone integration, the nonlinear eigenproblem for the electronic ground state in Kohn-Sham Density Functional Theory (DFT) can be written as <cit.>

    𝒟    =     g(ℋ,μ,σ) = ( 1 + exp(ℋ-μℐ/σ)   )^-1 , 
    ℋ    =     -1/2∇^2 + V_xc + ϕ + 𝒱_nl ,
 
where 𝒟 is the density matrix; g is the Fermi-Dirac function; μ is the Fermi level, which is determined by solving for the constraint on the total number of electrons, 
i.e., 2 Tr(𝒟) = N_e; 

σ = k_B T is the smearing, where k_B is Boltzmann's constant and T is the electronic temperature[The electronic temperature/smearing is typically set to be equal to the ionic temperature in QMD simulations, particularly for those performed at high temperature.]; ℋ is the Hamiltonian; V_xc is the exchange-correlation potential; 𝒱_nl is the nonlocal pseudopotential; and ϕ is the electrostatic potential, the solution to the Poisson equation <cit.>

    -1/4 π∇^2 ϕ(,) = ρ_𝒟() + b(,)
     
subject to periodic boundary conditions. Above, 

    ρ_𝒟() = 2 𝒟(,)

is the electron density and b = ∑_I b_I is the total pseudocharge density of the nuclei <cit.>, where b_I is the pseudocharge density of the I^th nucleus and the index I extends over all atoms in ^3.

Once the electronic ground state has been determined, the free energy can be written as <cit.>[The repulsive energy correction for overlapping pseudocharges <cit.> has been explicitly included in the free energy expression. This term plays a particularly important role in high-temperature simulations since the ions get significantly closer compared to ambient temperature. ] 

    ℱ ()     =     2Tr(𝒟ℋ)  + E_xc(ρ_𝒟) - ∫_Ω V_xc(ρ_𝒟()) ρ_𝒟()  d + 1/2∫_Ω (b(,)-ρ_𝒟()) ϕ(,)  d
        +     1/2∫_Ω( b̃(,) + b(,) ) V_c (,)  d - 1/2∑_I∫_Ωb̃_I(,_I) Ṽ_I(,_I)  d
        +     2 σTr( 𝒟log𝒟 + (ℐ-𝒟) log (ℐ-𝒟) )  ,

where E_xc is the exchange-correlation energy; b̃ = ∑_Ib̃_I is the total reference pseudocharge density,  where b̃_I is the reference pseudocharge density of the I^th nucleus and the summation index I runs over all atoms in ^3; V_c = ∑_I V_c,I = ∑_I(Ṽ_I-V_I ), where Ṽ_I and V_I are the potentials generated by b̃_I and b_I, respectively; and ℐ is the identity operator. The first term is the band structure energy (E_band) and the last term is the energy associated with the electronic entropy (E_ent). 

In order to update the positions of the ions during the Born-Oppenheimer quantum molecular dynamics (QMD) simulations, the Hellmann-Feynman force on the I^th nucleus can be written as
<cit.>[For the reasons mentioned in footnote <ref>, the force corresponding to the repulsive energy correction for overlapping pseudocharges <cit.> has been explicitly included in the atomic force expression. Note that this term has been significantly simplified from its original form by utilizing the relation ∫∇ f(|-_I|) g(|-_I|)  d=0 for spherically symmetric functions f and g, which in the present case are b_I, b_I', V_I, and V_I'.]

    𝐟_I    =    ∑_I'∫_Ω∇ b_I'(,_I') ϕ(,)  d +  1/2∑_I'∫_Ω[ (b̃(,)+b(,) ) ∇ V_c,I'
        +    ( ∇b̃_I'(,_I') + ∇ b_I'(,_I') ) V_c(,)   ]  d - 4 ∑_I'Tr(𝒱_nl^I'∇𝒟)  ,

where the summation index I' runs over the I^th atom and its periodic images, and 𝒱_nl,I' is the nonlocal pseudopotential associated with the I'^th atom. The first two terms together constitute the local component of the force (𝐟_I^l) and the last term is the nonlocal component of the force (𝐟_I^nl). 

Within a real-space representation, the density matrix 𝒟 has exponential decay for insulators as well as metallic systems at finite electronic temperature/smearing <cit.>. This decay in the density matrix is exploited by linear-scaling methods through truncation within the calculations to enable 𝒪(N) computation of the electron density, energy, and atomic forces. In doing so, it has been observed that there is exponential convergence in the energy and forces with the size of the truncation region for insulating <cit.> as well as metallic systems at finite temperature <cit.>. Note that in the above description for DFT, we have employed a local reformulation of the electrostatics to enable 𝒪(N) scaling for the complete Kohn-Sham problem. 



§ FORMULATION AND IMPLEMENTATION OF SQDFT
 

SQDFT is a large-scale parallel implementation of the Clenshaw-Curtis Spectral Quadrature (SQ) method[The Clenshaw-Curtis variant of SQ is chosen here because it is more efficient compared to its Gauss counterpart <cit.>, particularly in the computation of the nonlocal component of the forces <cit.>.] <cit.> for 𝒪(N) Kohn-Sham density functional calculations. 
In this approach, all quantities of interest (in discrete form) are expressed as bilinear forms or sums of bilinear forms, that are then approximated by spatially localized quadrature rules. The method is identically applicable to insulating and metallic systems and well suited to massively parallel computation. Furthermore, the SQ method becomes more efficient as temperature is increased, since electronic interactions become more localized and the representation of the Fermi-Dirac function becomes more compact <cit.>. 

We employ the infinite-cell version of the Clenshaw-Curtis SQ method, wherein the results corresponding to the infinite crystal are obtained without recourse to Brillouin zone integration or large supercells <cit.>. Specifically, rather than employ Bloch boundary conditions for the orbitals on Ω, zero-Dirichlet boundary conditions are prescribed at infinity, and the relevant components  of the density matrix for spatial points within Ω are calculated by utilizing the potential within the truncation region surrounding that point <cit.>. Periodic boundary conditions are retained for the electrostatic potential. Indeed, the infinite-cell SQ approach is equivalent to the standard Γ-point SQ calculation when the size of the truncation region is smaller than the size of the domain, a situation common in large-scale DFT simulations, particularly those at high temperature. 

We utilize a high-order finite-difference discretization in order to exploit the locality of electronic interactions in real space, enable systematic convergence, and facilitate large-scale parallel implementation. We solve the fixed-point problem in Eq. <ref> using the self consistent field (SCF) method <cit.>[In SQDFT, we perform a fixed-point iteration with respect to the effective potential: V_eff = V_xc + ϕ .], whose convergence is accelerated using the Periodic Pulay mixing scheme <cit.>, a technique that significantly outperforms the well-established Anderson/Pulay mixing <cit.>. We solve the Poisson problem in Eq. <ref> using the Alternating Anderson-Richardson (AAR) method <cit.>, an approach that outperforms the conjugate gradient method <cit.> in the context of large-scale parallel computations <cit.>. We perform NVE (microcanonical) simulations using the leapfrog method and NVT (canonical) simulations using the Verlet algorithm with the Nose-Hoover thermostat. At each MD step, we extrapolate the electron density using information from the previous two steps <cit.> to reduce SCF iterations. We parallelize the calculations using domain decomposition, with the communication between processors handled via the Message Passing Interface (MPI) <cit.>.

In Fig. <ref>, we outline the key steps in a QMD simulation. We employ the Clenshaw-Curtis SQ method for calculating the electron density in each SCF iteration as well as for computing the energy and atomic forces (nonlocal component) once the electronic ground state has been determined. In the sections below, 

we discuss the calculation of the electron density, energy, and atomic forces in SQDFT. For additional details on the formulation and implementation of the pseudocharges, AAJ/AAR method, and Periodic Pulay mixing scheme, we refer the reader to the relevant previous works  <cit.>.




 §.§ Finite-difference discretization
 
We consider a cubical domain Ω and discretize it with a uniform grid of spacing h, resulting in N_d = n^3 finite-difference nodes, the collection of which is referred to as K_Ω. We parallelize the calculations by employing domain decomposition, i.e., we partition the domain into cubes of equal size such that Ω = ⋃_p=1^N_PΩ_p, where N_P is the total number of processors and Ω_p denotes the domain local to the p^th processor, each of which contains N_d^p = N_d/N_P finite-difference nodes. We refer to the collection of finite-difference nodes belonging to the p^th processor as K_Ω^p, where K_Ω = ⋃_p=1^N_P K_Ω^p with K_Ω^p ∩ K_Ω^q = ∅ if p ≠ q. 

We approximate the Laplacian arising in the Hamiltonian and the local electrostatic reformulation using the central finite-difference approximation:

    ∇^2_h f |^(i,j,k)   ≈   ∑_d=0^n_o w_d(f^(i+d,j,k) + f^(i-d,j,k) + f^(i,j+d,k) + f^(i,j-d,k) + f^(i,j,k+d) + f^(i,j,k-d))  , 
    
    w_0     =     - 1/h^2∑_q=1^n_o1/q^2 ,  
    
    w_d     =    2 (-1)^d+1/h^2 d^2(n_o!)^2/(n_o-d)! (n_o+d)! ,    d=1, 2, …, n_o  ,

where f^(i,j,k) denotes the value of the function f at the node indexed by (i,j,k) and 2n_o is the order of the approximation. Similarly, we approximate the gradient operator arising in the atomic forces using central finite-differences:

    ∇_h f |^(i,j,k)   ≈   ∑_d=1^n_ow̃_d ( ( f^(i+d,j,k) - f^(i-d,j,k)) 𝐞̂_1 + ( f^(i,j+d,k) - f^(i,j-d,k)) 𝐞̂_2  + ( f^(i,j,k+d) - f^(i,j,k-d)) 𝐞̂_3 )  , 
    w̃_d     =    (-1)^d+1/h d(n_o!)^2/(n_o-d)! (n_o+d)! ,    d=1, 2, …, n_o  ,

where 𝐞̂_1, 𝐞̂_2, and 𝐞̂_3 are the unit vectors along the edges of Ω. We enforce periodic boundary conditions by mapping any index that does not correspond to a node in the finite-difference grid to its periodic image within Ω. We approximate the spatial integrals arising in the Hamiltonian (projectors of the nonlocal pseudopotential), energy, and atomic forces using the trapezoidal rule:

    ∫_Ω f()  d≈ h^3  ∑_i, j, k=1^n f^(i,j,k) .

Even though this is a low order quadrature scheme, it has been chosen to ensure that the discrete free energy is consistent with the discrete Kohn-Sham equations, i.e., the calculated electronic ground state corresponds to the minimum of the free energy within the finite-difference discretization <cit.>. 

In accordance with the nearsightedness principle <cit.>, we define the region of influence for any finite-difference node as the cube of side 2R_cut centered at that node. A cube is chosen rather than a sphere due to its simplicity and efficiency within the finite-difference implementation. The parameter R_cut corresponds to the distance beyond which electronic interactions are ignored, i.e., theoretically R_cut corresponds to the truncation radius for the density matrix. Within the above described finite-difference framework, we define the nodal Hamiltonian _q ∈^N_s × N_s for any node q ∈ K_Ω as the restriction of the Hamiltonian to its region of influence <cit.>, where N_s = (2R_cut/h + 1)^3 is the number of finite-difference nodes within the region of influence. Similarly, 𝐰_q ∈^N_s × 1, ∇_h,q∈^N_s × N_s, and 𝐕_nl,q^I ∈^N_s × N_s represent the restriction of the standard basis vector, gradient matrix, and nonlocal pseudopotential matrix of the I^th atom to the region of influence, respectively. It is important to note that _q, ∇_h,q, and 𝐕_nl,q^I are not explicitly determined/stored in SQDFT, rather their multiplication with a vector is directly computed in a matrix-free way.




 §.§ Electron density
 
In each iteration of the SCF method, the electron density (Eq. <ref>) needs to be computed at all the finite-difference nodes. In SQDFT, the electron density at the q ∈ K_Ω^p node in the p^th processor is calculated using the relations <cit.>

    ρ_q     =    2/h^3∑_j=0^n_pl c_q^j ρ_q^j  , 
    
    c_q^j     =    2/π∫_-1^1g(r,μ̂_q,σ̂_q) T_j(r)/√(1-r^2) dr ,    j=0, 1, …, n_pl , 
    ρ_q^j     =    𝐰_q^T 𝐭_q^j ,    j=0, 1, …, n_pl/2
    ρ_q^1 - 2(𝐭_q^j+1/2)^T 𝐭_q^j-1/2 ,    j=n_pl/2+1, n_pl/2+3, …, n_pl-1 
    ρ_q^0 - 2(𝐭_q^j/2)^T 𝐭_q^j/2 ,    j=n_pl/2+2, n_pl/2+4, …, n_pl

where n_pl is the order of the Clenshaw-Curtis quadrature, chosen here to be a multiple of four for simplicity; T_j denotes the Chebyshev polynomial of degree j; c_q^j is the coefficient of T_j in the polynomial expansion of the Fermi-Dirac function, with the value of c_q^0 half of that given in the expression; μ̂_q=(μ-χ_q)/ζ_q is the scaled and shifted Fermi energy, where χ_q = (λ_q^max+λ_q^min)/2 and ζ_q = (λ_q^max-λ_q^min)/2, with λ_q^max and λ_q^min denoting the maximum and minimum eigenvalues of _q, respectively; σ̂_q = σ/ζ_q is the scaled smearing; and 𝐭_q^j = T_j(𝐇̂_q) 𝐰_q ∈^N_s × 1, which is determined using the following iteration, obtained as a consequence of the three term recurrence relation of Chebyshev polynomials:


    𝐭^i+1_q    =     2𝐇̂_q 𝐭^i_q-𝐭^i-1_q  ,    i=1, 2, …, n_pl/2
    𝐭_q^1     =    𝐇̂_q 𝐰_q  ,  𝐭_q^0 = 𝐰_q  ,

where _q = (_q-χ_q 𝐈)/ζ_q is the scaled and shifted nodal Hamiltonian whose spectrum lies in the interval [-1,1]. Note that the iteration in Eq. <ref> proceeds only up to n_pl/2 rather than n_pl, since we have employed the product property of Chebyshev polynomials[2 T_j(r) T_k(r) = T_j+k(r) + T_|j-k|(r). Note that in previous work where the expression for the electron density has been derived <cit.>, this property has not been utilized, and therefore the iteration in Eq. <ref> proceeds up to n_pl. The present approach reduces the cost for the computation of ρ_q^j by a factor of two.]
to directly calculate ρ_q^j for j=n_pl/2 +1, n_pl/2+2, …, n_pl, as described by Eq. <ref>. Since the values of ρ_q^j are independent of the Fermi level μ, they are first computed and stored. Next, μ is determined by satisfying the constraint on the total number of electrons:

    2 ∑_p=1^N_P∑_q∈ K_p∑_j=0^n_pl c_q^j ρ_q^j = N_e  ,

where the values of c_q^j are given by Eq. <ref>. Finally, for the c_q^j corresponding to the Fermi level μ, the electron density is calculated via Eq. <ref>. 

Since the simulation domain Ω is cubical and we have employed a uniform finite-difference grid with uniform domain decomposition, the layout of  effective potential values V_eff required from neighboring processors as part of the nodal Hamiltonians for grid points in each processor is identical. To accomplish this, we utilize the MPI command  <cit.> to communicate the required values of V_eff (i.e., those within the region of influence for grid points in each processor) between processors. Doing so reduces the number of MPI related calls that would otherwise be required in every matrix-vector multiplication. After the communication is complete, for every finite-difference node q∈ K_Ω^p, we first calculate λ_q^max and λ_q^min—maximum and minimum eigenvalues of _q, respectively—using the Lanczos method <cit.>. Next, we perform the recursive iteration in Eq. <ref>, the results of which are used to calculate ρ_q^j as given in Eq. <ref>. The matrix-vector multiplications required as part of the Lanczos method and the recursive iteration in Eq. <ref> are performed in a matrix-free manner. The Chebyshev coefficients c_q^j are calculated using the discrete orthogonality of the Chebyshev polynomials <cit.>. The Fermi level is determined using Brent's method <cit.>, with Newton-Raphson's method becoming the preferred choice as the temperature is increased. While doing so, the global communication between processors is handled by using the  command. 

During the electron density calculation, the memory costs are dominated by the storage of three vectors during the recursive iteration in Eq. <ref> (t_q^i+1, t_q^i, and t_q^i-1) and the storage of ρ_q^j for the calculation of the electron density in Eq. <ref>. Therefore, the memory costs per processor scale as 𝒪(3 N_s + n_pl N_d^p). The computational costs are dominated by the the matrix-vector products in the Lanczos iteration and the recursive iteration in Eq. <ref>. Therefore, the computational cost per processor scales as 𝒪(n_lancz N_s N_d^p + 1/2 n_pl N_s N_d^p ), where n_lancz is the number of Lanczos iterations required for determining λ_q^max and λ_q^min. Since N_s, n_pl, and n_lancz remain independent of system size, the overall memory and computational costs scale linearly with the number of finite-difference nodes in the domain Ω, and therefore 𝒪(N) with respect to the number of atoms. 





 §.§ Free energy
 
Once the electronic ground state has been determined, i.e., the SCF iteration has converged, the free energy (Eq. <ref>) needs to be calculated. 



  
Band structure energy The band structure energy in the Clenshaw-Curtis SQ method <cit.> takes the following form in parallel computations 

    E_band    =     2∑_p=1^N_P∑_q∈ K_p∑_j=0^n_pl (χ_q c_q^j + ζ_q d_q^j) ρ_q^j , 
    
    d_q^j     =    2/π∫_-1^1r g(r,μ̂_q,σ̂_q) T_j(r)/√(1-r^2) dr ,
 
where d_q^j is the coefficient of T_j in the polynomial expansion of the band structure energy function (i.e., r g(r)), with the value of d_q^0 half of that given in the expression. In addition, c_q^j and ρ_q^j are as given in Eqs. <ref> and <ref>, respectively. 



  
Electronic entropy energy The electronic entropy energy in the Clenshaw-Curtis SQ approach <cit.> takes the following form in parallel computations 

    S     =      2 σ∑_p=1^N_P∑_q∈ K_p∑_j=0^n_pl e_q^j ρ_q^j , 
    
    e_q^j     =    2/π∫_-1^1g(r,μ̂_q,σ̂_q) log g(r,μ̂_q,σ̂_q) + (1-g(r,μ̂_q,σ̂_q)) log (1-g(r,μ̂_q,σ̂_q)) T_j(r)/√(1-r^2) dr ,

where e_q^j is the coefficient of T_j in the polynomial expansion of the electronic entropy energy function (i.e., g(r) log g(r) + (1-g(r)) log (1-g(r))), with the value of e_q^0 half of that given in the expression. Again, ρ_q^j is as given in Eq. <ref>. 



  
Free energy The free energy of the system in SQDFT is computed as

    ℱ()     =     h^3 ∑_p=1^N_P∑_q∈ K_p( 2/h^3∑_j=0^n_pl (χ_q c_q^j + ζ_q d_q^j) ρ_q^j +  ε_xc(ρ_q) ρ_q -   V_xc(ρ_p) ρ_q + 1/2  (b_q - ρ_q) ϕ_q  
            + 1/2  (b̃_q + b_q) V_c,q  - 1/2∑_I ∈ D_p^bb̃_I,qṼ_I,q + 2 σ/h^3∑_j=0^n_pl e_q^j ρ_q^j)  ,

where the spatial integrals in Eq. <ref> have been approximated using the trapezoidal rule in Eq. <ref> and D_p^b is the set of all atoms (considering all atoms in ^3) whose pseudocharges have overlap with the processor domain Ω_p. We note that the exchange correlation energy E_xc has been modeled using the Local Density Approximation (LDA) <cit.>, wherein ϵ_xc(ρ) is the sum of the exchange and correlation energy per particle of a uniform electron gas of density ρ. 

During the free energy calculation, the values of c_q^j and ρ^j_q determined as part of the electron density computation in the last SCF iteration are directly utilized. The Chebyshev coefficients d_q^j and e_q^j are determined using the discrete orthogonality of the Chebyshev polynomials <cit.>. One  command is utilized to simultaneously compute all the components of the energy. The computational cost per processor scales as 𝒪(n_plN_d^p), which translates to an overall scaling of 𝒪(N) with respect to the number of atoms. 






 §.§ Atomic forces
 
In order to update the positions of the atoms during the course of the QMD simulation, the Hellmann-Feynman forces on the nuclei (Eq. <ref>) need to be computed.  


  
Local component
The local component of the force has the following discrete form in parallel computations

    𝐟_I^l = h^3 ∑_p=1^N_P∑_I' ∈ D_p,I^b∑_q ∈ K_p( ∇_h b_I'|_q ϕ_q + 1/2(b̃_q+b_q) ∇_h V_c,I'|_q + 1/2(∇_h b̃_I'|_q + ∇_h b_I'|_q ) V_c,q)

where the integrals in Eq. <ref> have been approximated using the trapezoidal rule (Eqn. <ref>) and D_p,I^b is the set of the I^th atom and its images whose pseudocharges have overlap with the processor domain Ω_p. 



  
Nonlocal component
The nonlocal component of the force in the SQ approach <cit.> takes the following form in parallel computations

    𝐟_I^nl = -4 ∑_p=1^N_P∑_I' ∈ D_p,I^c∑_q ∈ K_p_q^T  𝐕_nl,q^I'∇_h,q( ∑_j=0^n_pl c_q^j 𝐭_q^j)

where D_p,I^c is the set of the I^th atom and its images whose nonlocal projectors have overlap with the processor domain Ω_p, c_q^j is the coefficient of T_j in the polynomial expansion of the Fermi-Dirac function (Eq. <ref>); and 𝐭_q^j is determined using the recurrence relation:

    𝐭^i+1_q    =     2𝐇̂_q 𝐭^i_q-𝐭^i-1_q  ,    i=1, 2, …, n_pl
    𝐭_q^1     =    𝐇̂_q 𝐰_q  ,  𝐭_q^0 = 𝐰_q  ,

Note that unlike the iteration in Eq. <ref> which proceeds up to n_pl/2, the above iteration proceeds up to n_pl since some of the off-diagonal components of the density matrix are needed for the calculation of the nonlocal force in Eq. <ref>, i.e., ρ^j_q are not sufficient, rather 𝐭_q^j are required. 



  
Total atomic force
The atomic force in SQDFT is then calculated as

    𝐟_I     =     h^3 ∑_p=1^N_P∑_I' ∈ D_p,I^b∑_q ∈ K_p( ∇_h b_I'|_q ϕ_q + 1/2(b̃_q+b_q) ∇_h V_c,I'|_q + 1/2(∇_h b̃_I'|_q + ∇_h b_I'|_q ) V_c,q) 
        -     4 ∑_p=1^N_P∑_I' ∈ D_p,I^c∑_q ∈ K_p_q^T  𝐕_nl,q^I'∇_h,q( ∑_j=0^n_pl c_q^j 𝐭_q^j)  .


During the computation of the atomic forces, the values of c_q^j that were determined as part of the electron density calculation during the last SCF iteration are directly utilized. As mentioned previously, the nodal gradient matrix  ∇_h,q is not generated/stored explicitly, rather its product with 𝐭_q^j is calculated in matrix-free fashion. The total atomic force on all atoms is simultaneously computed using a single . Since the storage per processor scales as 𝒪(3 N_s) and the computational effort per processor scales as 𝒪(n_pl N_s N_d^p) + 𝒪(3 N_s)[This corresponds to the  calculation of the nonlocal component of the force, which is the dominant cost in the atomic force calculation in SQDFT. Note that 𝒪(3 N_s) arises due to the three matrix-vector products arising in the multiplication with the finite-difference gradient.], the overall  storage and computational cost scales as 𝒪(N) with respect to the number of atoms.




  
Relation to classical Fermi Operator Expansion (FOE)
The Clenshaw-Curtis SQ method bears some resemblance to the classical Fermi Operator Expansion (FOE) <cit.> in that both techniques use Chebyshev polynomials as the underlying basis for expanding the Fermi-Dirac function of a matrix. However, the key underlying difference is that the matrix in FOE corresponds to the Hamiltonian, whereas in SQ it corresponds to the nodal Hamiltonian. Therefore, truncation is automatically included within the SQ method and the key operation is reduced to local sparse matrix-vector products, as opposed to the global sparse matrix-matrix products in the FOE method. This 

makes SQ more efficient since (i) the Fermi level calculation does not require an outer loop, and (ii) inter-processor communication is needed just once per SCF iteration, unlike the FOE where it is required for every matrix-matrix multiplication. 

Finally, SQ also requires significantly less storage compared to FOE, making it especially well suited for modern high performance computing platforms.




§ RESULTS AND DISCUSSION
 
In this section, we demonstrate the accuracy, efficiency, and scaling of SQDFT in Kohn-Sham Density Functional Theory (DFT) calculations at high temperature. In all simulations, we employ a twelfth-order finite-difference discretization, norm-conserving Troullier-Martins pseudopotentials <cit.>, and the Local Density Approximation (LDA) <cit.> with the Perdew-Wang parametrization <cit.> of the correlation energy calculated by Ceperley-Alder <cit.>. To ensure accuracy of the standard 3s3p pseudopotentials employed, temperatures were limited to T ≲ 80000 K, where 2p states can be treated as fully occupied.[Calculations with deeper 2s2p3s3p pseudopotentials show 2p occupation of 5.9991 at 80000 K.
] 
To demonstrate applicability to metals and insulators alike, we consider two systems: (i) aluminum, a prototypical metal, and (ii) lithium hydride, a prototypical insulator. We compare the results obtained by SQDFT to benchmarks obtained by the finite-difference code SPARC <cit.> and planewave code ABINIT <cit.>, both of which solve the Kohn-Sham problem via diagonalization. 




 §.§ Accuracy and convergence
 
We first study the accuracy of SQDFT, i.e., we verify the convergence of computed energies and atomic forces with respect to key SQ parameters (i.e., quadrature order n_pl and truncation radius R_cut) as well as spatial discretization (i.e., mesh-size h). As representative systems, we consider a 32-atom cell of aluminum at the equilibrium lattice constant of 7.78 Bohr and a 64-atom cell of lithium hydride at the equilibrium lattice constant of 7.37 Bohr, with all atoms randomly displaced by up to 15 % of the equilibrium interatomic distance. We select the electronic temperature/smearing to be  σ=4 eV.

First, we verify the convergence of SQDFT energies and forces with respect to n_pl and R_cut in Fig. <ref>, with the reference diagonalization answers obtained by SPARC at the same mesh-size and a 4 × 4 × 4 Monkhorst-Pack grid for Brillouin zone integration.[Unlike standard codes, SQDFT can obtain the infinite crystal result without recourse to Brillouin zone integration.] We choose mesh-sizes of h=0.7780 and h=0.5264 Bohr for the aluminum and lithium hydride systems, respectively. It is clear that SQDFT obtains exponential convergence in the energy and atomic forces with respect to both parameters, in agreement with previous studies <cit.>. In particular, { n_pl,R_cut}∼{28,6} and { n_pl,R_cut}∼{40,6} are sufficient to obtain chemical accuracy in both the energy and forces for the aluminum and lithium hydride systems, respectively. Importantly, these values further reduce as the smearing/temperature is increased <cit.>, which makes SQDFT particularly attractive for high-temperature simulations. Note that neither the energy nor the atomic forces are variational with respect to n_pl and R_cut, hence the non-monotonic convergence in Fig. <ref>. 

 

Next, we verify the convergence of SQDFT energies and forces with mesh-size h to those computed by the established  planewave code ABINIT. To do so, we utilize n_pl=160 and R_cut = 10 Bohr in SQDFT, which are sufficient to put the associated errors well below the mesh errors of interest (see Fig. <ref>). In ABINIT, we employ a planewave cutoff of 50 Ha and a 4 × 4 × 4 Monkhorst-Pack grid for Brillouin zone integration, which results in energy and forces that are converged to within 10^-6 Ha/atom and 10^-6 Ha/Bohr, respectively. As shown in Fig. <ref>, both the energy and atomic forces in SQDFT converge rapidly and systematically, with chemical accuracy readily obtained. Notably, we see that energies and forces converge at comparable rates, without the need for additional measures such as double-grid <cit.> or high-order integration <cit.>. Therefore, accurate forces are easily obtained, as needed for structural relaxations and molecular dynamics simulations.

 




 §.§ Scaling and performance
 
We now study the scaling and performance of SQDFT on large-scale parallel computers with up to tens of thousands of processors. Specifically, we investigate the strong and weak scaling of SQDFT for aluminum and lithium hydride systems on the  and  supercomputers at the  Lawrence Livermore National Laboratory (LLNL) <cit.>. In all calculations, we employ a smearing of σ=4 eV and utilize (i) h=0.7780 Bohr and { n_pl,R_cut} = {28,6.224  Bohr} for the aluminum systems, and (ii) h=0.5264 Bohr and { n_pl,R_cut} = {40,6.387  Bohr} for the lithium hydride systems. These parameters are sufficient to obtain chemical accuracy of 0.001 Ha/atom and 0.001 Ha/Bohr in the energy and atomic forces, respectively, as demonstrated in the previous section.  

First, we perform a strong scaling study for a 2048-atom aluminum system on  and a 1728-atom lithium hydride system on , with atoms randomly displaced in both systems. For aluminum, the number of processors on  is varied from 64 to 8000. For lithium hydride, the number of processors on  is varied from 125 to 27000. The wall times per SCF iteration so obtained are presented in Fig. <ref>. Relative to the smallest number of processors, on the largest number of processors SQDFT achieves 97 % parallel efficiency for aluminum on  and 95 % for lithium hydride on . It is clear that SQDFT demonstrates excellent strong scaling. Notably, the wall time per SCF iteration for systems containing ∼ 2000 atoms can be reduced to less than 5 seconds.

Next, we perform a weak scaling study for aluminum and lithium hydride. For aluminum on , we increase the system size from 32 to 6912 atoms, while increasing the number of processors from 64 to 13824, maintaining two processors per atom for all systems. For lithium hydride on , we increase the system size from 8 to 10648 atoms, while increasing the number of processors from 27 to 35937, maintaining ∼ 3.4 processors per atom for all systems. The systems are generated by replicating 4-atom and 8-atom unit cells of aluminum and lithium hydride, respectively, with one atom in each unit cell randomly perturbed. We present the results so obtained in Fig. <ref>. We find the scaling with system size for aluminum on  to be 𝒪(N^1.00) and the scaling for lithium hydride on  to be 𝒪(N^1.01).

It is clear that SQDFT demonstrates excellent weak scaling, with near perfect 𝒪(N) scaling with respect to system size in practical calculations. 

Overall, the excellent strong and weak scaling of SQDFT up to tens of thousands of processors makes it possible to perform high-temperature Kohn-Sham molecular dynamics simulations at large length and time scales, as we show below.   

 
 



 §.§ High-temperature quantum molecular dynamics
 
We now consider the accuracy and efficiency of SQDFT in high-temperature Born-Oppenheimer quantum molecular dynamics (QMD) simulations. Such simulations are a cornerstone of modern warm dense matter theory <cit.>, providing equation-of-state and shock-compression predictions of unprecedented accuracy, up to temperatures of ∼ 100 eV and pressures of 100s of Mbar; see, e.g.,  <cit.>. At temperatures above ∼ 100 eV, conventional Kohn-Sham methods become prohibitively expensive, and so alternative methods such as OFMD <cit.> and, more recently, path integral Monte Carlo (PIMC) <cit.> have been employed to reach temperatures of 1000s of eV and higher. With sufficiently deep potentials, however, the new SQ methodology makes possible Kohn-Sham MD at temperatures of 1000s of eV as well.


As a representative example, we choose an 864-atom aluminum system and perform a 0.15 ps NVE QMD simulation with time step of 0.1 fs.[Due to the high velocities of the ions in high-temperature simulations, the time step has been chosen significantly smaller than in ambient calculations <cit.>.] We use a mesh-size h=0.7780 Bohr, quadrature order n_pl=28, truncation radius R_cut=6.224 Bohr, initial ionic temperature T = 116045 K, initial atomic positions close to perfect FCC crystal with one atom in each 4-atom unit cell randomly displaced by the same amount, and initial velocities randomly assigned based on the Maxwell-Boltzmann distribution. We integrate the equations of motion using the Leapfrog method <cit.>. The values of n_pl and R_cut have been chosen so as to put the associated errors close to an order of magnitude lower than the discretization error, which is ∼ 0.001 Ha/atom and ∼ 0.001 Ha/Bohr in the energy and forces, respectively. At each MD step, we set the electronic temperature equal to the ionic temperature, e.g, σ=10 eV at the start of the simulation. 

We perform the simulation on 3375 processors on  to obtain a wall clock time of ∼ 30 seconds per QMD step[Since the computational cost of SQDFT reduces with temperature, the wall time will further reduce as the temperature is increased.]. In Fig. <ref>, we plot the variation of  the total energy and temperature of the system over the course of the simulation. We observe that the temperature settles after ∼ 30 fs, subsequent to which the mean and standard deviation of the total energy are -3.3451 and 4.7 × 10^-4 Ha/atom, respectively. In addition, the drift in total energy as obtained from a linear fit is ∼ 2.7 × 10^-4 Ha/atom-ps. 


SQDFT thus shows excellent energy conservation, consistent with the accurate atomic forces obtained. 
We also plot the calculated radial distribution function 

in Fig. <ref>, which is in agreement with previous studies <cit.>. 


 


 




§ CONCLUDING REMARKS
 


We presented SQDFT: a large-scale parallel implementation of the Spectral Quadrature (SQ) method for 𝒪(N) Kohn-Sham density functional theory calculations at high temperature.
Specifically, we developed an efficient and scalable finite-difference implementation of the infinite-cell Clenshaw-Curtis SQ approach, in which results for the infinite crystal are obtained by expressing quantities of interest as bilinear forms or sums of bilinear forms, that are then approximated by spatially localized Clenshaw-Curtis quadrature rules. We demonstrated the accuracy of SQDFT by showing systematic convergence of energies and atomic forces with respect to quadrature order and truncation radius to reference diagonalization results, and convergence with mesh spacing to established planewave results, for both metallic and insulating systems. In all cases, chemical accuracy was readily obtained. We demonstrated excellent strong and weak parallel scaling on computer systems consisting of tens of thousands of processors, with near perfect 𝒪(N) scaling with system size, and wall clock times as low as a few seconds per SCF iteration for insulating and metallic systems of ∼ 2000 atoms. Finally, we verified the accuracy and efficiency of SQDFT in large-scale quantum molecular dynamics (QMD) simulations at high temperature, demonstrating excellent energy conservation and QMD step times of ∼ 30 seconds for an 864-atom aluminum system at ∼ 80000 K.

In the present work, we have focused on high-temperature Kohn-Sham DFT calculations. 

However, the SQ method is applicable at lower temperatures as well, with larger prefactor, as we show in Appendix <ref>. A possible approach to reduce this prefactor is to generate a localized orthonormal reduced basis (e.g., <cit.>), subsequent to which the SQ method is applied to the finite-difference Hamiltonian projected into this basis. This is indeed a promising path to 𝒪(N) DFT calculations of metals and insulators at ambient conditions which the authors are pursuing presently.






§ ACKNOWLEDGEMENTS

This work was supported in part by the National Science Foundation (Grant number 1333500), and performed in part under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. Early support from the Exascale Co-design Center for Materials in Extreme Environments supported by Office of Science Advanced Scientific Computing Research Program, and subsequent support from the Laboratory Directed Research and Development program at the Lawrence Livermore National Laboratory is gratefully acknowledged.








Appendix





§ AMBIENT TEMPERATURE KOHN-SHAM CALCULATIONS
 
Though the focus of the present work has been Kohn-Sham calculations at high temperature, SQDFT can also be utilized at ambient temperature, albeit with a larger prefactor. This increase in prefactor may, however, be mitigated by the excellent parallel scaling of SQDFT on large-scale parallel computers. 
In order to demonstrate this, we consider an 864-atom randomly perturbed aluminum system with smearing σ = 0.27 eV, as typical in calculations of metallic systems at ambient conditions in order to facilitate self-consistent convergence <cit.>. We utilize h=0.7780 Bohr and { n_pl,R_cut} = {320,18.672  Bohr}, which are sufficient to obtain chemical accuracy of 0.001 Ha/atom and 0.001 Ha/Bohr in the energy and atomic forces, respectively. We perform the calculations on , where the number of processors is varied from 1000 to 27000, the results of which are presented in Fig. <ref>. Relative to 1000 processors, the efficiency of SQDFT on 8000 processors is larger than 98 %, but on 27000 processors, the efficiency drops to 51 %. 

The reduced efficiency at this temperature at the largest processor counts arises due to the increased communications required for the larger nodal Hamiltonians (R_cut = 18.672  Bohr) relative to the computational work per processor. 
However, SQDFT is still able to achieve wall times of less than a minute per SCF iteration, which demonstrates its ability to perform large-scale Kohn-Sham quantum molecular dynamics simulations (QMD) even at ambient temperature, given sufficient number of processors.   



 



ReferenceStyle

84
natexlab#1#1
[#1],#1
[Hohenberg and Kohn(1964)]Hohenberg
authorP. Hohenberg, authorW. Kohn,
  journalPhysical Review volume136
  (year1964) pagesB864–B871.

[Kohn and Sham(1965)]Kohn1965
authorW. Kohn, authorL. J. Sham,
  journalPhysical Review volume140
  (year1965) pagesA1133–A1138.

[Goedecker(1999)]Goedecker
authorS. Goedecker, journalRev. Mod. Phys.
  volume71 (year1999) pages1085–1123.

[Bowler and Miyazaki(2012)]Bowler2012
authorD. R. Bowler, authorT. Miyazaki,
  journalReports on Progress in Physics volume75
  (year2012) pages036503.

[Aarons et al.(2016)Aarons, Sarwar, Thompsett, and
  Skylaris]aarons2016perspective
authorJ. Aarons, authorM. Sarwar,
  authorD. Thompsett, authorC.-K. Skylaris,
  journalThe Journal of Chemical Physics volume145
  (year2016) pages220901.

[Marx and Hutter(2009)]marx2009ab
authorD. Marx, authorJ. Hutter, titleAb
  initio Molecular Dynamics: Basic Theory and Advanced Methods,
  publisherCambridge University Press, year2009.

[Kresse and Hafner(1993)]kresse1993ab
authorG. Kresse, authorJ. Hafner,
  journalPhysical Review B volume47
  (year1993) pages558.

[Goedecker(1998)]goedecker1998decay
authorS. Goedecker, journalPhysical Review B
  volume58 (year1998) pages3501.

[Ismail-Beigi and Arias(1999)]ismail1999locality
authorS. Ismail-Beigi, authorT. Arias,
  journalPhysical review letters volume82
  (year1999) pages2127.

[Zhang and Drabold(2001)]zhang2001properties
authorX. Zhang, authorD. Drabold,
  journalPhysical Review B volume63
  (year2001) pages233109.

[Taraskin et al.(2002)Taraskin, Fry, Zhang, Drabold, and
  Elliott]taraskin2002spatial
authorS. Taraskin, authorP. Fry,
  authorX. Zhang, authorD. Drabold,
  authorS. Elliott, journalPhysical Review B
  volume66 (year2002) pages233101.

[Benzi et al.(2013)Benzi, Boito, and Razouk]benzi2013decay
authorM. Benzi, authorP. Boito,
  authorN. Razouk, journalSIAM Review
  volume55 (year2013) pages3–64.

[Soler et al.(2002)Soler, Artacho, Gale, Garcia, Junquera, Ordejon,
  and Sanchez-Portal]SIESTA
authorJ. M. Soler, authorE. Artacho,
  authorJ. D. Gale, authorA. Garcia,
  authorJ. Junquera, authorP. Ordejon,
  authorD. Sanchez-Portal, journalJ. Phys.: Condes.
  Matter volume14 (year2002)
  pages2745–2779.

[SIESTA: www.icmab.es/siesta(7 05)]SIESTAweb
SIESTA: www.icmab.es/siesta, yearaccessed 2017-07-05.

[Gillan et al.(2007)Gillan, Bowler, Torralba, and
  Miyazaki]Conquestref
authorM. J. Gillan, authorD. R. Bowler,
  authorA. S. Torralba, authorT. Miyazaki,
  journalComput. Phys. Commun. volume177
  (year2007) pages14–18.

[Conquest: www.order-n.org(7 05)]Conquestweb
Conquest: www.order-n.org, yearaccessed 2017-07-05.

[Skylaris et al.(2005)Skylaris, Haynes, Mostofi, and Payne]ONETEP
authorC. K. Skylaris, authorP. D. Haynes,
  authorA. A. Mostofi, authorM. C. Payne,
  journalJ. Chem. Phys. volume122
  (year2005).

[ONETEP: www.onetep.org(7 05)]ONETEPweb
ONETEP: www.onetep.org, yearaccessed 2017-07-05.

[Tsuchida(2007)]FEMTECKref
authorE. Tsuchida, journalJ. Phys. Soc. Jpn.
  volume76 (year2007).

[Osei-Kuffuor and Fattebert(2014)]MGMolref
authorD. Osei-Kuffuor, authorJ.-L. Fattebert,
  journalPhys. Rev. Lett. volume112
  (year2014).

[Mohr et al.(2014)Mohr, Ratcliff, Boulanger, Genovese, Caliste,
  Deutsch, and Goedecker]BigDFTref
authorS. Mohr, authorL. E. Ratcliff,
  authorP. Boulanger, authorL. Genovese,
  authorD. Caliste, authorT. Deutsch,
  authorS. Goedecker, journalJ. Chem. Phys.
  volume140 (year2014).

[BigDFT: bigdft.org(7 05)]BigDFTweb
BigDFT: bigdft.org, yearaccessed 2017-07-05.

[OpenMX: www.openmx-square.org(7 05)]OpenMXweb
OpenMX: www.openmx-square.org, yearaccessed 2017-07-05.

[Bock et al.(2014)Bock, Challacombe, Gan, Henkelman, Nemeth,
  Niklasson, Odell, Schwegler, Tymczak, and Weber]FreeONref
authorN. Bock, authorM. Challacombe,
  authorC. K. Gan, authorG. Henkelman,
  authorK. Nemeth, authorA. M. N. Niklasson,
  authorA. Odell, authorE. Schwegler,
  authorC. J. Tymczak, authorV. Weber,
  titleFreeON, year2014. noteLos
  Alamos National Laboratory (LA-CC 01-2; LA-CC-04-086), Copyright
  University of California.

[FreeON: www.openhub.net/p/freeon(7 05)]FreeONweb
FreeON: www.openhub.net/p/freeon, yearaccessed 2017-07-05.

[Ruiz-Serrano et al.(2012)Ruiz-Serrano, Hine, and
  Skylaris]RuiHinSky12
authorA. Ruiz-Serrano, authorN. D. M. Hine,
  authorC.-K. Skylaris, journalJ. Chem. Phys.
  volume136 (year2012).

[Graziani et al.(2014)Graziani, Desjarlais, Redmer, and
  Tricky]gradesred2014
editorF. Graziani, editorM. P. Desjarlais,
  editorR. Redmer, editorS. B. Tricky (Eds.),
  titleFrontiers and Challenges in Warm Dense Matter, Lecture Notes
  in Computational Science and Engineering, publisherSpringer,
  year2014.

[Graziani et al.(2012)Graziani, Batista, Benedict, Castor, Chen, Chen,
  Fichtl, Glosli, Grabowski, Graf, Hau-Riege, Hazi, Khairallah, Krauss,
  Langdon, London, Markmann, Murillo, Richards, Scott, Shepherd, Stanton,
  Streitz, Surh, Weisheit, and Whitley]grabasben2011
authorF. R. Graziani, authorV. S. Batista,
  authorL. X. Benedict, authorJ. I. Castor,
  authorH. Chen, authorS. N. Chen,
  authorC. A. Fichtl, authorJ. N. Glosli,
  authorP. E. Grabowski, authorA. T. Graf,
  authorS. P. Hau-Riege, authorA. U. Hazi,
  authorS. A. Khairallah, authorL. Krauss,
  authorA. B. Langdon, authorR. A. London,
  authorA. Markmann, authorM. S. Murillo,
  authorD. F. Richards, authorH. A. Scott,
  authorR. Shepherd, authorL. G. Stanton,
  authorF. H. Streitz, authorM. P. Surh,
  authorJ. C. Weisheit, authorH. D. Whitley,
  journalHigh Energy Density Physics volume8
  (year2012) pages105–131.

[Renaudin et al.(2003)Renaudin, Blancard, Clérouin, Faussurier,
  Noiret, and Recoules]renaudin2003aluminum
authorP. Renaudin, authorC. Blancard,
  authorJ. Clérouin, authorG. Faussurier,
  authorP. Noiret, authorV. Recoules,
  journalPhysical review letters volume91
  (year2003) pages075002.

[Dharma-Wardana(2006)]dharma2006static
authorM. Dharma-Wardana, journalPhysical Review E
  volume73 (year2006) pages036401.

[Ernstorfer et al.(2009)Ernstorfer, Harb, Hebeisen, Sciaini,
  Dartigalongue, and Miller]ernstorfer2009formation
authorR. Ernstorfer, authorM. Harb,
  authorC. T. Hebeisen, authorG. Sciaini,
  authorT. Dartigalongue, authorR. D. Miller,
  journalScience volume323 (year2009)
  pages1033–1037.

[White et al.(2013)White, Richardson, Crowley, Pattison, Harris, and
  Gregori]white2013orbital
authorT. White, authorS. Richardson,
  authorB. Crowley, authorL. Pattison,
  authorJ. Harris, authorG. Gregori,
  journalPhysical review letters volume111
  (year2013) pages175002.

[Lambert et al.(2006)Lambert, Clerouin, and Zerah]lamclegil2006
authorF. Lambert, authorJ. Clerouin,
  authorG. Zerah, journalPhys. Rev. E
  volume73 (year2006) pages016403.

[Zhang et al.(2016)Zhang, Wang, Kang, Zhang, and
  He]zhang2016extended
authorS. Zhang, authorH. Wang,
  authorW. Kang, authorP. Zhang,
  authorX. He, journalPhysics of Plasmas
  volume23 (year2016) pages042707.

[Cangi and Pribram-Jones(2015)]cangi2015efficient
authorA. Cangi, authorA. Pribram-Jones,
  journalPhysical Review B volume92
  (year2015) pages161113.

[Hu et al.(2016)Hu, Militzer, Collins, Driver, and
  Kress]humilkre2016
authorS. X. Hu, authorB. Militzer,
  authorL. A. Collins, authorK. P. Driver,
  authorJ. D. Kress, journalPhys. Rev. B
  volume94 (year2016).

[Suryanarayana(2013)]suryanarayana2013spectral
authorP. Suryanarayana, journalChemical Physics Letters
  volume584 (year2013) pages182–187.

[Pratapa et al.(2016)Pratapa, Suryanarayana, and
  Pask]pratapa2016spectral
authorP. P. Pratapa, authorP. Suryanarayana,
  authorJ. E. Pask, journalComputer Physics
  Communications volume200 (year2016)
  pages96–107.

[Prodan and Kohn(2005)]prodan2005nearsightedness
authorE. Prodan, authorW. Kohn,
  journalProceedings of the National Academy of Sciences of the
  United States of America volume102 (year2005)
  pages11635–11638.

[Anantharaman and Cancès(2009)]anantharaman2009existence
authorA. Anantharaman, authorE. Cancès,
  journalAnnales de l'Institut Henri Poincare (C) Non Linear
  Analysis volume26 (year2009)
  pages2425–2455.

[Pask and Sterne(2005)]Pask2005
authorJ. E. Pask, authorP. A. Sterne,
  journalPhys. Rev. B volume71 (year2005)
  pages113101.

[Gavini et al.(2007)Gavini, Knap, Bhattacharya, and
  Ortiz]Gavini2007
authorV. Gavini, authorJ. Knap,
  authorK. Bhattacharya, authorM. Ortiz,
  journalJournal of the Mechanics and Physics of Solids
  volume55 (year2007) pages669 – 696.

[Suryanarayana et al.(2010)Suryanarayana, Gavini, Blesgen,
  Bhattacharya, and Ortiz]Phanish2010
authorP. Suryanarayana, authorV. Gavini,
  authorT. Blesgen, authorK. Bhattacharya,
  authorM. Ortiz, journalJournal of the Mechanics and
  Physics of Solids volume58 (year2010)
  pages256 – 280.

[Suryanarayana et al.(2011)Suryanarayana, Bhattacharya, and
  Ortiz]Phanish2011
authorP. Suryanarayana, authorK. Bhattacharya,
  authorM. Ortiz, journalJournal of Computational
  Physics volume230 (year2011) pages5226 –
  5238.

[Suryanarayana and Phanish(2014)]Suryanarayana2014524
authorP. Suryanarayana, authorD. Phanish,
  journalJournal of Computational Physics volume275
  (year2014) pages524 – 538.

[Ghosh and Suryanarayana(2016)]ghosh2016higher
authorS. Ghosh, authorP. Suryanarayana,
  journalJournal of Computational Physics volume307
  (year2016) pages634–652.

[Bowler et al.(2002)Bowler, Miyazaki, and Gillan]bowler2002recent
authorD. Bowler, authorT. Miyazaki,
  authorM. Gillan, journalJournal of Physics: Condensed
  Matter volume14 (year2002) pages2781.

[Skylaris and Haynes(2007)]skylaris2007achieving
authorC.-K. Skylaris, authorP. D. Haynes,
  journalThe Journal of chemical physics volume127
  (year2007) pages164712.

[Suryanarayana(2017)]suryanarayana2017nearsightedness
authorP. Suryanarayana, journalChemical Physics Letters
  volume679 (year2017) pages146–151.

[Suryanarayana et al.(2013)Suryanarayana, Bhattacharya, and
  Ortiz]Phanish2012
authorP. Suryanarayana, authorK. Bhattacharya,
  authorM. Ortiz, journalJournal of the Mechanics and
  Physics of Solids volume61 (year2013)
  pages38 – 60.

[Lin and Yang(2013)]lin2013elliptic
authorL. Lin, authorC. Yang, journalSIAM
  Journal on Scientific Computing volume35 (year2013)
  pagesS277–S298.

[Banerjee et al.(2016)Banerjee, Suryanarayana, and
  Pask]banerjee2016periodic
authorA. S. Banerjee, authorP. Suryanarayana,
  authorJ. E. Pask, journalChemical Physics Letters
  volume647 (year2016) pages31–35.

[Anderson(1965)]anderson1965iterative
authorD. G. Anderson, journalJournal of the ACM (JACM)
  volume12 (year1965) pages547–560.

[Pulay(1980)]pulay1980convergence
authorP. Pulay, journalChemical Physics Letters
  volume73 (year1980) pages393–398.

[Pratapa et al.(2016)Pratapa, Suryanarayana, and
  Pask]pratapa2016anderson
authorP. P. Pratapa, authorP. Suryanarayana,
  authorJ. E. Pask, journalJournal of Computational
  Physics volume306 (year2016)
  pages43–54.

[Suryanarayana et al.(2016)Suryanarayana, Pratapa, and
  Pask]suryanarayana2016alternating
authorP. Suryanarayana, authorP. P. Pratapa,
  authorJ. E. Pask, journalarXiv preprint
  arXiv:1606.08740  (year2016).

[Shewchuk(1994)]Shewchuk1994
authorJ. R. Shewchuk, titleAn introduction to the
  conjugate gradient method without the agonizing pain, year1994.

[Alfè(1999)]alfe1999ab
authorD. Alfè, journalComputer Physics
  Communications volume118 (year1999)
  pages31–33.

[Gropp et al.(1999)Gropp, Lusk, and Skjellum]gropp1999using
authorW. Gropp, authorE. Lusk,
  authorA. Skjellum, titleUsing MPI: portable parallel
  programming with the message-passing interface, volume volume1,
  publisherMIT press, year1999.

[Ghosh and Suryanarayana(2017a)]ghosh2016sparc2
authorS. Ghosh, authorP. Suryanarayana,
  journalComputer Physics Communications volume216
  (year2017a) pages109 – 125.

[Ghosh and Suryanarayana(2017b)]ghosh2016sparc1
authorS. Ghosh, authorP. Suryanarayana,
  journalComputer Physics Communications volume212
  (year2017b) pages189 – 204.

[Gropp et al.(2014)Gropp, Hoefler, Thakur, and Lusk]gropp2014using
authorW. Gropp, authorT. Hoefler,
  authorR. Thakur, authorE. Lusk, titleUsing
  advanced MPI: Modern features of the message-passing interface,
  publisherMIT Press, year2014.

[Lanczos(1950)]lanczos1950iteration
authorC. Lanczos, titleAn iteration method for the
  solution of the eigenvalue problem of linear differential and integral
  operators, publisherUnited States Governm. Press Office Los
  Angeles, CA, year1950.

[Pratapa and Suryanarayana(2016)]pratapa2016numerically
authorP. P. Pratapa, authorP. Suryanarayana,
  journalMechanics Research Communications volume78
  (year2016) pages27–33.

[Gil et al.(2007)Gil, Segura, and Temme]gil2007numerical
authorA. Gil, authorJ. Segura, authorN. M.
  Temme, titleNumerical methods for special functions,
  publisherSIAM, year2007.

[Press(2007)]press2007numerical
authorW. H. Press, titleNumerical recipes 3rd edition:
  The art of scientific computing, publisherCambridge university
  press, year2007.

[Goedecker and Teter(1995)]goedecker1995tight
authorS. Goedecker, authorM. Teter,
  journalPhysical Review B volume51
  (year1995) pages9455.

[Goedecker and Colombo(1994)]goedecker1994efficient
authorS. Goedecker, authorL. Colombo,
  journalPhysical Review Letters volume73
  (year1994) pages122.

[Troullier and Martins(1991)]Troullier
authorN. Troullier, authorJ. L. Martins,
  journalPhysical Review B volume43
  (year1991) pages1993–2006.

[Perdew and Wang(1992)]perdew1992accurate
authorJ. P. Perdew, authorY. Wang,
  journalPhysical Review B volume45
  (year1992) pages13244.

[Ceperley and Alder(1980)]Ceperley1980
authorD. M. Ceperley, authorB. J. Alder,
  journalPhys. Rev. Lett. volume45
  (year1980) pages566–569.

[Gonze et al.(2002)Gonze, Beuken, Caracas, Detraux, Fuchs, Rignanese,
  Sindic, Verstraete, Zerah, Jollet, Torrent, Roy, Mikami, Ghosez, Raty, and
  Allan]ABINIT
authorX. Gonze, authorJ. M. Beuken,
  authorR. Caracas, authorF. Detraux,
  authorM. Fuchs, authorG. M. Rignanese,
  authorL. Sindic, authorM. Verstraete,
  authorG. Zerah, authorF. Jollet,
  authorM. Torrent, authorA. Roy,
  authorM. Mikami, authorP. Ghosez,
  authorJ. Y. Raty, authorD. C. Allan,
  journalComputational Materials Science volume25
  (year2002) pages478–492(15).

[Ono and Hirose(1999)]OnoHir99
authorT. Ono, authorK. Hirose, journalPhys.
  Rev. Lett. volume82 (year1999)
  pages5016–5019.

[Bobbitt et al.(2015)Bobbitt, Schofield, Lena, and
  Chelikowsky]BobSchChe15
authorN. S. Bobbitt, authorG. Schofield,
  authorC. Lena, authorJ. R. Chelikowsky,
  journalPhys. Chem. Chem. Phys.  (year2015).
  noteDOI: 10.1039/c5cp02561c.

[Lawrence Livermore National Laboratory (LLNL) high performance
  computing systems: https://computation.llnl.gov/computing/machine-catalog(6
  27)]LLNLwebMachines
Lawrence Livermore National Laboratory (LLNL) high performance computing
  systems: https://computation.llnl.gov/computing/machine-catalog,
  yearaccessed 2017-06-27.

[Zhang et al.(2016)Zhang, Driver, Soubiran, and
  Militzer]zhadrimil2016
authorS. Zhang, authorK. P. Driver,
  authorF. Soubiran, authorB. Militzer,
  journalHigh Energy Density Physics volume21
  (year2016) pages16–19.

[Zhang et al.(2017)Zhang, Driver, Soubiran, and
  Militzer]zhadrimil2017b
authorS. Zhang, authorK. P. Driver,
  authorF. Soubiran, authorB. Militzer,
  journalJ. Chem. Phys. volume146
  (year2017).

[Shamp et al.(2017)Shamp, Zurek, Ogitsu, Fratanduono, and
  Hamel]shazurham2017
authorA. Shamp, authorE. Zurek,
  authorT. Ogitsu, authorD. E. Fratanduono,
  authorS. Hamel, journalPhys. Rev. B
  volume95 (year2017).

[Driver et al.(2017)Driver, Soubiran, Zhang, and
  Militzer]drisoumil2017
authorK. P. Driver, authorF. Soubiran,
  authorS. Zhang, authorB. Militzer,
  journalHigh Energy Density Physics volume23
  (year2017) pages81–89.

[Zhang et al.(2017)Zhang, Driver, Soubiran, and
  Militzer]zhadrimil2017
authorS. Zhang, authorK. P. Driver,
  authorF. Soubiran, authorB. Militzer,
  journalPhys. Rev. E volume96
  (year2017).

[Rapaport(2004)]rapaport2004art
authorD. C. Rapaport, titleThe art of molecular dynamics
  simulation, publisherCambridge university press,
  year2004.

[Lin et al.(2012)Lin, Lu, Ying, and Weinan]lin2012adaptive
authorL. Lin, authorJ. Lu, authorL. Ying,
  authorE. Weinan, journalJournal of Computational
  Physics volume231 (year2012)
  pages2140–2154.

[Zhang et al.(2017)Zhang, Lin, Hu, Yang, and Pask]zhang2017adaptive
authorG. Zhang, authorL. Lin, authorW. Hu,
  authorC. Yang, authorJ. E. Pask,
  journalJournal of Computational Physics volume335
  (year2017) pages426–443.

[Kresse and Furthmüller(1996)]VASP
authorG. Kresse, authorJ. Furthmüller,
  journalPhysical Review B volume54
  (year1996) pages11169–11186.




