{'fraction_non_alphanumeric': 0.03184112828772909, 'fraction_numerical': 0.007758942426682905, 'mean_word_length': 4.581624821839711, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Large-scale ab initio simulations based on systematically improvable atomic basis

We present a first-principles computer code package (ABACUS) that is based on density functional theory and numerical atomic basis sets. Theoretical foundations and numerical techniques used in the code are described, with focus on the accuracy and transferability of the hierarchical atomic basis sets as generated using a scheme proposed by Chen, Guo and He [J. Phys.:Condens. Matter \textbf{22}, 445501 (2010)]. Benchmark results are presented for a variety of systems include molecules, solids, surfaces, and defects. All results show that the ABACUS package with its associated atomic basis sets is an efficient and reliable tool for simulating both small and large-scale materials.


I. INTRODUCTION
The density functional theory 1,2 (DFT) based first-principles methods are becoming increasingly important in the research fields of condensed matter physics, material sciences, chemistry, and biology. With the rapid development of supercomputers and the advances of numerical algorithms, nowadays it is possible to study the electronic, structural and dynamical properties of complicated physical systems containing thousands of atoms using DFT.
In these cases, the efficiency of widely used plane wave (PW) basis is largely limited, because of its extended nature. Instead, local bases, such as atomic orbitals, are the better choices.
Atomic orbitals have several advantages as basis sets for the ab initio electronic structure calculations in the Kohn-Sham scheme. 1,2 First, the basis size of atomic orbitals is much smaller compared to other basis sets, such as PW or real-space mesh. Second, the atomic orbitals are strictly localized and therefore can be combined with either the so-called linear scaling algorithms 3  While the analytical Gaussian-type orbitals have been well established for ab initio calculations in the quantum chemistry community for decades, the numerically tabulated atomic orbitals are getting more and more popular in the computational physics community. Several first-principles codes based on the numerical atomic orbitals have been developed in recent years, e.g., SIESTA, 6 OpenMX, 7 FHI-aims, 8 to name just a few, which aim at large-scale DFT calculations by exploiting the compactness and locality of numerical atomic orbitals.
However, the numerical atomic orbitals must be constructed very carefully to ensure both good accuracy and transferability. Furthermore, it would be highly desirable if the quality of the basis sets can be systematically improved in an unbiased way. Recently, some of us [Chen, Guo, and He (CGH)] proposed a new scheme 9,10 to construct systematically improvable optimized atomic basis sets for DFT calculations. Based on the CGH procedure for basis set generations, we have developed a DFT package from scratch, named Atomic-orbital Based Ab-initio Computation at UStc (ABACUS) here in the Key Laboratory of Quantum Infor-mation, University of Science and Technology of China (USTC). In the ABACUS package, besides the primary option of using numerical atomic orbitals as basis functions, PW basis can also be employed as an alternative choice. This dual basis feature is very useful for accuracy and consistency checks in benchmark calculations. For both basis set choices, the package uses normal conserving pseudopotential in the Unified Pseudopotential Format (UPF) that has been used in Quantum ESPRESSO. 11 The UPF pseudopotentials can be generated from the Opium package. 12 Regarding the exchange-correlation functionals, we have implemented the local (spin) density approximation [L(S)DA], and the generalized gradient approximation (GGA) as constructed by Perdew, Burke, and Ernzerhof (PBE).
In addition, semi-empirical van der Waals (vdW) corrected DFT scheme as proposed by Grimme (DFT-D2) 13 has also been implemented. Other advanced functionals such as hybrid functionals are currently under development and will be reported in a later work. At the level of LDA and GGA, the ABACUS package can do typical electronic structure calculations, structure relaxations, and molecular dynamics.
In this paper, we first describe the main features of the ABACUS package, as well as the major techniques that are used to implement DFT algorithms with atomic basis sets. In a previous study, 9 the CGH orbitals have been demonstrated to be accurate and transferable for the group IV and group III-V semiconductors. Here, we extend the tested systems to a larger range of elements, including the alkali elements, 3d transition metals, group VI and group VII elements, with focus on the structural and electronic properties of molecules, solids, surfaces, and defects. The results demonstrate that ABACUS with the CGH orbitals are highly reliable for both finite and extended systems. In particular, the basis set at the level of double-ζ plus polarization function (DZP) is an excellent choice to compromise between accuracy and computational cost, and can be safely used in production calculations in most situations.
The rest of paper is organized as follows. In Sec. II, we introduce the basic algorithms and numerical techniques. In Sec. III, we will demonstrate the performance of the ABA-CUS package, focus on the accuracy of the atomic orbitals generated using the CGH scheme, for a variety of benchmark systems. Finally, we summarize our work in Sec. IV.

II. METHODS
In this section, we first briefly recapitulate the basic formulation of solving Kohn-Sham equations in atomic basis (Sec. II A) to set up the stage. This is followed by a description of the main techniques used in ABACUS. Topics to be covered include the generation of the CGH atomic orbitals (Sec. II B), the construction of Hamiltonian matrix and overlap matrix (Sec. II C), the solvers for Kohn-Sham equations (Sec. II D), and finally the total energy and force calculations (Sec. II E).

A. The Kohn-Sham equation in atomic basis
The central task in DFT calculations is to solve the Kohn-Sham equation, 1,2 H KS Ψ n (r) = ǫ n Ψ n (r) , where ǫ n and Ψ n (r) are the Kohn-Sham eigenvalues and eigenfunctions for state n. Hartree atomic unit (e = = m e = 1) is used here and throughtout the paper. The Kohn-Sham HamiltonianĤ KS can be written as, whereT = − 1 2 ∇ 2 ,V ext (r),V H [ρ(r)], andV xc [ρ(r)] are the kinetic energy operator, the external potential, the Hartree potential, and the exchange-correlation potential, respectively. The Kohn-Sham HamiltonianĤ KS thus depends on the electron density ρ(r), which can be determined from the occupied Kohn-Sham orbitals Here for simplicity we assume that the system is spin-degenerate, and hence the spin index is omitted. Extending the algorithm described here to the spin-polarized case is straightforward and has been implemented in ABACUS.
Norm-conserving pseudopotentials are used to describe the ion-electron interactions. The external potentialV ext (r) in Eq. (2) contains the summation of the ion-electron potentials of all atoms plus, when they exist, applied external potentials. Therefore (in the absence of the applied external potential), wherev ps αi is a norm-conserving pseudopotential 14 for the i-th atom of element type α, and τ αi is the atomic coordinate in the cell R. The pseudopotential can split into a local part of the potentialv L α and separable fully non-local potentials 15vN L α , v ps α =v L α +v N L α .
The applied external potentials, e.g., electric fields, can be easily added to the local part of the potential, while the non-local pseudopotential can be written as, where |χ αlmn are non-local projectors, with l, m, n being the angular momentum, the magnetic momentum, and the multiplicity of projectors, respectively. In Eq. (6), l max and n max are the maximal angular momentum and the maximal multiplicity of projectors for each angular momentum channel, respectively.
The Kohn-Sham equation is usually solved within certain basis sets. The ABA-CUS package offers two choices of basis sets: the PW basis set and the atomic basis set.
The advantage to do so is that the results obtained using atomic basis sets can be directly compared to those obtained from PW basis sets for small systems, and thus provides valuable benchmarks for the former. This will be clearly seen in Sec. III where the benchmark results for a variety of systems are presented. However, since the PW algorithm has been well developed and documented, here we only focus on the algorithms of the atomic-basis implementation.
Without losing generality, we consider crystalline systems under periodic boundary conditions. The Kohn-Sham eigenfunctions in Eq. (1) then become Bloch orbitals which, within atom-centered basis set, can be expanded as, where φ µ (r − τ αi − R) are the atomic orbitals centering on the i-th atom of type α in the unit cell R. The orbital index µ is a compact one, µ = {α, i, l, m, ζ} with l being the angular momentum, m the magnetic quantum number, and ζ the number of atomic orbitals for a given l. Here n and k are the band index and Bloch wave vector, and c nµ,k are the Kohn-Sham eigen-coefficients. Finally N is the number of unit cells in the Born-von-Karmen supercell under the periodic boundary conditions. Using Eq. (7), the electron density within atom-centered basis sets can be computed as where f nk is the Fermi occupation factor, and N k is the number of k points in the Brillouin zone (BZ) sampling, which is typically equivalent to the number of real-space unit cells N in the Bloch summation. ρ µν (R) in Eq. 8 is the density matrix in real space, defined as Please note that in the last line of Eq. 8, we have assumed, without losing generality, the atomic orbitals to be real, i.e., φ * µ = φ µ . Given the expansion of the Kohn-Sham states in terms of atomic orbitals in Eq. (7), the Kohn-Sham equation Eq. (1) becomes a generalized eigenvalue problem, where H(k), S(k) and c k are the Hamiltonian matrix, overlap matrix and eigenvectors at a given k point, respectively. E k is a diagonal matrix whose entries are the Kohn-Sham eigenenergies. To obtain the Hamiltonian matrix H(k), we first calculate where µ, ν are atomic orbital indices within one unit cell, and φ µR =φ µ (r − R − τ αi ), φ ν0 =φ ν (r − τ βj ). The Hamiltonian matrix at a given k point can be obtained as, Similarly, the overlap matrix at a given k point is obtained as, where, The construction of H(k) and S(k), as well as solving Eq. (10) take most of the computational time. These two aspects will be discussed in more details in Secs II C and II D.
In many cases, when the investigated unit cell is large enough, a single Γ-point in the BZ is enough to get converged results. In these cases, both the Hamiltonian and overlap matrices are real symmetric matrices. In the ABACUS package, we treat the Γ-point only calculations separately to improve the efficiency.

B. Systematically improvable atomic basis sets
Before going into the construction processes of H µν (R) and S µν (R), here we introduce the CGH atomic orbitals that are used in ABACUS. The quality of atomic basis is essential to obtain accurate results. Unlike PW basis, with which the quality of the calculations can be systematically improved by simply increasing the PW energy cutoff, the way to generate high-quality atomic basis functions is much more complicated. In the last decades, considerable efforts have been devoted to developing high quality atomic orbitals. 7,16-19 ABA-CUS adopts a scheme proposed by CGH 9 to generate systematically improvable, optimized atomic basis sets.
An atomic basis function can be written as a radial function multiplied by spherical harmonics (in practice we use solid spherical harmonic functions, which are real functions), where the indices l, m, and ζ have the usual meanings of angular momentum quantum number, magnetic quantum number, and the multiplicity of the orbitals for l. One usually needs more than one radial functions for each angular momentum to improve the quality and transferability of the atomic basis sets.
In the CGH scheme, the radial function f lζ (r) is expanded in terms of a set of spherical Bessel functions (SBFs), with the coefficients of the SBFs yet to be determined, i.e., where j l (qr c ) is the SBF with radius cutoff r c . The possible q values are chosen such that j l (qr c )=0. A kinetic energy cutoff is chosen to determine the maximal value of q, and thus the number of SBFs. We set the SBFs to be strictly zero beyond the radius cutoff r c . In fact, they have been used directly as short-ranged basis set in first-principles calculations. 20,21 Because (almost) any function within the radius cutoff r c can be represented as a linear combination of SBFs, this gives us a large number of degrees of freedom for optimizing the atomic basis set.
To obtain optimized atomic basis, we vary the coefficients of the SBFs to minimize the spillage between the atomic basis set and a set of selected reference systems. The spillage S is a positive number defined as the difference between the Hilbert space spanned the atomic basis set and the wave functions calculated by PW basis, whereP is a projector spanned by the atomic orbitals, Here, S −1 (k) is the inverse of overlap matrix S(k) between numerical atomic orbitals. The spillage has been proposed before to measure the quality of a set of atomic basis. 17,22,23 We then use simulated annealing method to determine the coefficients that minimize the spillage. Sometimes the numerical orbitals obtained from this procedure have unphysical oscillations, which may lower the transferability of the basis set. To eliminate these unphysical oscillations, we further minimize the kinetic energy of each atomic orbital while keeping the spillage almost unchanged. More details of this procedure can be found in Ref. 9.
The reference systems used in the basis-generation procedure are very important to ensure the transferability of numerical atomic orbitals. The CGH scheme allows the users to choose freely the target systems to generate high-quality basis sets for different purposes. From our experience, homonulclear dimers are very good reference systems 16,24 for generating transferable atomic orbitals for general purposes. To avoid any possible bias of the basis set towards certain geometrical structure, an average over dimers of several different bond lengths (compressed or elongated) is taken as the target in the optimization procedure. We provide scripts to generate the atomic bases using diatomic molecules as reference systems.
The scripts set up PW calculations provided in the code for the dimers at various bond lengths. We remark that after the atomic basis sets are generated, for consistency, one must use the same pseudopotential and energy cutoff in later atomic orbitals based calculations as those used in the basis generation.
The CGH scheme is very flexible and easy to implement. One can choose freely the angular momentum of the orbitals, and the multiplicity of the radial functions for each angular momentum. All atomic orbitals are generated from the same procedure and criteria.
These orbitals form a sequence of hierarchial basis sets, which have a systematic convergence behavior towards the PW reference. Furthermore, without any assumptions of the shapes of radial functions f lζ (r), in principle we can get the fully optimized radial functions. As will be shown in Sec. III, the atomic orbitals generated in this way indeed show excellent accuracy and transferability for various systems.

C. Hamiltonian and overlap matrices construction
As mentioned above, in order to construct H(k) and S(k) at given k points using Eqs. (12) and (13), we need to first calculate H µν (R) and S µν (R). During the processes, we take the full advantage of the short-range nature of the atomic orbitals, i.e., only the matrix elements whose corresponding atomic orbitals have non-zero overlaps are evaluated. This is because each matrix element could be written as an integral in real-space grids, if this integral involves two spatially well-separated basis functions that are not overlapping, then the result of the integral should be zero. This feature leads to a sparse matrix and then O(N) scaling of the number of integrals, which is a significant advantage compared to PW based methods. Now we briefly discuss how each term in H µν (R) and S µν (R) is calculated. As is clear from Eq. (11), the Hamiltonian matrix has several components, which are computed by two different techniques, namely the two-center integral technique and the grid integral technique, respectively. First, the kinetic energy matrix T µν (R) = φ µR |T |φ ν0 , the non-local pseudopotential matrix V N L µν (R) = φ µR |V N L |φ ν0 , as well as the overlap matrix S µν (R) = φ µR |φ ν0 can be efficiently calculated by the two-center integral technique, 25 which has been described thoroughly in Ref. 6. The two-center integrals can be split into two parts: a one-dimensional integral over radial functions and an angular integral involving spherical harmonic functions. The radial integrals are tabulated for a wide range of distances between two orbitals once for all. The value for two orbitals within a reasonable distance can then be interpolated from the table. The angular integral involving spherical harmonic functions leads to quantities of so-called the Gaunt coefficients, which can also be easily calculated.
Therefore, the two center integrals can be evaluated very efficiently. 6 Further details on the two-center integral technique are given in Appendix A 1.
The matrix elements of the local potentials, V loc µν (R) = φ µR |V loc |φ ν0 , with are evaluated on a uniform grid in real space. Here V L (r) is the sum of all local pseudopotentials. We first evaluate the local potential V loc (r) on each grid point: the local pseudopotentials and Hartree potentials are calculated using techniques adapted from PW basis, i.e., they are first calculated in the reciprocal space, and then Fourier transformed to the real-space grid. The exchange-correlation potential can be directly evaluated on the real-space grid. Once we have V loc (r), the matrix elements V loc µν (R) are directly summed over the real-space grid. More details on the grid-based integral technique for the local potentials are given in Appendix A 2. The grid integrals are one of the most time consuming parts in the algorithms based on atomic orbitals. However, the computation efforts of the grid integrals only scale linearly with the system size, and can be easily parallelized. Recently, Lin et al. [27][28][29] developed the PEXSI technique, which provides an alternative way for solving the Kohn-Sham problem without using a diagonalization procedure. Com-pared to linear scaling approach, PEXSI does not rely on the nearsightedness principle either to truncate density matrix elements. In a Γ-point calculation, the basic idea of PEXSI can be illustrated as follows. Denote by M the number of atomic orbitals, Φ(r) = [φ 1 (r), · · · , φ M (r)] the collection of all atomic orbitals in the real space, andγ(r, r ′ ) the single particle density matrix in the real space. Then the PEXSI approach first expandsγ(r, r ′ ) using a P -term pole expansion asγ Here H, S, Γ are the Hamiltonian matrix, the overlap matrix and the single particle density matrix represented under the atomic orbital basis set Φ, respectively. ǫ F is the chemical potential or Fermi energy. The complex shifts {z l } and weights {ω ρ l } are determined through a simple semi-analytic formula, and takes negligible amount of time to compute. The number of terms of the pole expansion is proportional to log(β∆E), where β is the inverse of temperature and ∆E is the spectral radius, which can be approximated by the largest eigenvalue of the (H, S) matrix pencil. The logarithmic scaling makes the pole expansion a highly efficient approach to expand the Fermi operator.
At first it may seem that the entire Green's function-like object [H − (z l + ǫ F )S] −1 needs to be computed. However, if we target the electron density ρ(r) =γ(r, r), then only are actually needed. A selected inversion algorithm can be used to efficiently compute these selected elements of the Green's function, and therefore the entire electron density. The computational cost of the PEXSI technique scales at most as O(N 2 ). The actual computational cost depends on the dimensionality of the system: the cost for quasi-1D systems such as nanotubes is O(N) i.e. linear scaling; for quasi-twodimensional systems such as graphene and surfaces (slabs) the cost is O(N 1.5 ); for general three-dimensional bulk systems the cost is O(N 2 ). 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. This method is not only applicable to the efficient computation of electron density, but also to other physical quantities such as free energy, atomic forces, density of states and local density of states.
All these quantities can be obtained without computing any eigenvalues or eigenvectors. For instance, the atomic force for atom i in species α can be computed as where the first part is independent of PEXSI algorithm and will be discussed in the next subsection. The second term in Eq. 21 depends on the energy density matrix, which is written as This matrix is given again by pole expansions with the same poles as those used for computing the charge density, with different weights {ω E l }. For more detailed information we refer readers to Refs. 29 and 30. In order to use the PEXSI technique for multiple k-point calculations, we need to work with the Green's function of a non-Hermitian Hamiltonian that is only structurally symmetric. The massively parallel selected inversion method for non-Hermitian but structurally symmetric matrices are currently under development, and will be integrated into ABACUS in the future to perform large scale electron structure calculations with multiple k-point sampling.
Compared to existing techniques, the PEXSI method has some notable features: 1) The efficiency of the PEXSI technique does not depend on the existence of a finite Highest Occupied Molecular Orbital (HOMO)-Lowest Unoccupied Molecular Orbital (LUMO) gap, and can be accurately applied to general materials systems including small gapped systems and metallic systems. The method remains accurate at low temperatures. 2) 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 use of 10, 000 ∼ 100, 000 processors on high performance machines. 3) As a Fermi operator expansion based method, PEXSI allows the use of a hybrid scheme that combines density of states estimation based on Sylvester's law of inertia with Newton's method to obtain the chemical potential. This is a highly efficient and robust approach with respect to the initial guess of the chemical potential, and is independent of the presence of gap states. 4) 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, and this requirement is satisfied when atomic orbitals are used to discretize the Kohn-Sham Hamiltonian. The sequential version of PEXSI has been demonstrated before with ABACUS , 29 and the massively parallel version of PEXSI is recently integrated with SIESTA 30 for studying large scale systems with more than 10, 000 atoms with insulating and metallic characters on more than 10, 000 processors. The parallel PEXSI method is to be integrated with ABACUS .

E. Total energy and force calculations
Once the Kohn-Sham equation is solved, one can obtain the total energy of the system using the Harris functional, 31 where E band is the Kohn-Sham band energy, which is the summation over occupied Kohn-Sham orbital energies. E H , E xc and E II are the Hartree energy, the exchange-correlation energy, and the Coulomb energy between ions respectively. The second term in the above equation is the so-called double-counting energy arising from the Hartree and exchangecorrelation potential, which have been included in the band energy term. If the Kohn-Sham equation is solved by matrix diagonalization, then the band energy is the summation of the Kohn-Sham eigenvalues ǫ nk of all occupied bands, i.e., Alternatively, if the PEXSI method is chosen to be the Kohn-Sham equation solver, the band energy is calculated as (only for Γ-point now) where ρ is the density matrix and H is the Hamiltonian matrix. The Hartree and exchangecorrelation energies are calculated on a uniform real-space grid, and E II are calculated by the Ewald summation technique. 32 The forces acting on the ions are given by the derivative of the total energy with respect to the atomic coordinates. Analytical expressions for forces computed with atomic basis sets are more sophisticated than those with the PW basis sets. Besides the Feynman-Hellmann forces, Pulay forces 33 due to the change of the atomic basis sets during a structural relaxation should also be considered.
Therefore, we rewrite the total energy as the sum of two parts: E KS is the electronic part of the total energy in Eq. 23, while E II is the energy due to the Coulomb interactions between ions. Then the total force experienced by the i-ion of type α After some derivations (more details are given in Appendix A 3), we arrive at with Note that µ = {α, i, l, m, ζ} is the compact index for the atomic orbitals. The first term of the above equation is the Feynman-Hellmann force, 34 whereas the rest two terms yield the so-called Pulay forces. 33 Following Ref. 6, one can prove that the term related to ∂ρµν (R) where, is the element of "energy density matrix". In the above equation, ǫ nk is the band energy for band n at wave vector k. This term arises because the atomic orbitals are not orthogonal.
Similar to the total energy evaluation, different force terms are also evaluated using different techniques to maximize the efficiency of force calculations. Specifically, due to the long range tail of local pseudopotential in real space, it is better to calculate it in reciprocal space using PW basis set, and this is how Feynman-Hellmann force associated with the local pseudopotentials is implemented. Another advantage of this implementation is that the derivative of the local pseudopotential can be easily done in reciprocal space. By taking advantage of the short-range character of the non-local pseudopotential operators and atomic orbitals, the force terms including the Feynman-Hellmann force arising from the non-local pseudopotential operator, the Pulay forces arising from the kinetic energy operator, as well as the non-orthogonal forces, are calculated using two-center integral techniques. Here the number of orbitals on each angular momentum channel is (2l + 1), where l is the angular momentum quantum number. The polarization functions refer to orbitals that have higher l than the maximal one used in a SZ basis set. Specifically, in a SZ basis set that contains only one s orbital, the p orbitals are referred to the polarization functions; in a SZ basis set that has both s and p orbitals, the d orbitals are indicated as the polarization functions. Finally, for transition metals which use all s, p, d orbitals in a SZ basis set, the f orbitals are the polarization functions.
Two more parameters are needed to define a CGH atomic orbital. First, because CGH atomic orbitals are generated by an optimization procedure that is based on the results from PW calculations of target systems (typically diatomic molecules here), thus the generated set of atomic orbitals depend on a specific energy cutoff (E cut ) used in PW calculations.
Second, all CGH atomic orbitals are enforced to be strictly localized within a radius R cut , beyond which the atomic orbitals are set to be exactly zero. Table I lists both parameters for 24 elements that will be used in the followings to test physical properties of systems. Compared to a PW basis set which can be systematically increased to reach arbitrary accuracy in a calculation, the accuracy of the existing LCAO basis sets are known to be difficult to improve systematically. However, our construction strategy for numerical atomic orbitals described in Sec. II B can in principle guarantees a systematic convergence towards the PW accuracy for the target system. We note that, for general systems, the convergence behavior of an atomic basis set should be checked a posteriori. This is carried out separately for molecules in Sec. III B and for solids in Sec. III C. Among all these tests, first of all, let us look into one numerical issue that is very common in atomic-orbital based calculations -the so-called eggbox effect 37 -which occurs when the integrals of matrix elements are evaluated on a finite, uniformly spaced real-space grid.

A. Eggbox effect
The eggbox effect refers to the artificial rippling of the ground-state total energy as a function of the atomic displacements relative to an uniform real-space grid points. 37,38 Specifically, it arises from the numerical errors in evaluating the integrals of a Hamiltonian operator with respect to the local orbitals on a finite uniform real-space grid. This effect is completely artificial but considerably complicates the calculations of forces acting on atoms and phonon dispersions. Naturally, the denser the real-space grid is -corresponding to a higher energy cutoff in the reciprocal space, the smaller in magnitude that the rippling of the ground-state total energy will be. Conversely, if an atomic orbital in the real space is designed in a cautious way that the high-energy components of its Fourier transform are suppressed, have not been well understood, but these only occur at a smaller energy scale and have not caused any problems so far. As will be shown in Sec.III D, accurate structural relaxations can be carried out without further corrections for the eggbox effect.

B. Molecules
Having the eggbox effect under control, we now look into the convergence quality of atomic basis set for small molecules. Similar to the eggbox test cases, a cubic cell with a side length 20 Bohr is chosen here to avoid the artificial interactions between a molecule and its images.
Taking N 2 as an example, we plot its ground-state energy versus its bond distance in Fig. 2. The hierarchical basis sets using in this calculation range from SZ to QZTP. The reference is obtained from PW calculations. As illustrated in In Fig. 2, as the size of the atomic basis set increases, the total energy of the N 2 molecule converges systematically towards the PW limit.
From these curves, one can deduce the equilibrium bond length, atomization energy, and vibrational frequency of the N 2 molecule. These quantities obtained with atomic basis sets for a selected molecular set can be used to validate the convergence quality of the localized basis sets, when compared to the corresponding results obtained from PW calculations.
In the followings, we present the benchmark results for bond lengths, atomization energies, and vibrational frequencies for 11 chemically bonded diatomic molecules. Followed by the interaction energies of the S22 molecular test set, 40 obtained by the DFT-D method 13 as implemented in ABACUS. The first test is used to validate our methods for chemically bonded dimers while the second one is for weakly bonded systems.

Bond lengths
The bond length of a molecule is an important quantity. In Table II Table II, all the calculated equilibrium bond lengths systematically approach the corresponding PW results with increased atomic orbital basis sets. Specifically, the mean absolute error (MAE) for the DZP basis set is 0.018Å for both LDA and PBE calculations.

As shown in
This accuracy is sufficiently good for most practical purposes. When going beyond DZP to TZDP and QZTP basis sets, the MAEs become even smaller; at the QZTP level, the MAE is only 0.004Å for LDA and 0.003Å for PBE.
For alkali-metal elements, the construction of high-quality localized atomic orbitals 6,45 is highly challenging because these orbitals tend to be very diffusive and have a longer tail than the atomic orbitals of other elements. However, by using CGH orbitals, we found that a rather satisfactory description of the molecular bonding involving alkali metal atoms can be achieved, as can be seen from the examples of Na 2 and LiH in Table II, if a large cutoff radius, i.e., between 10 Bohr and 12 Bohr, is used. The same observation holds for alkali metal elements in bulk materials, as will be shown in Sec. III C.
The experimental values in Table II

Atomization energies
The atomization energy refers to the energy cost to split a molecule into individual atoms.
It is an important property in thermochemistry. The benchmark for the atomization energy are done using the same diatomic molecular set and hierarchal atomic orbitals as used for bond-length tests. The results are presented in Table III. For open-shell atoms, the spinpolarized, symmetry-broken solution usually has the lowest ground-state energy, and is thus taken here as the reference of calculated atomization energy. As can be seen from Table III, for all molecules, the atomization energies obtained with the atomic basis sets converge systematically towards the PW limit.
Quantitatively, however, the accuracy of atomization energies obtained with atomic basis sets is not so spectacular as the one of geometrical properties. In particular, the atomization energies from SZ basis are indeed too small; the corresponding MAE is over 1 eV. The MAE is improved by about a factor of 2 when the basis set goes from SZ to DZ, but is still far from satisfactory. A dramatic improvement is achieved at the DZP level where a MAE around 0.1 eV is an acceptable accuracy for most practical purposes. Nevertheless, to reach the "chemical accuracy" (1 kcal/mol ≈ 0.043 eV), one has to go further in the hierarchy of atomic basis sets. The "chemical accuracy" is almost reached with the TZDP basis (0.034 eV for LDA and 0.041 eV for PBE) and well reached with QZTP (0.022 eV for LDA and 0.026 eV for PBE). Here, the so-called basis-set error is much smaller than the errors of the energy functionals. For this diatomic molecular test set listed in Table III, on average PBE overbinds by 0.36 eV and LDA overbinds by 0.75 eV by comparing the converged PW results to the experimental ones.

Vibrational frequency
For a given molecule, the equilibrium bond length and the atomization energy reflect the position and depth of the minimum in its potential energy surface, respectively, and the vibrational frequency probes the curvature of the potential energy surface around the equilibrium geometry. The vibrational frequency can be measured directly by experiment and is a powerful probe of structural and bonding characteristics of a molecule. In

Weak interaction energy
In the previous tests, we have demonstrated the convergence behavior of our basis sets for chemically bonded diatomic molecules. In a sense, the good performance of these atomic basis sets for the tested molecules is not surprising, since the homo-nuclear diatomic molecules are the target systems when generating atomic orbitals. Now we turn to the study of bigger and more weakly bonded molecules -the S22 test set. 40 This test set contains 22 weakly interacting molecular complexes that include hydrogen bonding, dispersion interaction, and mixed bonding types. Since its inception, the S22 test set has been widely used to benchmark computational methods that deal with van der Waals (vdW) interactions.
Each member in the S22 test set is a molecular dimer that contains two monomers interacting with each other. The interaction energy is defined as the difference between the total energy of the dimer and the sum of the total energies of the two individual monomers in their fully relaxed geometries. In this work, the interaction energies of the S22 molecules are calculated using the DFT-D2 method of Grimme, 46 specifically PBE-D2, as recently implemented in ABACUS.
The calculations for S22 molecules are still done with the supercell approach with cubic boxes that are sufficiently large (up to 50 Bohr side length) to avoid artificial interactions between the molecule and it the periodic images. Fig. 3 presents the interaction energy differences of the S22 molecules calculated from two methods. The first is PBE-D2 with

Cohesive energies
The cohesive energies of the tested solids at their equilibrium lattice constants are presented in Table VI. Again, the general trend that the LCAO results systematically approach the PW limit can be observed. It is also interesting to point out that the LCAO basis sets show even better convergence behavior for cohesive energies of solids than atomization energies of molecules. For example, at the DZP level, the MAE is only 0.04 eV for the cohesive energies of solids, compared to 0.09 eV for the atomization energies of molecules. The same behavior holds for higher levels of basis sets such as TZDP and QZTP, where the MAEs of 0.01 eV for solids are also twice smaller than those of molecules.

Bulk Moduli
The bulk modulus is another key property of solids, reflecting the variation of the groundstate energy with respect to the unit cell volume around the equilibrium state. In Table VII  In the previous tests for molecules and solids, we have established the reliability of the ABACUS package and the accuracy of its associated LCAO basis sets. Here we test a more "challenging" problem -the reconstruction of the Si(100) surface. This surface is tech-    p(2 × 2) reconstructions, as shown in Fig. 4(a-c). It is well known that the Si atoms on the surface layer form dimers to lower the energy of the system by removing one of the two dangling bonds. 60 In the p(2 × 1)s reconstruction structure, the two atoms on the top layer come closer symmetrically, with bond length between them becomes slightly shorter than the nearest neighbor distance in the Si bulk. In the p(2 × 1)a case, the dimers buckle out of the Si(100) surface. Finally, in the p(2 × 2) case, the buckled dimers change their orientations alternatively (see Fig. 4(c)). Note that another c(2 × 4) reconstruction exists where the adjacent buckled dimer rows orientate oppositely. However, the energy lowering of this reconstruction is almost identical to that of p(2 × 2) when using DFT with LDA, 58 thus this fourth reconstruction structure is not considered in this work.
Next we describe the computational setup of our simulations. ABACUS is used with both LCAO DZP basis set and PW basis set. Also, the Quantum ESPRESSO (QE) package 11 is used, which serves as an independent check for the ABACUS package. The same pseudopotential for Si is used for all three simulations. We model the Si(100) surface by using a repeated slab which contains 12 atomic layers, because the structural distortion below the surface layer extends 4-5 layers into the bulk as noticed before. 61 The central two layers are fixed during the structural relaxation, and only atoms in the outermost five layers on each side are allowed to relax. The conjugate gradient algorithm is used for structural relaxation with a force threshold of 0.01 eV/Å. The slabs are separated by a vacuum of 30Å thick to avoid artificial interactions between neighboring slabs. When computing the energy differences of the p(2 × 1)s and p(2 × 1)a reconstructions with respect to the ideal surface, the (2×1) unit cell in the x-y plane is used, while in the case of the p(2 × 2) reconstruction, the (2×2) unit cell is used for both relaxed and ideal surfaces. In all calculations, a (6 × 6 × 1) k-point mesh is used.  Next we examine the relaxed structures obtained by the three methods. We choose the p(2 × 1) reconstruction because its structure distortion is most pronounced among the three.
Specifically, we look at the pentagon pattern formed by the atoms from the top three layers of the p(2 × 1) structure, as illustrated in Fig. 4(d). In Table IX,    that are often difficult and expensive to carry out. Since usually a large supercell is required to model defects in solids, the LCAO technique, which scales favorably with the system size, is a preferable choice for simulating the electronic structures of defects.
As a specific example, we employ the ABACUS package to study the group III-V Specifically, they are 1 × 1 × 2, 1 × 2 × 2, 2 × 2 × 2, 2 × 2 × 4, and 2 × 4 × 4 k-point meshes for system sizes ranging from 256 atoms to 16 atoms, respectively. The LDA is chosen to be exchange-correlation functional. The internal geometry of each supercell is fully relaxed. x increases. We note that, the geometries of the supercells in this work are chosen in a special way, i.e., by doubling the supercell size successively along the x, y, and z directions. Thus the exact behavior of the plot in Fig. 5 might slightly changes if the shape of the supercells is chosen differently. However, this should not alter the general trend we obtained. We also find that the closing of the band gap happens when x is larger than 0.125.

Formation energy of GaAs
The next relevant issue is the stability of the GaAs 1−x N x alloy. A key quantity here, for small x, is the formation energy of an N defect, which is defined as E f = E(Ga n As n−1 N) − E(Ga n As n ) + µ(As) − µ(N), where E(Ga n As n ) is the energy of a Ga n As n supercell containing n GaAs formula units, E(Ga n As n−1 N) is the energy of the above supercell but with one As atom replaced by an N impurity atom. µ(As) and µ(N) are the chemical potentials for the As atom and the N atom, respectively. The actual value of the atomic chemical potential depends on the chemical conditions in the experiment.
Here we consider two situations: the As-rich & N-rich condition and the As-poor & N-rich condition. In the first one, µ(As) = µ rich (As) is chosen to be the energy per atom in the yellow arsenic crystal form. In the second case, µ poor (As) is given by adding the formation enthalpy of bulk GaAs to µ rich (As). The formation enthalpy of bulk GaAs bulk is taken here as the atomization energy of GaAs per atom at zero temperature (2.00 eV). Finally, µ rich (N) is chosen to be one half of the total energy of the N 2 molecule. Obviously, the corresponding N concentration in the supercell Ga n As n−1 N is given by x = 1/n. We would like to point out that, by using Eq. (31), the N impurities are regarded as point defects, and the possible interactions between neighboring N defects are neglected.
In Fig. 6, the formation energy of an N defect is plotted as a function of the N concentration x for both As-rich (upper panel) and As-poor (lower panel) conditions. Calculations are done using the ABACUS package with two sets of LCAO bases: 1) DZP for all elements; 2) DZP for Ga and As, and TZDP for N. The results are compared to those obtained from the ABACUS/PW calculations, and those obtained by the Vienna Ab-initio Simulation Package (VASP) package with projected augmented wave method. 69,70 Fig. 6 shows that under the As-rich condition, the N defect formation energies of GaAs 1−x N x are positive for all N compositions, indicating the very low probability to dope N into GaAs under such condition. However, under the As-poor condition, the formation energies become negative for small N concentrations, which is consistent with the experimental finding that GaAs 1−x N x alloy can only be formed for narrow composition range near the endpoints. 63 The formation energy curves from ABACUS/DZP calculations follow the same trend as the ones from ABACUS/PW calculations for both As-rich and As-poor cases in Fig. 6.
However, there are noticeable differences of calculated formation energies between the two methods. Specifically, the formation energies associated with the DZP basis set are underestimated by 0.2 to 0.3 eV compared to the ones from PW calculations. This is due to the fact that µ(N) is calculated by half of the N 2 energy (cf. Eq. 31), and the N 2 total energy calculated with the DZP basis still has an appreciable difference from the PW reference result, as shown in Fig 2. In fact, by just increasing the N basis set from DZP to TZDP, the corresponding formation energy difference is largely improved to within 0.1 eV compared to PW results for all cases, as shown in Fig. 6. Finally, comparing the ABACUS/PW results and the VASP results reveals that there exists an appreciable difference between the norm-conserving psedudopotenital treatment and the projector augmented wave method.
Despite these differences, the general trend of the dependence of the formation energy on the N concentration is well reproduced within all calculations.
In conclusion, the study of GaAs 1−x N x alloy illustrates that the ABACUS package with its LCAO basis sets can be used for reliable defect calculations. We would also like to note that, for a thorough and faithful treatment of the phase stability problem of the GaAs 1−x N x alloy, one needs to consider the influence of finite temperatures and include the entropy effect, but this goes beyond the scope of the present work.

IV. SUMMARY
To summarize, in this paper we introduce a comprehensive first-principles package, named ABACUS , in which both plane waves and efficient localized atomic orbitals can be used for electronic-structure calculations. In particular, we present the mathematical foundation and numerical techniques behind the atomic-orbital-based implementation within this package.
The performance and reliability of the ABACUS package were benchmarked for a variety of systems containing molecules, solids, surfaces and defects. Furthermore, we show that the hierarchial atomic basis sets generated with the CGH scheme allows for a systematic convergence towards the plane-wave accuracy, and the DZP basis set offers an excellent compromise between accuracy and the computational load, and can be used in production calculations for most purposes. The package is currently under active development, with more features and functionalities are being implemented. With all these efforts, we expect the ABACUS package will become a powerful and reliable tool for simulations of large-scale materials.
Because S lµmµ,lν mν ,lm (R), G lµmµ,lν mν ,lm and T lµmµ,lν mν ,lm (R) are all independent of the coordinates of atoms, they can be tabulated at the beginning of the DFT calculations once for all. For any given distance between two atoms, the corresponding overlap matrix elements and kinetic energy matrix elements can be calculated efficiently by interpolation method.

Grid-based techniques
The Hamiltonian matrix elements V loc µν are evaluated on a uniform real space grid, with both atomic orbitals and local potentials presented on each grid point. The local potential is, where V L (r) = R αi v L α (r − τ αi − R) is the summation of all the local pseudopotentials for i-th atom of element type α. Plane wave basis and Fast Fourier Transform (FFT) techniques are used here to efficiently evaluate V L (r) and V H (r) on the grid. Because the local pseudopotential v L α (r) has a fairly long tail in real space, it is inefficient to calculate V L (r) directly on a real space grid. Therefore, V L (G) is first calculated in reciprocal space where S α (G) = i e −iG·τ αi is the structure factor. An FFT is carried out to bring V L (G) back to real space. From our tests, this construction processes of V L (r) only take a small portion of total computational time, even for systems containing thousands of atoms. This is different from the method used in Ref. 6, where the short-ranged neutral atom potentials are used. Using the same set of plane wave basis, the Hartree potential is also first evaluated in reciprocal space and then be brought back to real space by using an FFT. The full formula of Hartree potential is

Force calculations
As shown in the main text, the forces evaluated from the basis of atomic orbitals have four contributions, which are the Feynman-Hellmann force, the Pulay force, the force due to nonorthogonality of the atomic orbitals, and the Ewald force due to the Coulomb interactions between the ions.
The Ewald force can be calculated analytically using the Ewlad summation techniques, 32 and therefore not discussed here. We discuss here the techniques to evaluate the rest three terms as follows.
1. Feynman-Hellmann force, In Hamiltonian H, only the local and non-local pseudopotentials explicitly depend on the coordinates of ions. Thus we can further break F F H αi into two terms related to pseudopotentials. The first one is related to the local pseudopotential, and can be evaluated in reciprocal space, This method has been shown to be accurate and fast. 71 The second term involves the nonlocal pseudopotential, A non-local pseudopotential projector χ αilmn is also a one dimensional numerical orbital like atomic orbital. Therefore, the above equation can be evaluated efficiently using two-center integral technique introduced before.

Pulay force
The existence of Pulay force is due to the fact that the basis set is not complete, The kinetic energy term and the non-local pseudopotential term in F pulay αi can also be calculated by the two center integral technique, whereas the local potential term is evaluated by the grid integral technique. These steps are similar to the calculations of Hamiltonian matrix, except one needs to substitute φ µ and φ ν with dφµ dτ αi and dφν dτ αi . To evaluate the derivatives of an atomic orbital with respect to the coordinates of the atoms, one needs to deal with both radial atomic orbital and spherical harmonic function parts. The radial part can be calculated numerically on a one dimensional grid while the second part can be obtained analytically by using the real spherical harmonic functions.
3. The force arises from the fact that atomic orbitals are not orthogonal, where E µν (R) = 1 N k nk f nk ǫ nk c * µn,k c nν,k e −ik·R (A17) is the element of "energy density matrix". 6 The derivative of overlap matrix with respect to atomic coordinates is, where τ µ and τ ν are the atomic coordinates for orbital µ and ν, and D = R + τ ν − τ µ , is the distance between two orbitals. The further expansion of this term is In order to make Y lm (D) analytical at the origin, here it is multiplied by D l . d dD D −l S lµmµ,lν mν ,lm (D) can be calculated numerically using interpolation method while d dD Y lm (D)D l can be calculated analytically by using real spherical harmonic functions.