Highly accurate local basis sets for large-scale DFT calculations in conquest

Over the last thirty years, density functional theory (DFT) has emerged as the leading electronic structure modelling technique, used in fields as diverse as biochemistry and electronic engineering, alongside physics, chemistry and materials science. The ability to model the atomic and electronic structure of molecules, liquids, nanoparticles and solids has become a key part of the scientific method.

However, almost all simulations are performed on a very restricted number of atoms: typically a few hundred, and very rarely beyond a thousand. There are various reasons for the restriction, the most often-cited of which is scaling: the asymptotic scaling with the number of atoms of standard DFT is cubic in time and quadratic in memory. While there has been considerable effort in developing linear scaling approaches to DFT^1,2) which have been demonstrated to be capable of calculations on millions of atoms,^3,4) these introduce approximations (typically integrating over energy, so that eigenstates are not easily available, though can be calculated within a range.⁵⁾) There is therefore a significant gap in the sizes of systems that can be modelled exactly, lying in the range of thousands to tens of thousands of atoms.

There are numerous examples of simulations that require this size of system, for instance: realistic doping in semiconductors; biological molecules with or without explicit water; nanostructures; compounds with very dilute compositions; and large-scale defects such as dislocations. There are, of course, many possible approaches to large-scale electronic structure calculations. Approximate methods, such as density functional tight binding^6,7) are effective, but introduce some empiricism and/or fitting, reduce basis sets and remove explicit calculation of certain terms. Full DFT calculations can be made both more efficient and to scale better in parallel by judicious choice of basis sets. The largest example of a full DFT calculation used a real-space finite-difference approach to perform calculations on 10 000+ atom systems⁸⁾ with a demonstration that 100 000 atoms are possible.⁹⁾

Localised basis sets are a common choice for efficient, large-scale DFT approaches. Using such a basis, it is easy to form sparse matrices, and to diagonalise the Hamiltonian exactly. Examples of these approaches include the pseudo-atomic orbitals (PAOs) found in the widely-used siesta code¹⁰⁾ and the Gaussian basis sets in the CP2K code.¹¹⁾

Our focus in this paper is the large-scale DFT code, conquest, which is capable of exact calculations on systems with several thousand atoms^5,12) and linear scaling calculations on millions of atoms.^3,4) conquest uses a basis set of support functions (SFs) to represent both the Hamiltonian and either the density matrix or the Kohn–Sham eigenstates. These SFs can be further represented by other functions: either a systematic basis of blip functions;¹³⁾ or numerical PAOs¹⁴⁾ (though conquest often uses a 1:1 mapping of PAOs to SFs). conquest can solve for the electronic ground state using three key methods: linear scaling DFT^3,15,16) which we will not discuss further here; multi-site SFs,^17–20) where a small number of SFs are made from the PAOs from several atoms (hence multi-site), with the Hamiltonian solved by exact diagonalisation; and the primitive PAO basis set (one PAO is one SF), also solved using diagonalisation.

PAOs are convenient and easily generated, but are not systematic in their convergence. In this paper, we will describe our approach to generating basis sets of varying sizes, and will demonstrate that it is possible to produce relatively modest basis sets that closely reproduce converged plane wave results using the same pseudopotentials.

The conquest code has been described in detail elsewhere;^{15,16,20–22)} nevertheless, some details of the implementation may be useful to assist the reader, and we summarise these here.

The central part of the conquest code relates to the creation of matrices and their use in finding the ground state of the system being studied. When using PAOs as the basis set (as in this paper), we form the overlap matrix elements, and the kinetic and non-local pseudopotential contributions to the Hamiltonian matrix elements using the high accuracy approach pioneered by the siesta code.¹⁰⁾ The remaining parts of the Hamiltonian matrix elements (from the local part of the pseudopotential as well as the Hartree and exchange-correlation potentials) are calculated by integration on a uniform spatial grid.

The ground state of the system is found using either exact diagonalisation of the general eigenproblem (using Scalapack routines) or a linear scaling approach.¹⁵⁾ For exact diagonalisation, Bloch's theorem is applied and a Monkhorst–Pack mesh of points in reciprocal space is generated. For linear scaling, which is purely real-space, this is not needed. Forces²²⁾ and stresses are calculated as exact derivatives of the energy (including Pulay forces and stresses from the movement of the basis functions).

2.1. Pseudopotentials

All basis sets have advantages and disadvantages; the well-known advantages of plane waves (simplicity and systematic convergence) apply for relatively small systems, where the efficiency of fast Fourier transforms (FFTs) also applies; for large systems, however, the poor parallel scaling of FFTs becomes significant. Moreover, plane waves, being solutions for the free electron, are poorly adapted to atomic calculations (leading to the use of pseudopotentials) and more generally to the solid state (leading to large numbers of basis functions).

There are approaches to systematic basis sets with localised orbitals, where SFs (also known as non-orthogonal generalised Wannier functions, among other names) are defined by a radius and a grid spacing equivalent to a plane wave cutoff. Examples of these include: the blip functions used in conquest;¹³⁾ periodic sinc functions;^31,32) wavelets;³³⁾ and Lagrange functions.³⁴⁾ It is also possible to use the grid points directly as basis functions with real-space techniques³⁵⁾ including finite elements^36,37) and finite differences.³⁸⁾

There is, however, a simplicity and intuition which comes from using atomic orbitals, or (in the case of a pseudopotential calculation) PAOs, which consist of radial functions, normally tabulated on a fine radial mesh, multiplied by appropriate spherical harmonics. Integrals between basis functions can be simplified, with angular parts found analytically and radial grids so fine as to give effectively analytic results. The determination of the basis set then relies on three questions: how are the radial functions calculated?; how large an angular momentum is required?; and how many functions should there be for each angular momentum?

The question of angular momentum is in large part determined by the valence electrons being considered: these angular momenta must be represented in the basis set. It is then useful to consider polarisation functions; following the Siesta code¹⁰⁾ we can consider this as the effect of a local electric field on the highest valence shell. For elements up to the lanthanides, this ensures that the PAO basis set includes angular momenta up to l = 2; unless an atom contains valence f electrons, this is generally sufficient.

The radial functions are always found with some form of confinement (to ensure sparsity of matrices and efficiency). While this confinement can be a simple radial cutoff, it is not clear how to choose that cutoff: a smaller cutoff will give greater efficiency, but may not yield as good a basis function, particularly in materials where long-range or weak interactions are key. The Siesta code introduced the idea of a uniform energy shift for all PAOs¹⁰⁾ to allow for a consistent definition of confinement for different angular momenta. It also introduced a flexible confinement potential to allow optimisation of the radial functions.³⁹⁾

Of course, with more than one radial function per angular momentum it is possible to choose different radii, but the question of how to set these is still difficult. The siesta code uses a split norm idea for subsequent radial functions: using a smooth polynomial from the origin to a certain point, and beyond that point matching the original. OpenMX calculates a set of five or six excited states (solutions with increasing numbers of nodes) for each angular momentum in a hard confining potential⁴⁰⁾ and then combines these primitive functions into a number of composite functions (typically with two radial functions per valence state and one per polarisation state), fitted to dimers or diatomic molecules.^41,42) It is also possible to optimise confining potentials by fitting to the eigenstates of converged plane wave calculations.⁴³⁾ As well as PAOs, Gaussian orbitals are in use, for instance in CP2K;⁴⁴⁾ Gaussian orbitals are the primitive basis for the original method on which conquest's multi-site SFs are based.¹⁷⁾

In developing PAO basis sets for conquest, we have considered the fact that the basis functions need to give variational flexibility to the Kohn–Sham eigenstates, and therefore need to span a range of radii. The simplest way, then, is to start with the idea of energy shifts and to specify radii on the basis of one energy that will give a highly confined orbital and one that will give a loosely confined orbital. We have found, empirically, that energy shifts of 2 and 0.02 eV fulfil these criteria well, and have also found that for situations where we specify a third radial function, using an intermediate value of 0.2 eV is effective. In fact, for most systems this gives three relatively evenly spaced confinement radii.

We define one set of PAOs by this equal energy criterion: a basis set with one radial function uses 0.02 eV; a set with two functions uses 2 and 0.02 eV; and with three functions uses 2, 0.2 and 0.02 eV. We then define four basis sizes: minimal (one radial function, valence only or SZ); small (one radial function for valence and polarisation, or SZP); medium (two radial functions for valence and one for polarisation, or DZP); and large (three radial functions for all angular momenta, or TZTP) though in this paper we will not test the minimal basis.

We also consider an alternative, specifying equal radii for all angular momenta. We do this by finding the radii for all angular momenta for a given confinement energy, and taking their mean. This gives a uniform set of radii for all angular momenta. This option tends to be more efficient, as the overall radius of the largest SF is smaller than the equal energy construction; however, it can require a finer integration grid as the more diffuse functions are compressed. We will see below that the two approaches typically produce similar results, though the equal radii approach is maybe a little better overall.

Our concern in this section is to test the accuracy of different PAO basis sets against the converged plane wave result. For this purpose we use the PWSCF code from the QuantumEspresso suite,⁴⁵⁾ which reads the same pseudopotentials as conquest and allows a direct comparison. We note that we do not compare the results to experiments as this is not our aim: this type of comparison has been performed extensively elsewhere.

We have chosen various solid systems for our tests: three of the elemental semiconductors (carbon, silicon and germanium, giving an excellent sampling of different gap sizes); simple oxides (two polymorphs of SiO₂ and MgO); perovskite oxides, including materials with semi-core states (SrTiO₃, PbTiO₃ in the cubic phase and the distorted perovskite MgSiO₃); an example of a metallic system, non-magnetic bcc iron; and two weakly bonded systems (the ordinary form of ice, ice XI, which has hydrogen bonds; and a layered material, hexagonal boron nitride). These systems represent different forms of bonding and environment, and offer a good test of the different basis sets sizes and truncation approaches.

The calculations used either the PBE functional⁴⁶⁾ (Fe, SiO₂, MgO, MgSiO₃, ice XI and BN) or the PBEsol functional⁴⁷⁾ (C, Si, Ge, SrTiO₃ and PbTiO₃); the differences between the functionals are relatively small, though PBEsol is considered generally better for solids. The same functional has been used for all calculations of the same system (both QE and conquest). We use the LibXC library⁴⁸⁾ for PAO generation and conquest calculations.

conquest uses a regular grid for certain integrals, which is set using an energy cutoff. In these calculations the energy has been converged with respect to this grid. In plane wave calculations the charge density grid is generally taken to have an energy cutoff four times that of the plane waves, and we note that the conquest grid energy cutoff is often around four times the QE plane wave cutoff. Details of cutoffs and Brillouin zone sampling (using Monkhorst–Pack grids) are given in the table captions.

For the elemental semiconductors, shown in Table I, we see good agreement with DZP basis sets (differences up to 2%–3% of the bulk modulus and 0.5% of the lattice constant) and excellent agreement with TZTP basis sets (differences typically less than 1% of bulk modulus and 0.2% of lattice constant). For context, when plotting binding energy curves (energy versus volume), a difference of 1% in the bulk modulus can barely be distinguished (after correcting for differences in lattice parameters). We illustrate the worst agreement, for Ge, in Fig. 1 and the PAOs generated with equal radii (note that all the curves are plotted relative to the fitted minimum volume and energy). In this case, the TZTP basis set gives very slightly worse agreement than the DZP, which is almost indistinguishable from the fully converged plane wave result. The SZP is a little too soft, but is still respectable. There is little difference between the equal radii and equal energy basis sets, except at SZP where the equal energy approach is maybe a little better.

Zoom In Zoom Out Reset image size

The simple oxides, with data found in Table II, show similar performance to the elemental semiconductors. The TZTP basis sets are converged to within 1% of the bulk modulus for the equal radii, but to within 2% for the equal energies, while the equilibrium volumes are all with 1% (for both α-quartz and stishovite the unit cells are non-cubic, so we use equilibrium volume rather than lattice constant).

A more challenging test for default PAOs (rather than PAOs tuned or optimised to an environment) is to compare the equilibrium phase (alpha quartz) to a high pressure phase (stishovite). The coordination of the Si atoms changes from four to six between these two phases. We compare the energy difference per formula unit (i.e. per SiO₂ unit) for the basis sets to the converged plane wave result (−0.28 eV/unit). The TZTP basis sets give a good comparison (−0.20 eV/unit for equal energies and −0.22 eV/unit for equal radii). The DZP basis sets using perturbative polarisation (the default setting) give the correct ordering (−0.08 eV/unit for equal energies, and −0.05 eV/unit for equal radii) if a less accurate magnitude; however, generating l = 2 orbitals simply as excited eigenstates of the confined atom gives the incorrect ordering for both DZP (+0.15 eV/unit for equal energies and +0.19 eV/unit for equal radii) and SZP (+1.46 eV/unit for equal energies and +0.52 eV/unit for equal radii; these values are improved with perturbative polarisation but still have the incorrect sign). To test the dependence further, we added a second polarisation function (giving DZDP), which has an excellent result: −0.19 eV/unit for equal energies, and −0.20 eV/unit for equal radii (for perturbative polarisation functions; non-perturbative functions give a similar result). Clearly comparisons of stability of structures require radial flexibility in all angular momentum channels, and we would recommend use of at least DZDP whenever considering this kind of problem with default (i.e. not optimised) basis sets. For structural properties of the individual phases, however, the performance of the DZP basis sets is reasonable, though not quite as accurate as for the elemental semiconductors.

We note that, in general, significantly different environments will require reasonable basis set flexibility as we have just shown for silicon dioxide. This may become more of an issue when simulating large-scale problems, as conquest is designed to do, and it will be important to explore the effect of basis set size in large systems with areas of different structure. While it is possible to tune a basis set to reproduce specific physical quantities with respect to a reference structure, these basis sets will be less transferable than large, general basis sets, and will require careful use in extended systems.

The perovskite structures, shown in Table III, feature elements with semi-core states (4s and 4p states in Sr and 3s and 3p in Ti; the 5d states for Pb might be considered semi-core but are not in this case) which are described with a single radial function. MgSiO₃ combines the elements seen in the simple oxides (we note that the pseudopotential in use here includes the 2s and 2p states, and that they are treated as semi-core). The performance is excellent for PbTiO₃ and good for SrTiO₃ (though the bulk modulus for TZTP with equal radii is surprisingly inaccurate) and the equilibrium volume for both these materials is very close to the plane wave result. It is possible that treating the Sr 4p states as valence, with further radial functions, might improve the performance. The bulk modulus for MgSiO₃ is significantly worse than either MgO or SiO₂, with errors of nearly 5%. Nevertheless, these results give confidence in the default basis sets.

Metallic bonding is very different to the covalent and ionic bonding studied thus far, but the performance of the default basis sets for non-magnetic bcc Fe, shown in Table IV, is excellent (we chose non-magnetic Fe simply for convenience; conquest is capable of spin-polarised operation as simply as non-spin-polarised, but a magnetic calculation would provide no more of a test of the PAO basis set than a non-magnetic one). In this case, both TZTP basis sets reproduce the plane wave results (the equal energy case has errors of 0.5% while the equal radii 2% in bulk modulus) while the DZP give excellent results. The equal energy SZP is still reasonable, though the equal radii is a little inaccurate. (For the Fe atom, both 3s and 3p states are included as semi-core states.)

Weakly bound systems offer a larger challenge to local orbital basis sets (putting to one side the issues that DFT has with these systems, which will be the same for any other basis set). We start with hydrogen bonding, considering the optimum unit cell for ice XI, the ordinary form of ice, shown in Table V. With a TZTP basis set, the three parameters are all accurate, comparing to the PW result to better than 1%, while DZP and SZP are much less accurate. This is the first case we have found where the use of a default (non-optimised) DZP basis set might cause a significant error. It is interesting to see how effective the TZTP basis sets are, even for a difficult system like this.

Finally we turn to the layered material, boron nitride, which has dispersion interactions between layers, where the range of the PAOs may play a key role in the DFT part of the interaction; for the conquest calculations, counterpoise corrections were used⁴⁹⁾ to account for basis set superposition errors. Data for the minimum distance between layers, and the resulting interaction energy is shown in Table VI. We report results for both standard PBE and PBE with semi-empirical dispersion corrections;⁵⁰⁾ again, we note that we are testing the accuracy of our default PAO basis sets against converged PW calculations (and not the accuracy of DFT).

Good accuracy in comparison to the PW results is achieved with DZP and TZTP for the equal radii and equal energy basis sets for both approaches, though the equal energy DZP basis is perhaps a little high in the PBE interaction energy. The SZP basis set is less accurate, though it is interesting that the inter-layer spacing is rather close with the PBE-D2 approach, suggesting that the dispersion corrections mask the basis set errors.

We have introduced an approach to the construction of PAO basis sets for the conquest code based on increasing energies (spaced by factors of ten). The radial functions are either based on a geometric increase in energy (equal energies) or on equal radii based on the average of the energies for all angular momenta (equal radii). These are both implemented in the conquest PAO generation code which will be made available with the conquest code on its release.

We have compared the performance of small (SZP), medium (DZP) and large (TZTP) basis sets with exact diagonalisation to converged plane wave results using the same functionals and pseudopotentials. Overall, the large (TZTP) basis sets are shown to reproduce the plane wave results to better than 1% in bulk modulus and often within 0.1% of the lattice constant for a variety of solids. Except for weakly bonded systems (particularly ice) we found that the medium (DZP) basis sets are nearly as accurate as the large basis sets and will typically be 5–10 times faster (simply from the reduction in matrix size). In general, both approaches (equal energies and equal radii) produce reliable results, though the equal radii approach gives slightly better results overall, and with smaller overall radii will give better efficiency. Alongside results from codes such as FHI-aims,⁵¹⁾ these results demonstrate that local orbital codes can be as accurate as converged plane wave codes.

This work was supported by World Premier International Research Centre Initiative (WPI Initiative) on Materials Nanoarchitectonics (MANA), New Energy and Industrial Technology Development Organization of Japan (NEDO) Grant (P16010), "Exploratory Challenge on Post-K computer" by MEXT, and JSPS Grant-in-Aid for Scientific Research (18H01143, 17H05224). The authors are grateful for computational support from the UK Materials and Molecular Modelling Hub, which is partially funded by EPSRC (EP/P020194), for which access was obtained via the UKCP consortium and funded by EPSRC grant ref EP/P022561/1. They also acknowledge computational support from the UK national high-performance computing service, ARCHER, for which access was obtained via the UKCP consortium and funded by EPSRC grant ref EP/K013564/1.