You are an expert at summarizing long articles. Proceed to summarize the following text: 
kohn - sham density functional theory ( ksdft ) is the most widely used framework for electronic - structure calculations , and plays an important role in the analysis of electronic , structural and optical properties of molecules , solids and nano - structures . the efficiency of ksdft depends largely on the computational cost associated with the evaluation of the electron charge density for a given potential within a self - consistent field ( scf ) iteration . the most straightforward way to perform such an evaluation is to partially diagonalize the kohn - sham hamiltonian by computing a set of eigenvectors corresponding to the algebraically smallest eigenvalues of the hamiltonian . the complexity of this approach is @xmath1 , where @xmath2 is the number of electrons in the atomistic system of interest . as the number of atoms or electrons in the system increases , the cost of diagonalization becomes prohibitively expensive . although linear scaling algorithms  @xcite are attractive alternatives for improving the efficiency of ksdft , they rely on using the nearsightedness principle  @xcite , which asserts that the density perturbation induced by a local change in the external potential decays away from where the perturbation is applied . one can then truncate elements of the density matrix away from the diagonal . such truncation can be in practice applied only to insulating systems whose density matrix elements decay exponentially away from the diagonal , but not to metallic systems at low temperature , for which the density matrix elements decay only algebraically away from the diagonal . the recently developed pole expansion and selected inversion ( pexsi ) method  @xcite provides an alternative way for solving the kohn - sham problem without using a diagonalization procedure , and without invoking the nearsightedness principle to truncate density matrix elements . compared to existing techniques , the pexsi method has a few salient features : 1 ) pexsi expresses physical quantities such as electron density , free energy , atomic forces , density of states and local density of states in terms of a spectral projector which is evaluated without computing any eigenvalues or eigenvectors . 2 ) the computational cost of the pexsi technique scales at most as @xmath3 ) . the actual computational cost depends on the dimensionality of the system : the cost for quasi-1d systems such as nanotubes is @xmath4 _ i.e. _ linear scaling ; for quasi-2d systems such as graphene and surfaces ( slabs ) is @xmath5 ; for general 3d bulk systems is @xmath6 . this favorable scaling hinges on the sparse character of the hamiltonian and overlap matrices , but not on any fundamental assumption about the localization properties of the single particle density matrix . 3 ) the pexsi technique can be accurately applied to general materials system including small gapped systems and metallic systems , and remains accurate at low temperatures . 4 ) the pexsi method has a two - level parallelism structure and is by design highly scalable . the recently developed massively parallel pexsi technique can make efficient usage of @xmath7 processors on high performance machines . 5 ) pexsi can be controlled with a few input parameters , and can act nearly as a black - box substitution of the diagonalization procedure commonly used in electronic structure calculations .    in order to benefit from the pexsi method , the hamiltonian and overlap matrices must be sparse . this requirement is satisfied if localized discretization is used for representing the kohn - sham hamiltonian by a finite sized matrix . examples of localized discretization include numerical atomic orbitals  @xcite , gaussian type orbitals  @xcite , the finite difference  @xcite and finite element  @xcite methods , adaptive curvilinear coordinates  @xcite , optimized nonorthogonal orbitals  @xcite and adaptive local basis functions  @xcite . in contrast , the plane - wave basis set is not localized and therefore can not directly benefit from the pexsi method . even though they are formally localized , the number of degrees of freedom per atom associated with methods such as the finite difference and the finite element is usually much larger than that associated with other methods such as numerical atomic orbitals , leading to an increase of the preconstant factor in the computational cost . therefore the finite difference and the finite element methods may not benefit as much from the pexsi technique as those based on numerical atomic orbitals . we note also that the use of hybrid functionals with an orbital - based exact exchange term  @xcite may significantly impact the sparsity pattern of the hamiltonian matrix and increase the computational cost .    in previous work  @xcite , the applicability of the pexsi method was demonstrated for accelerating atomic orbital based electronic structure calculations . with the sequential implementation of the pexsi method , it was possible to perform electronic structure calculations accurately for a nanotube containing 10,000 atoms discretized by a single-@xmath8 ( sz , minimal ) basis , and to perform geometry optimization for a nanotube that contains more than @xmath9 atoms with a double-@xmath8 plus polarization ( dzp ) basis . however , the sequential implementation does not benefit from the inherent parallelism in the pexsi method , and therefore leads to limited or no improvement for general electronic structure calculations . the contribution of this paper is twofold : 1 ) we present the siesta - pexsi method , which combines the siesta method  @xcite based on numerical atomic orbitals and the recently developed massively parallel pexsi method  @xcite . the siesta - pexsi method can be efficiently scalable to more than @xmath0 processors . we provide performance data for a range of systems , including strong and weak scaling characteristics , and illustrate the crossover points beyond which the new approach is more efficient than diagonalization . the accuracy of the result obtained from the siesta - pexsi method is nearly indistinguishable from the result obtained from the diagonalization method . 2 ) we develop a hybrid scheme of density of states estimation and newton s method to obtain the chemical potential . we demonstrate that the scheme is highly efficient and robust with respect to the initial guess , with or without the presence of gap states . the siesta - pexsi approach has been implemented as a new solver in siesta , with built - in heuristics that balance efficiency and accuracy , but at the same time offering full control by the user . this paper is organized as follows . in section [ sec : theory ] , we describe the massively parallel pexsi technique , and how to integrate it with the siesta method . we also present a new method to update the chemical potential . in section  [ sec : numerical ] , we report the performance of the siesta - pexsi method on several problems .    throughout the paper , we use @xmath10 to denote the imaginary part of a complex matrix @xmath11 . we use @xmath12 to denote the discretized hamiltonian matrix and the corresponding overlap matrix obtained from a basis set @xmath13 such as numerical atomic orbitals . similarly @xmath14 denotes the single particle density matrix operator , and the corresponding electron density is denoted by @xmath15 . the matrix @xmath16 denotes the single particle density matrix represented in the @xmath13 basis . in pexsi , @xmath16 and related matrices are approximated by a finite @xmath17-term pole expansion , denoted by @xmath18 respectively . however , to simplify notation , we will drop the subscript @xmath17 and simply use @xmath19 to denote the approximated matrices unless otherwise noted . the ground - state electron charge density @xmath15 of an atomistic system can be obtained from the self - consistent solution to the kohn - sham equations @xmath20 \psi_i(x ) = \psi_i(x ) \varepsilon_i , \label{kseqs}\ ] ] where @xmath21 is the kohn - sham hamiltonian that depends on @xmath15 , @xmath22 are the kohn - sham orbitals which in turn determine the charge density by @xmath23 with occupation numbers @xmath24 that can be chosen according to the fermi - dirac distribution function @xmath25 where @xmath26 is the chemical potential chosen to ensure that @xmath27 and @xmath28 is the inverse of the temperature , i.e. , @xmath29 with @xmath30 being the boltzmann constant .    the most straightforward method to solve the kohn - sham problem is to expand the orbitals @xmath31 as a linear combination of a finite number of basis functions @xmath32 , and thus recast as a ( generalized ) eigenvalue problem within an iterative procedure to achieve self - consistency in the charge density . the computational complexity of this approach is @xmath33 , where @xmath34 is the number of basis functions and is generally proportional to the number of electrons @xmath2 or atoms in the system to be studied . this approach becomes prohibitively expensive when the size of the system increases . formally , the electronic - structure problem can be recast in terms of the one - particle density matrix defined by @xmath35 with @xmath26 chosen so that @xmath36 = n_e$ ] . @xmath37 can thus be evaluated without the need for diagonalization , if the fermi function is approximated by a linear combination of a number of simpler functions . this is the fermi operator expansion ( foe ) method  @xcite .    while most of the foe schemes require as many as @xmath38 or @xmath39 terms of simple functions ( with @xmath40 the spectrum width ) , the recently developed pole expansion  @xcite is particularly promising since it requires only @xmath41 terms of simple rational functions . the pole expansion has the analytic expression @xmath42 we refer readers to ref .   for more details . the complex shifts @xmath43 and weights @xmath44 are determined only by @xmath45 and the number of poles @xmath17 . all quantities in the pole expansion are known explicitly and their calculation takes negligible amount of time . following the derivation in ref .  , we can use ( [ eqn : polerho ] ) to approximate the single particle density matrix @xmath37 by its @xmath17-term pole expansion , denoted by @xmath46 as    @xmath47    where @xmath48 $ ] is a collective vector of basis functions , @xmath49 , @xmath50 , and @xmath16 is an @xmath51 matrix represented in terms of the @xmath13 basis . to simplify our notation , we will drop the subscript @xmath17 originating from the @xmath17-term pole expansion approximation unless otherwise noted . it would seem that the need to carry out @xmath17 matrix inversions in ( [ eqn : gammapole ] ) would mean that the computational complexity of this approach is still close to the @xmath33 scaling of diagonalization . however , what is really needed in practice is just the electron density in real space , that is @xmath52 when the basis functions @xmath53 are compactly supported in real space , the product of two functions @xmath54 and @xmath55 will be zero when they do not overlap . these @xmath56 pairs can be excluded from the summation in eq .  . consequently , we _ only _ need @xmath57 such that @xmath58 in eq .  . this set of @xmath57 s is a subset of @xmath59 . to obtain these _ selected elements _ , we need to compute the corresponding elements of @xmath60 for all @xmath61 . we emphasize that we compute the selected elements of the density matrix because only these elements are needed to compute physical quantities such as charge density , energy and forces , due to the localized character of the basis set . the computed selected elements of the density matrix are accurate , and should be regarded as if we performed a conventional @xmath33 calculation first , and then only kept the corresponding selected elements of the density matrix . in principle we could retrieve any matrix element of the density matrix , simply by enlarging the set of `` selected elements '' . this process is fundamentally different from the usage of `` near - sightedness '' approximation , which throws away the information of the density matrix beyond the truncation region . the recently developed selected inversion method  @xcite provides an efficient way of computing the selected elements of an inverse matrix . for a ( complex ) symmetric matrix of the form @xmath62 , the selected inversion algorithm first constructs an @xmath63 factorization of @xmath11 , where @xmath64 is a block lower triangular matrix called the cholesky factor , and @xmath65 is a block diagonal matrix . the computational scaling of the selected inversion algorithm is only proportional to the number of nonzero elements in the cholesky factor @xmath64 , which is @xmath66 for quasi-1d systems , @xmath67 for quasi-2d systems , and @xmath68 for 3d bulk systems , thus achieving universal asymptotic improvement over the diagonalization method for systems of all dimensions . it should be noted that the selected inversion algorithm is an _ exact _ method for computing selected elements of @xmath69 if exact arithmetic is employed , and in practice the only source of error originates from the roundoff error . in particular , the selected inversion algorithm does not rely on any localization property of @xmath69 . in addition to computing the charge density at a reduced computational complexity in each scf iteration , we can also use this pole - expansion and selected - inversion ( pexsi ) technique to compute the free energy and the atomic forces efficiently without diagonalizing the kohn - sham hamiltonian . following the derivation in ref . , the relevant expressions are : @xmath70 + \mu n_e    - \frac12 \iint \frac{{\rho}(x )    { \rho}(y)}{{\left\lvertx - y\right\rvert } } { \,\mathrm{d}}x { \,\mathrm{d}}y\\    &    + e_{{{\mathrm{xc}}}}[{\rho } ]    - \int v_{{{\mathrm{xc}}}}[{\rho } ] { \rho}(x ) { \,\mathrm{d}}x ,    \end{split }    \label{eqn : helmholtzpole}\ ] ] @xmath71        + { \mathrm{tr}}\left [ \gamma^e \frac{\partial s}{\partial r_i } \right ] . \label{eqn : forcepole}\ ] ] where the energy and free - energy density matrices @xmath72 and @xmath73 are given by pole expansions with _ the same poles _ as those used for computing the charge density : @xmath74 since the trace terms in   and   require only the @xmath75th entries of @xmath76 for @xmath75 satisfying @xmath77 or @xmath78 , the needed elements of the energy - density matrices can be computed without additional complexity . in addition to its favorable asymptotic complexity , the pexsi method is also inherently more scalable than the standard approach based on matrix diagonalization when it is implemented on a parallel computer . the parallelism in pexsi exists at two levels . first , the selected inversions associated with different poles ( usually on the order of @xmath79 ) are completely independent . second , each selected inversion itself can be parallelized by using the parallel selected inversion method called pselinv  @xcite . a parallel selected inversion consists of the following steps :    1 .   the rows and columns of the matrices @xmath80 and @xmath81 are reordered to reduce the number of nonzeros in the triangular factor of the @xmath63 decomposition of @xmath82 . 2 .   a parallel symbolic factorization of @xmath82 is performed to identify the location of the nonzero matrix elements in @xmath64 . the numerical @xmath63 decomposition ( or equivalent @xmath83 decomposition ) of @xmath84 is performed . the desired selected elements of @xmath85 are computed from @xmath64 and @xmath65 . step 1 can be performed in parallel by using the parmetis  @xcite or the pt - scotch  @xcite software packages . its cost is much smaller compared to the numerical factorization , and only needs to be done once per scf cycle . although for symmetric matrices only @xmath86 factorization is needed , the pexsi package currently use the superlu_dist software package  @xcite to perform steps 2 and 3 in parallel with @xmath83 factorization . the @xmath83 factorization contains equivalent information as in @xmath86 factorization but can be twice as expensive in the worst case due to the lack of usage of the symmetry property of the matrix . pexsi and pselinv has a independent data structure , and allows to be interfaced with other sparse direct solvers . the cost of symbolic factorization in step 2 is usually much lower than the numerical factorization . the numerical factorization procedure can be described in terms of the traversal of a tree called the _ elimination tree_. each node of the tree corresponds to a block of continuous columns of @xmath87 . a node @xmath88 is the parent of a node @xmath89 if and only if @xmath90 in superlu_dist , each node is distributed among a subset of processors block cyclically . the traversal of the elimination tree proceeds from the leaves towards the root . the update of @xmath91 is performed in parallel on a subset of processors assigned to @xmath91 and its children nodes , and the main operations involved in the update are a number of dense matrix - matrix multiplications . in addition to parallelism within the update of a supernode , additional concurrency can be exploited in the traversal of different branches of the elimination tree for updating different supernodes . all of these techniques can be used in step 4 to compute selected elements of @xmath92 . in this step , the elimination tree is traversed from the root down towards the leaves . as each node is traversed , selected elements of @xmath92 within the columns that are mapped to that node are computed through a number of dense matrix - matrix multiplications . communication is needed among processors that are mapped to the node and its ancestors . multiple nodes belonging to different branches of the elimination tree can be traversed simultaneously if the update of these nodes do not involve communications with the same ancestor .    for sufficiently large problems , superlu_dist can achieve substantial speedup in the numerical factorization when hundreds to thousands of processors are used . similar or better speedup factors can be observed when selected elements of @xmath92 are computed from the distributed @xmath64 and @xmath65 factors . it is shown that for large matrices a single selected inversion can scale to @xmath93 cores , and if @xmath94 poles are used in the pole expansion of the fermi - dirac function , then total number of computational cores that can be efficiently utilized by the parallel pexsi method is @xmath95 . we will give some more concrete examples in section  [ sec : numerical ] to demonstrate the performance of our implementation of the parallel pexsi algorithm . compared to pexsi , it is generally difficult to efficiently use that many computational cores in a dense matrix calculation algorithm such as those implemented in scalapack . the main bottleneck of the computation using scalapack is the reduction of a dense matrix to a tridiagonal matrix and the back transformation of the eigenvectors of a tridiagonal matrix , which are inherently sequential . moreover , the cost of the diagonalization method scales cubically with respect to the matrix dimension . this limits the size of the matrix that can be handled by scalapack , as well as the parallel scalability on massively parallel computers . although some recent progress has been made to make these transformations more efficient  @xcite , making diagonalization scalable on more than tens of thousands of cores remains a challenging task . the chemical potential @xmath26 required in the pole expansion is not known a priori . its value must in general be determined so that the total number of electrons is appropriate : @xmath96 = { \mathrm{tr}}[\gamma \phi^t \phi ] = { \mathrm{tr}}[\gamma s ] . \label{eqn : nmu}\ ] ] as the right - hand - side is a non - decreasing function of @xmath26 , the chemical potential can be efficiently obtained by newton s method , maybe combined with extra safeguards such as bisection . the required derivative @xmath97 can be computed with very little extra cost using the pole expansion of the derivative of the fermi - dirac distribution @xmath98 , which can be constructed , for the reasons illustrated in ref .  , by using the same shifts @xmath61 as those in : @xmath99 ,      \label{}\ ] ] @xmath100 when newton s method is used , the convergence of @xmath26 is rapid near the correct chemical potential . however , the standard newton s method may not be robust enough when the initial guess is far away from the correct chemical potential . it may give , for example , too large a correction when @xmath101 is close to zero , as when @xmath26 is near the edge or in the middle of a band gap . one way to overcome the above difficulty is to use an approximation to the function @xmath102 to narrow down the region in which the correct @xmath26 must lie . this function can be seen effectively as a ( temperature smeared ) cumulative density of states , counting the number of eigenvalues in the interval @xmath103 . we can evaluate @xmath104 , its zero - temperature limit , without computing any eigenvalues of @xmath105 . instead , we perform a matrix decomposition of the shifted matrix @xmath106 , where @xmath64 is unit lower triangular and @xmath65 is diagonal . it follows from sylvester s law of inertia  @xcite , which states that the inertia ( the number of negative , zero and positive eigenvalues ) of a real symmetric matrix does not change under a congruent transform , that @xmath65 has the same inertia as that of @xmath107 . hence , we can obtain @xmath104 by simply counting the number of negative entries in @xmath65 . note that the matrix decomposition @xmath108 can be computed efficiently by using a sparse @xmath86 or @xmath83 factorization in real arithmetic . it requires fewer floating point operations than the complex arithmetic direct sparse factorization used in pexsi .    to estimate @xmath109 for a finite @xmath28 , we use the identity @xmath110 and perform an integration by parts to obtain @xmath111 the integral in can be evaluated numerically by sampling @xmath112 and @xmath104 at a number of quadrature points @xmath113 and performing a weighted sum of @xmath114 for @xmath115 . the @xmath116 evaluations of @xmath117 can be performed simultaneously using the @xmath63 factorization - based inertia counting procedure described above , with @xmath116 groups of processors . since the derivative of the fermi - dirac function is sharply peaked , we can approximate @xmath109 in a given interval by sampling @xmath104 in a slightly wider interval . fig . [ fig : ne ] shows the number of electrons at zero temperature @xmath104 obtained from inertia counting procedure , and the interpolated finite temperature profile @xmath118 at @xmath119 k for a dna system ( with finite gap ) and a sih system ( with zero gap ) near the fermi energy . while the finite temperature smearing effect is negligible for insulators at @xmath119 k , it is more pronounced for metals and leads to a more smooth @xmath118 which is suitable for applying newton s method to find the chemical potential . on the other hand , the inertia counting procedure obtains the global profile of the cumulative density of states , and does not suffer from the problem of being trapped in intermediate band gaps . this approximation to the function @xmath118 is then available for use in determining the approximate placement of the true @xmath26 by simple root finding : @xmath120 . in practice , we need to start from an interval which is large enough to contain the chemical potential . for systems with a gap the interval should also contain the gap edges . in this case , two roots are sought : @xmath121 and @xmath122 , where @xmath123 is small , say 0.1 , and so @xmath124 and @xmath125 will be estimates of the band edges from which @xmath26 can be determined as @xmath126 . the fidelity of @xmath118 , and thus the quality of the estimation of @xmath26 , depends on the density with which the interval can be sampled , which with a fixed number @xmath116 of sampling points for @xmath104 will increase as the interval is narrowed . in practice , the above procedure is repeated with progressively smaller intervals until either the @xmath26 estimate stabilizes or the size of the interval is small enough . unless the original interval is very large , 2 or 3 inertia counts are enough to provide an adequate @xmath26 from which to start the pexsi process , and the number of inertia counts can be one or even zero ( i.e. the inertia counting procedure is turned off ) in subsequent scf iterations . newton s method in the solver then takes over for the final refining of @xmath26 . the hybrid procedure of inertia counting and newton s method is efficient and robust for both insulating and metallic systems . this will be further demonstrated in section  [ sec : numerical ] . the density of states @xmath127 is defined as @xmath128 the factor of two comes from spin degeneracy . as advanced in the previous section , the cumulative density of states ( cdos ) @xmath129 is exactly the function @xmath104 discussed there . to evaluate @xmath127 , we sample @xmath104 at a set of @xmath130 using the inertia counting procedure , and use an appropriate interpolation scheme to approximate @xmath104 for other values of @xmath131 . the approximation to the dos is then obtained by numerical differentiation . we remark that inertia counting is not the only diagonalization - free method for computing the dos . other methods that make use of matrix vector multiplications only are also possible ( see _ e.g._a recent review  @xcite ) . another physical quantity that can be easily approximated via selected inversion is the localized density of states ( ldos ) , defined as @xmath132 and representing in fact the contribution to the charge density of the states with eigenvalues in the vicinity of @xmath131 , as filtered by the @xmath133 function . when a finite basis set @xmath13 is used , the ldos can be represented as @xmath134 where @xmath135 is a matrix of the coefficients of the expansion in @xmath13 of the orbitals @xmath136 , and @xmath137 is a diagonal matrix of the eigenvalues . it follows from the sokhotski - plemelj formula  @xcite @xmath138 that @xmath139 where symbol @xmath140 in stands for the cauchy principal value .    combining eq .   and , we obtain the following alternative expression for the ldos : @xmath141^{-1 }      \phi^{t}(x ) . \end{split }      \label{eqn : ldosgreen}\ ] ] note that eq .   allows us to compute the ldos without using any eigenvalue or eigenvector . in practice , we take @xmath142 to be a small positive number , and the local dos can be approximated by @xmath143^{-1}_{ij } ,      \label{eqn : ldosgreen}\ ] ] which is similar to the computation of the electron density . again only the selected elements of the matrix @xmath144^{-1}$ ] are needed for the computation of the ldos for each @xmath131 , and the selected elements can be obtained efficiently by the pselinv procedure . siesta is a density functional theory code which uses finite - range atomic orbitals to discretize the kohn - sham problem  @xcite , and thus handles internally sparse h , s , and single - particle density matrices . siesta is therefore well suited to implement the pexsi method along the lines explained above . the pexsi method can be directly integrated into siesta as a new kind of electronic - structure solver . conceptually , the interface between siesta and pexsi is straightforward . the existing siesta framework takes care of setting up the basis set and of constructing the sparse @xmath80 and @xmath81 matrices at each iteration of the self - consistent - field ( scf ) cycle . @xmath80 and @xmath81 are passed to the pexsi module , which returns the density matrix @xmath16 and , optionally , the energy - density matrix @xmath145 ( needed for the calculation of forces ) and the @xmath73 matrix that can be used to estimate the electronic entropy . siesta then computes the charge density to generate a new hamiltonian to continue the cycle , until convergence is achieved . energies and forces are computed as needed . the details of the interface are controlled by a number of parameters which provide flexibility to the user , especially in regard to the bracketing and tolerances involved in the determination of the chemical potential , and in the context of parallel computation . we list some of the more relevant parameters in table  [ tab : param ] , but we should note that even finer control is possible by other , more specialized parameters . .main siesta - pexsi parameters [ cols="^,^",options="header " , ]     another way to study the weak scaling of the siesta - pexsi approach is to examine how the total computational cost , which is the product of the wallclock time and the number of processors ( i.e. , the first two columns for each method in table  [ tab : exampleweakscaling ] ) , change with respect to the problem size . for perfect weak scaling , the change in cost should match that predicted by the computational complexity of the sequential algorithm . we plot the cost involved in solving each instance of the dna , c - bn and h@xmath146o problems in fig . [ fig : efficiencyweakscaling ] , which shows that the data for each series falls approximately on a line . the ( fitted ) slope of the cost for dna , which includes both factorization and selected inversion , is observed to be approximately 1.3 , higher than the @xmath147 ideal scaling for quasi-1d systems . for c - bn the slope of the fitting line is approximately 1.7 , which is a bit larger than the @xmath148 ideal scaling for quasi-2d systems . the slope of the line for the h@xmath146o series is approximately 2.18 , again slightly larger than the expected @xmath149 scaling for 3d systems . the observed degradation in parallel scaling ( relatively larger for dna ) comes about because neither selected inversion nor factorization ( which has the same asymptotic complexity ) can scale perfectly due to the complicated data communication and task dependency involved . it should be noted also that the ideal weak - scaling data in table [ tab : exampleweakscaling ] depends on the size and processor count chosen for the smallest instance of each series . we have already mentioned that all systems considered in our timings have been run with complete parallelization over poles . this limits the gains from increasing processor counts to the actual efficiency of selected - inversion and factorization within a pole , as a function of @xmath150 . roughly , the computational load for a given instance will depend on the number of non - zeros in the cholesky factor @xmath64 , which is plotted as a function of matrix size in fig .  [ fig : cholesky_nnz ] . it can be seen that c - bn has the highest such non - zero density , closely followed by h@xmath146o . this is because c - bn has a more densely packed structure than h@xmath151o , and thus a larger prefactor , but there is a clear trend for a faster increase in the number of non - zeros in the cholesky factor for h@xmath151o , in agreement with the asymptotic complexity of 3d and quasi-2d systems , respectively . the dna system has significantly fewer non - zero values in l for a given size of h , and so it can use a relatively smaller number of processors per pole efficiently .     for the systems shown on fig . [ fig : efficiencyweakscaling ] , as a function of the dimension of @xmath80 . ] we can also see from fig . [ fig : efficiencyweakscaling ] that the crossover point at which siesta - pexsi becomes more efficient ( in terms of cost ) than the standard diagonalization procedure in siesta is around matrix dimension @xmath34=7000 ( @xmath152700 atoms ) for the dna problem , @xmath34=50,000 ( @xmath1524,000 atoms ) for the c - bn problem and @xmath34=40,000 ( @xmath1525,000 atoms ) for the h@xmath146o problem . the actual orderings of the crossover points depend on the dimensionality ( through the power of n in the scaling ) , and also on a prefactor which is system dependent . since more processors can be used efficiently , the advantage of siesta - pexsi in time - to - solution performance is even more pronounced than the benefit from the reduction of total work load , and is already evident for smaller systems , as can be seen in table [ tab : exampleweakscaling ] . the conclusions regarding the performance comparison between pexsi and diagonalization , being based on the asymptotic scaling of the algorithm , remain largely valid even taking into account the efficiency improvements that can been obtained by refactoring some of the internals of scalapack  @xcite . number of processors ) for the first scf step for the h@xmath146o , c - bn , and dna examples . the points correspond to the various problem sizes as listed in table  [ tab : examplefeatures ] . the times for diagonalization show some variations due to differences in the construction of @xmath80 and the building of the @xmath16 matrix , and are represented by the gray stripe . the configurations and results can be found in more detail in table [ tab : exampleweakscaling ] . ]      as we discussed in section  [ subsec : chemicalpotential ] , the current implementation of siesta - pexsi uses a combination of an inertia counting bracketing technique and newton s method to determine the chemical potential @xmath26 that satisfies the nonlinear equation @xmath153 , where @xmath154 is defined in . in this section , we demonstrate the effectiveness of this approach for both metallic and insulating systems . the metallic example is the sih system with @xmath155 electrons . since our calculation is spin - restricted , the odd number of electrons is solely controlled through the fractional occupation of the h - derived state in the gap . any small perturbation of the chemical potential can change the number of electrons , and the perturbation is temperature dependent ( we use 300 k ) . despite the relatively small system size , this can be considered as the most difficult case for finding the chemical potential without having access to the eigenvalues . [ fig : convsih ] illustrates the number of inertia counting cycles and pexsi solver steps needed to obtain a chemical potential that satisfies @xmath156 for the converged sih system when using 60 poles . the calculation starts from a wide initial @xmath26 search interval of @xmath157 ev , and uses an adaptive electron tolerance which starts at a coarse @xmath158=0.1 and progressively tightens towards the target value of @xmath159 as the deviation from self - consistency in the h matrix elements moves towards its tolerance target of @xmath160 ry . the number of inertia counting steps needed to narrow down the interval that contains the true chemical potential is roughly @xmath161 in the first few scf steps . after the @xmath162th iteration , the inertia counting procedure is turned off . the number of pexsi solver steps needed is at most @xmath163 in the first few scf steps , and it becomes 2 after the @xmath164th scf step . [ fig : convsih ] also shows the absolute error @xmath165 at each scf step , which is always below the set tolerance . our test for insulating systems is the dna-1 system with @xmath166 electrons , at @xmath119 k and using 50 poles . [ fig : convdna1 ] illustrates the number of inertia counting and pexsi solver steps needed to reduce @xmath167 to a level below a final tolerance of @xmath159 , starting at a coarser tolerance of @xmath168 . the calculation also starts from a wide initial interval of @xmath157 ev . usually two inertia counting steps are needed in the first few scf iterations . the second iteration requires one more due to a large jump in h. after the @xmath162th iteration , the inertia counting procedure is completely turned off . throughout the scf cycle , which ends when the deviation from self - consistency in h is below @xmath160 ry , just one pexsi solver call is needed to provide a solution with a number of electrons within the tolerance . this is because @xmath26 is kept within the gap by the bracketing procedure .          in section [ subsec : dos ] , we described how to use siesta - pexsi to obtain spectral information of the atomistic system such as dos and ldos without computing eigenvalues and eigenvectors of @xmath105 . in fig . [ fig : dos ] ( a ) we plot the dos obtained from the standard diagonalization method and the inertia counting procedure implemented in siesta - pexsi for the c - bn@xmath169 system near the fermi level ( `` e@xmath170 '' ) . we observe nearly perfect agreement between the dos curves obtained from these two approaches . [ fig : dos ] ( b ) shows the dos near the fermi level for the dna-25 system using both the inertia counting method and diagonalization . the usage of pexsi for evaluating dos without obtaining eigenvalues is also significant for systems at large size . for the dna-25 system , using 64 @xmath150 and 3200 processors to evaluate @xmath171 points in the dos near the fermi energy ( each group of 64 processors evaluate @xmath172 points of cumulative dos with inertia counting ) only takes @xmath173 sec , while diagonalization using the same number of processors takes @xmath174 sec . it should be noted that the performance of pexsi can further be improved by simply using more cores : after the scf converges , one can restart the calculation to compute the dos using the converged density matrix , and use @xmath175 processors to compute @xmath171 points in the dos in parallel . in such case , the wall clock can be reduced to @xmath176 sec . [ fig : ldos ] presents the ldos for the sih system for an energy interval of width 0.4 ev around the fermi level , showing the state due to the added h atom in the bulk si system . we have combined the pole expansion and selected inversion ( pexsi ) technique with the siesta method for kohn - sham density functional theory ( kdsft ) calculation . the resulting siesta - pexsi method can efficiently use more than @xmath0 processors , and is particularly suitable for performing large scale _ ab initio _ materials simulation on high performance parallel computers . the siesta - pexsi method does not compute eigenvalues or eigenvectors of the kohn - sham hamiltonian , and its accuracy is fully comparable to that obtained from the standard matrix diagonalization based siesta calculation for general systems , including insulating , semi - metallic , and metallic systems at low temperature . the current implementation of siesta - pexsi does not yet support spin - polarized systems , but this can be achieved with minor changes to the code . furthermore , the code does not support * k*-point sampling , as in principle it is not needed for large - enough systems . however , in some cases , as in graphene and similar semi - metallic systems , an appropriate computation of the spectral properties such as the dos does need * k*-point sampling even when very large supercells are used . we will address this problem in future work , which requires modifying pexsi to handle complex hermitian @xmath80 matrices . in such case the shifted matrix @xmath177 is no longer symmetric but is only structurally symmetric . this further development will also be useful for _ ab initio _ study of non - collinear magnetism and spin - orbit coupling effects . this work was partially supported by the laboratory directed research and development program of lawrence berkeley national laboratory under the u.s . department of energy contract number de - ac02 - 05ch11231 , by scientific discovery through advanced computing ( scidac ) program funded by u.s . department of energy , office of science , advanced scientific computing research and basic energy sciences , by the center for applied mathematics for energy research applications ( camera ) , which is a partnership between basic energy sciences and advanced scientific computing research at the u.s department of energy ( l. l. and c. y. ) , by the european community s seventh framework programme [ fp7/2007 - 2013 ] under the prace project grant agreement number 283493 ( g. h. ) , and by the spanish mineco through grants fis2009 - 12721-c04 - 03 , fis2012 - 37549-c05 - 05 and csd2007 - 00050 ( a. g. ) . 38ifxundefined [ 1 ] ifx#1 ifnum [ 1 ] # 1firstoftwo secondoftwo ifx [ 1 ] # 1firstoftwo secondoftwo `` `` # 1''''@noop [ 0]secondoftwosanitize@url [ 0 ] + 12$12  & 12#1212_12%12@startlink[1]@endlink[0]@bib@innerbibempty @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * , ( ) @noop * * , ( ) @noop * * ,   ( ) @noop ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * , ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) @noop * * ,   ( ) link:\doibase 10.1016/j.parco.2011.05.002 [ * * ,   ( ) ] @noop * * ,   ( ) @noop ( ) @noop _ _  ( , ,  )
we describe a scheme for efficient large - scale electronic - structure calculations based on the combination of the pole expansion and selected inversion ( pexsi ) technique with the siesta method , which uses numerical atomic orbitals within the kohn - sham density functional theory ( ksdft ) framework . 
 the pexsi technique can efficiently utilize the sparsity pattern of the hamiltonian and overlap matrices generated in siesta , and for large systems has a much lower computational complexity than that associated with the matrix diagonalization procedure . 
 the pexsi technique can be used to evaluate the electron density , free energy , atomic forces , density of states and local density of states without computing any eigenvalue or eigenvector of the kohn - sham hamiltonian . 
 it can achieve accuracy fully comparable to that obtained from a matrix diagonalization procedure for general systems , including metallic systems at low temperature . 
 the pexsi method is also highly scalable . with the recently developed massively parallel pexsi technique , we can make efficient use of more than @xmath0 processors on high performance machines . 
 we demonstrate the performance and accuracy of the siesta - pexsi method using several examples of large scale electronic structure calculations , including 1d , 2d and bulk problems with insulating , semi - metallic , and metallic character .