Results of the DeltaCodes GGA comparison exercise
FP Verification
Table of Contents
Introduction
The DeltaCodes project is an ongoing effort to test the mutual consistency of different implementations (code projects) of density-functional-theory; the first major success of the project has been the formulation of a standardised test set, consisting of 71 elemental crystals, together with the definition of a simple parameter for assessing the similarity of the results of different codes. These tests compare the equilibrium volume and bulk modulii of the different cases as calculated by different codes using the popular PBE generalised gradient approximation for exchange and correlation. A large number of codes have participated in the verification exercise, representing a considerable community engagement, and highlighting the growing importance of accuracy and reliability in the use of DFT in solid state physics. This large community effort is coordinated by Kurt Lejaeghere and Stefaan Cottenier DeltaCodes at Uni Ghent – details of the results may be found in the recent Science article.
The following presents the results of the full-potential lmf code for the DeltaCodes exercise. There are two data sets corresponding to the use of an linearized-muffin-tin-orbital with local-orbital (LMTO+LO) basis or a plane-wave plus LMTO basis (PMT documentation). The LMTO basis is highly suited for the treatment of dense solids, but is less satisfactory for highly open structures such as molecular crystals where the atom-centred basis functions are not always adequate for describing the wave function in the large interstitial region. Addition of plane waves to the basis fulfils this requirement satisfactorily, such that the full-potential lmf code may be considered “general purpose.”
Procedure
Selection of LMTO Basis Parameters
The highest l defining the basis, LMX, is taken from an internal table (by the blm setup program) and is generally given by the one more than the highest l channel occupied in the atomic case – such low values for the angular momentum expansion are only possible because of the efficient three-fold potential representation. Increasing LMX increases the size of the basis significantly, but is necessary if the highest accuracy is required. All calculations are based on a 2-kappa LMTO basis, consisting of two basis functions per l,m. This choice of basis is activated by the control flags HAM_AUTOBAS_LMTO=5 and HAM_AUTOBAS_MTO=2. The tails of the basis are Hankel functions with energies chosen automatically in the “GW-style,” specified by HAM_AUTOBAS_GW=1, with EH=-0.3 and EH2=-1.1Ryd. The smoothing of the Hankel functions at the atom centres is an important part of the basis definition and is specified by the RSHM parameters. The smoothing radii are determined automatically for the free atom case, depend upon the sphere radius and are different for different l. The automatically generated RSHM values are used without modification for all volumes.
Local Orbital Parameters
Extended local orbitals are included in the basis for all elements where the atomic state lies higher than -2.5Ryd or where the charge beyond the augmentation sphere is larger than 0.002. These conditions capture more semicore states than the default settings, specifically the semicore 3p states in Ca and Sc (where the atomic states lie at -2.06 and -2.47Ryd) and the transition metals Ru, Rh, Ir and Pt, whose the 4p and 5p states extend significantly into the interstitial.
case: Ca (atomic 3p @ -2.06Ryd)
V0[AA^3/at] B0[GPa] B1
none 43.28763 18.16842 3.507 (rather large)
PZ= 0 13.939 42.31310 17.57103 3.392
case: Sc (atomic 3p @ -2.47Ryd)
V0[AA^3/at] B0[GPa] B1
PZ= 0 0 4.4 25.00486 57.46736 3.567 (rather large)
PZ= 0 13.94 4.4 24.63446 54.95447 3.390
The transition metal block elements are significantly improved by the inclusion of high lying local orbitals for the d- channel; these improve the quality of the basis (essentially of the augmentation) with only modest increase in numerical expense. The special cases Hf, Ta and W are the only transition metals where no such high lying local orbitals are included, because they cause some instability in combination with the very low 4f semicore local orbitals; this leads to a small reduction of the quality of the basis for these elements (which is reflected in their DeltaCodes values). The importance of including both the high lying and semicore local orbitals is clear for the example of osmium:
case: Os (**REL=11**)
**LMX** **PZ** V0[AA^3/at] B0[GPa] B1
3 PZ= 0 15.942 14.34250 399.85174 4.904 (rather large)
3 PZ= 0 0 6.4 14.33800 402.36085 4.911 (rather large)
3 PZ= 0 15.942 6.4 14.29217 396.91855 4.847
4 none 14.33470 404.14945 4.899 (rather large)
4 PZ= 0 15.942 14.33686 401.40625 4.892 (rather large)
4 PZ= 0 15.942 6.4 14.28739 398.77032 4.842
By default, the PZ of extended semicore local orbitals is floated during iteration toward self-consistency, while P for the valence states of the same l is fixed. We observe that this procedure can lead to instability in some situations (notably Cr) and a reduction of accuracy in general. Consequently, the floating of all semicore local orbitals is deactivated using IDMOD=1 for the corresponding l.
Augmentation
Within the atomic spheres, the basis functions are augmented with solutions to the radial Schroedinger equation. These are solved for the valence states in the presence of a core charge density that is constant, having been derived in an atomic calculation, but which is allowed to extend into the interstitial. Scalar-relativistic effects are included in the solution of the radial problem. Although a number of augmented plane-wave all-electron codes in the DeltaCodes reports solve the radial Dirac equation for the core states, which includes also the spin-orbit interaction, the relatively small charge of the core that is typical in LAPW (a much larger set of states are included in the valence than is typical in LMTO) means that the effect of the inclusion of a fully-relativistic core is minimal. For the LMTO method, this is not quite true, and for consistency with the DeltaCodes tests, the core and the valence states are treated at the scalar-relativistic level only.
In the LMTO method, solutions of the radial problem, together with their energy derivatives, form a basis in terms of which the true solution may be expressed. Correspondingly, the individual augmentation functions need not be exact solutions to the muffin-tin potential and – in order to have an identical basis for all volumes – the augmentation functions evaluated at one volume may be used unchanged in calculations at other volumes. This is done using the –wratrho and –rdatrho command line arguments to lmf. The rationale of fixing the augmentation is to ensure that changes in the basis do not give rise to errors in E(V); in practice, fixing the basis in many cases slightly degrades the agreement with other codes or else give rise to numerical instabilities in the calculation. We find that using augmentation obtained at one volume for other volumes gives generally consistent results to the non-fixed case, while using augmentation functions calculated for the atomic problem (which corresponds in some respects to the pseudopotential approach) gives rather different results, indicating the importance of solving the radial problem for the augmentation simultaneously with the band problem.
Integration of the radial equations proceeds on a shifted logarithmic grid with A=0.01, increasing the mesh density.
case: Ag
E(A=0.01)-E(A=0.025)=0.3mRyd/atom
E(A=0.008)-E(A=0.01)=0.005mRyd/atom
Choosing Augmentation Sphere Radii
In order to maintain a constant quality of augmentation, sphere radii are fixed during variation of the volume. The augmentation sphere radii () are chosen to be 98% of the touching sphere radius at the equilibrium volume V_0. The scaling ensures that the spheres are all-but-touching also at 0.94V_0, as required for the tests, although the results are insensitive to small overlap. For some elements with very large lattice spacing, for instance the noble elements, this procedure gives unreasonably large , and the radial Schroedinger equation becomes difficult to solve exactly; in these cases, is restricted to 3.0bohr without significant loss of accuracy. For large , failure of the radial solver may happen during the self-consistency cycle. In practice, allowing the sphere radii to change and having touching-spheres at all volumes has only a small impact on the quality of the results (provided that the maximum ~3.0).
Determining GMAX and the k-Sampling
Using a trial k-point mesh (using a heuristic length rule, or using published meshes), the GMAX parameter is iterated upwards for each element at the equilibrium volume until the converged total energy changes by less than 1e-5Ryd. A regular k-point mesh is used in all the calculations with a number of points along each direction proportional to the length of corresponding reciprocal-lattice vector. The total number of k-points is decided by successively sampling, using the converged GMAX parameter, a range of values (increasing by a factor >1.15 at each step) until the change in total energy is less than 1e-5Ryd. In determining the GMAX and k-point meshes, an finer total energy convergence criteria of 2e-6Ryd is used. Scaling of the GMAX parameter with lattice constant, in order to ensure a constant number of grid points throughout variation of the volume, is unnecessary provided that GMAX is converged at V_0.
Inclusion of Plane Waves
The majority of the elements in the test set form close packed structures that are extremely well represented by the LMTO+LO basis; the remainder benefit substantially from the inclusion of plane waves in the basis. These replace the traditional use of empty spheres (floating orbitals) in the LMTO method by providing a set of non-atom-centred functions but which have the advantage of being simple and effective for arbitrarily complex geometries.
The value of NKABC and GMAX determined for the LMTO case is used. As the plane wave basis is made bigger, the order of matching at must be increased, in recent builds of lmf this is automatically taken care of, and the variable KMXA is updated in line with increases in PWEMAX.
Other Settings
To increase the precision of the calculations, the HAM_TOL criteria is set to 1e-16 and the Ewald summation tolerance ( EWALD_TOL ) is set to this value also. For evaluation of the volume dependence of the energy, the total energy is converged to better than 1e-5Ryd at each volume. The tetrahedron method (which is the default) is used for all k-integrations.
Choice of Reference
The mutual agreement of different all electron and PAW schemes demonstrated by the Science article shows that choosing one or another different code as a reference will lead to largely similar conclusions. In particular, if any calculation deviates strongly from the result of one code, it will also deviate considerably from the results of a second code. This is particularly true between the all-electron schemes with which the full-potential LMTO code should be compared. For simplicity and in accord with earlier incarnations of the DeltaCodes tests, the results are graphically represented here in comparison to the results of WIEN2k; a comparison with the other codes is straightforward using the data and reference information available from the DeltaCodes authors.
LMTO+LO Results
The performance of the LMTO+LO method is illustrated below for the close-packed elements, where the Delta values in comparison with the WIEN2k results are plotted. The computational parameters to reproduce these values are detailed LMTO settings, the fitted equation-of-state parameters can be read here, and the tabulated Delta values can be read here. A clear degradation of the LMTO results is clear for group V and group VI elements. For these cases, the inclusion of additional high lying (p) local orbitals has little effect.
PMT Results
The PMT method is very effective at remedying the shortcomings of the LMTO basis for open and intermediate packing systems, but is unnecessary for close packed (HCP, FCC, BCC) cases. Plane waves and LMTO basis functions are not orthogonal and, for close packed cases the overlap can make the basis set overcomplete – this causes real practical problems for the inclusion of plane waves in LMTO calculations and a consistent strategy is needed in defining when and how to use them.
PMT strategy 1:
PWEMAX, the plane wave energy cut-off, may be iterated upwards and the total energy converged in terms of this parameter. This can be done on top of the and basis chosen for the LMTO basis, such that the plane wave part can be considered as an extension to the LMTO method. For a modest convergence of the total energy (to 1mRyd) this may be possible with no further modification, although for the open systems the PW basis may become very large and expensive: for example nitrogen where touching spheres (the interatomic distance is 2.1bohr, so =1.05bohr) occupy only 2.4% of the cell volume, the convergence of the total energy with the plane-wave basis can be quite slow. The convergence of the equation-of-state parameters, however, proceeds much more quickly.
case: N
**PWEMAX**[Ryd] Etotal[Ryd]
LMTO: -875.591611
0.5 -875.834779
1.0 -875.949097
1.5 -876.006302
2.0 -876.035862
4.0 -876.086631
6.0 -876.114569
8.0 -876.134762
10.0 -876.150119
12.0 -876.161751
14.0 -876.170625
PMT strategy 2:
An alternative strategy is to progressively replace the LMTO basis with the PMT one, and removal the second set (and, possibly higher-l components) of LMTO basis functions. With a larger PW component, numerical difficulties emerge because of linear dependence in the basis requiring some care in the specific choice of setup: often it is helpful to freeze the augmentation functions and to use a considerably reduced augmentation radius. These steps can be followed to converge the total energy more precisely than is possible using the LMTO method alone, even with optimised basis parameters, but the importance of the different approaches, and their effect on numerical stability and expense, is quite system dependent; this added complexity is not generally necessary for physical properties such as E(V) but may be very valuable for response-function calculations () where the larger basis affords a much better representation of empty states.
PMT strategy 3:
Keeping the large and LMTO basis, the most pragmatic and effective use of the plane wave basis is to attempt a small series of different PW sets, with PWEMAX=0,2,5Ryd, until the desired property is well represented. To enable this, the basis can be checked for overcompleteness by diagonalising the overlap matrix and discarding basis functions corresponding to eigenvalues smaller than a given threshold (OVEPS). We find OVEPS=1e-10 is generally suitable, although it must be remembered that when OVEPS is in use, the basis is no longer exactly defined by PWEMAX and may change during the self-consistency cycle. In the figure, a cut-off of 2Ryd is used for all elements except N,C,O where a larger cut-off of 7.5Ryd is necessary.
Conclusions
Very satisfactory results can be obtained for close packed structures using the default LMTO basis set when local orbitals are included for semicore states and, in the case of the transition elements, high-lying local orbitals are included in the basis. Agreement with other codes presented in the DeltaCodes reports is shown using LMTO basis set parameters (particularly , RSMH, EH, and PZ) that are calculated entirely automatically without intervention and thus represent the accuracy that would be obtained more or less “out of the box”. We note, of course, that the LMTO basis is actually quite sophisticated and increased accuracy can be obtained by optimising the basis – the lmf code contains utilities to conduct such an optimisation, which should be investigated for each study independently. For open systems, the LMTO basis can be generally improved using the combined LMTO/plane wave approach: the addition of a small set of plane waves, with cut-off 2Ryd, is shown to correct the description of the more open non-metal structures without the need for empty spheres. Depending upon the particular study, the relative importance of the MTO and PW parts of the basis in describing the valence states can be changed, which is of particular benefit for calculations which rely upon calculation of response functions via a sum over states.