An efficient procedure for the development of optimized Projector Augmented Wave basis functions
aa r X i v : . [ c ond - m a t . m t r l - s c i ] A ug An efficient procedure for the development of optimized Projector Augmented Wavebasis functions
R. J. Snow, A. F. Wright, and C. Y. Fong Department of Physics, University of California, Davis, California 95616, USA Sandia National Laboratories, Albuquerque, New Mexico 87185-1415, USA (Dated: July 21, 2010)In the Projector Augmented Wave (PAW) method, a local potential, basis functions, and pro-jector functions form an All-Electron (AE) basis for valence wave functions in the application ofDensity Functional Theory (DFT). The construction of these potentials, basis functions and pro-jector functions for each element can be complex, and several codes capable of utilizing the PAWmethod have been otherwise prevented from its use by the lack of PAW basis sets for all atoms. Wehave developed a procedure that improves the ease and efficiency of construction of PAW basis sets.An evolutionary algorithm is used to optimize PAW basis sets to accurately reproduce scatteringproperties of the atom and which converge well with respect to the energy cutoff in a planewavebasis. We demonstrate the procedure for the case of Ga with the 4s, 4p, and 3d electrons treatedas valence. Calculations with this Ga PAW basis set are efficient and reproduce results of linearizedaugmented plane wave (LAPW) calculations. We also discuss the relationship between total energyconvergence with respect to the energy cutoff and the magnitude of the matching radius of the PAWset.
I. INTRODUCTION
The Projector Augmented Wave (PAW) method ofBl¨ochl provides an efficient, all-electron (AE) basisfor Density Functional Theory (DFT) calculationswithin the frozen-core approximation. The method isapplicable to all atoms, and the efficient treatment offirst-row and transition metal elements compared withnorm-conserving pseudopotential methods is improveddramatically . The PAW method requires the construc-tion of a local potential, atom-centered radial (pseudo-and AE) basis functions, and projector functions for eachatom. For PAW functions highly optimized with regardto computational efficiency, this construction can be com-plex, including a local potential and multiple sets of ba-sis and projector functions for each angular momentumchannel and principal quantum number. To acceleratethe construction of these PAW functions, and, further,to optimize with respect to the energy cutoff, we havedeveloped an efficient, partially automated procedure forsearching this high-dimensional parameter space.In section II we introduce the basic elements of thePAW and its associated parameters. Our procedure forobtaining efficient ab initio PAW functions is developedin section III. In section IV results of the constructionof the Gallium (Ga) PAW are given including logarith-mic derivatives, total energy convergence, optimizationresults, and lattice constant and bulk modulus calcu-lations. A primary key to obtaining optimal total en-ergy convergence is a large matching radius, and this isshown by the comparison of a series of Ga PAW basis setswith increasing matching radii. A moderate variation oflattice constant with the matching radius parameter isfound, and this will be discussed later. PAW basis setsconstructed in this way are expected to accurately and ef-ficiently calculate the properties of solids within the DFTframework and the frozen-core approximation.
Radius [Bohr] -15-10-50 P o t e n ti a l E n e r gy [ R y ] AE potentialPAW local potential
Gallium AE and Local Potentials R PAW
FIG. 1. (Color online) Ga AE and PAW local potentials.
II. ELEMENTS OF THE PAW CONSTRUCTION
The construction of PAW basis sets is similar to theconstruction of pseudopotentials. First, a reference con-figuration is chosen and electrons are divided into coreand valence electrons. A local potential is created whichis smooth within the PAW matching radius R P AW andmatches the AE potential at and beyond R P AW , as in fig-ure 1. Radial basis functions and dual projector functionsused as the basis for electronic wave functions are con-structed for each angular momentum channel such thata projection operation of ˜ p on a smooth (pseudo) basisfunction ˜ φ reproduces the AE basis function, φ , as in fig-ure 2. A single basis function per angular momentumchannel is often sufficient, but more regularly two func-tions per channel are required to approach a completebasis set , and in rare cases multiple basis functions inone angular momentum channel but with different prin-cipal quantum numbers may also be needed to representthe wave functions of semi-core electrons.We use the ’atompaw’ program to construct basis setsin an RRKJ (Rappe, Rabe, Kaxiras, and Joannopou- -2024 L=0 φφ ∼ p ~ b -2024-2024 L=1 -20240 0.5 1 1.5 2 2.5 3
Radius [Bohr] -6-4-2024
L=2
Gallium Projector and Basis Functions
FIG. 2. (Color online) Ga Projector and Basis Functions. Foreach angular momentum channel, two basis functions may beused. The first PAW AE basis function, pseudo basis function,and projector function are shown in column a and the secondset of functions is in column b. los) style , one of several available in the ’atompaw’program. The RRKJ scheme optimizes the kinetic en-ergy because this dominates the total energy conver-gence. With RRKJ pseudo- basis functions and projec-tors, each projector function is associated with an indi-vidual matching radius, R L,i , and an energy E L,i . Theenergy of the first pseudo basis function in each angu-lar momentum channel is taken as the eigenvalue of thecorresponding AE solution in the atomic system, leavingthree free construction parameters per angular momen-tum channel when two projectors are used.A typical PAW will have 3 ∗ N L + 2 construction pa-rameters, including energy and matching radius param-eters for the local potential, where N L is the numberof angular momentum channels used in the basis, usu-ally corresponding to the angular momentum channelsof the valence electrons. These parameters are chosen tooptimize the matching of logarithmic derivatives and thetotal energy convergence with respect to an energy cutoffof the planewave basis when using these PAW functionsto perform calculations of solids. Basis functions associ-ated with separate angular momentum channels appearto be independent with respect to logarithmic deriva-tives. However, in a crystalline calculation, these werefound to have a strong interdependence in their effect ontotal energy convergence and in the self-consistent en-ergy minimization. Therefore, with a typical resolutionof 0.05 Bohr in the matching radii and 0.1 Rydbergs inenergy, and total ranges of 1 Bohr in the radii and 10 Ry-dbergs in energies, this creates an entire parameter spaceof 20 N L +1 ∗ N L +1 total number of possible parame-ter sets, or about 1 . ∗ sets for elements with s, p,and d basis functions. The size of this parameter spaceprecludes a purely brute-force attempt to optimize PAWprojection basis sets. In addition to previous advice inthe literature and in various groups regarding the selec-tion of construction parameters, below and in section III we discuss further ways that this parameter space can bereduced and efficiently searched.The matching of logarithmic derivatives of the pseudo-and AE wave functions, evaluated with respect to theradius at R P AW as a function of energy, represents anequivalence of the scattering properties of the PAW basisset and the atom. The matching of logarithmic deriva-tives over a wide range of energies is one indication of atransferable basis set. In this ab initio construction ofPAW basis sets, the primary aim is to accurately repro-duce logarithmic derivatives , and we do not optimizethe sets by the matching of results of calculations withany observable.Traditionally PAW parameters are chosen by hand,logarithmic derivatives are inspected visually, and post-testing of basis sets for optimal total energy convergence,and eventually for completeness by comparing with otherAE results, is done individually for each constructedPAW. This may take many iterations to discover well-matching logarithmic derivatives which are also optimalfor total energy convergence. Also, much of parame-ter space invariably remains unexplored due to its pro-hibitive size. In this paper we propose an approach to re-duce the time and labor associated with producing PAWbasis sets and to discover otherwise unobtainable opti-mized PAWs. In our procedure, first logarithmic deriva-tives are given a visual inspection to define a limitedrange of possible construction parameters, thereby reduc-ing the available parameter space significantly. Then anevolutionary algorithm is used to further optimize boththe matching of logarithmic derivatives and the total en-ergy convergence of a PAW. In some cases, it is nec-essary to also optimize with respect to the number ofself-consistent iterations. III. METHOD
The method described in this paper for obtaining anoptimized PAW basis set consists of, first, the selectionof configurational parameters, including the division ofelectrons into core and valence electrons and the selec-tion of the basic types of local potentials and basis func-tions. Next, parameter ranges are reduced through avisual inspection of logarithmic derivative matching andbasis function smoothness. Then, an evolutionary al-gorithm further optimizes the PAW within the reducedparamter space. Lastly, a few tests are performed to en-sure the basic accuracy and tranferability of the PAW.Further testing of the PAW may be necessary dependingon the intended calculation .Initial parameter ranges can be restricted as follows. Alocal potential matching radius R P AW should be as largeas possible without resulting in sphere overlap in crys-talline, molecular, or molecular dynamics calculations ,in order to optimize the total energy convergence, as willbe discussed in section IV.5. This is usually a little lessthan half the nearest neighbor distance in a particularsystem, with room to allow for ionic motion or relax-2tion. In our case, a Troullier-Martins style local po-tential is used, and the energy parameter often rangesfrom -2 to 5 Rydbergs, while energy parameters for theRRKJ projector functions have ranged from -2 to 8 Ry-dbergs. After the local potential matching radius R P AW is determined, all individual matching radii R L,i are setequal to this. Then a range for the local potential energyis determined by viewing the logarithmic derivative forthe L = L max + 1 channel. With a suitable local po-tential energy chosen in this range, energies of the pro-jector functions are varied simultaneously to find regionswith good logarithmic derivatives, but with a secondarygoal of obtaining smooth basis functions as well. If well-matched logarithmic derivatives are unobtainable withonly these energy parameters, either the R P AW or the R L,i may then also be altered to improve the logarith-mic derivatives. Ranges of parameters which yield well-matched logarithmic derivatives are then used within theoptimization program. In the case of Ga, E s ranged from3.0 to 7.0 Rydbergs, E p ranged from 0.0 to 7.0 Rydbergs, E d ranged from 0.0 to 3.5 Rydbergs, and E loc rangedfrom 0.0 to 5.0 Rydbergs, leaving an effective parameterspace for Ga of about five million sets.In order to automate the process of optimizing log-arithmic derivatives and total energy convergence, wedefine fitness scores for both of these properties whichwill be used in an optimization algorithm. For this pur-pose a numerical comparison of logarithmic derivatives isevaluated for each potential PAW set, according to thesummation, P ( y i,AE − y i,P AW ) / ( exp ( − ABS ( dy/dx )),where in this expression y i,AE and y i,P AW are the all-electron and pseudo- wave function logarithmic deriva-tives as functions of energy, respectively. The purpose ofthe denominator here is merely to attenuate the effect ofthe difference between the AE and PAW values near adivergence in the logarithmic derivatives, as can be seenfor example in figure 3. Scores for each angular momen-tum channel are summed and typically normalized by anaverage ideal matching score.If logarithmic derivatives are not within a reasonabletolerance (of about 0.3), a penalty factor is given as thefitness score and returned to the optimization routinewithout further testing. Ghost states may be de-tected by the presence of a divergence in the PAW log-arithmic derivative where there is none in the AE loga-rithmic derivative, and this case is given a large penaltyfactor. If logarithmic derivatives are free of ghost statesand the matching of PAW and AE derivatives is withintolerance, crystalline tests at two successive energy cut-offs are performed, and their energy difference, normal-ized by an ideal difference (about 0.01 Ry), is used asa convergence score. Logarithmic derivative matchingis emphasized with stepped weighting, using a factor of10 for the weight function until within an ideal toler-ance (about 0.03) at which point the weight is reducedto 1, equivalent to the weight of the total energy conver-gence scores. Care must be taken so that the test withlarge energy cutoff is well-converged, in order to preventa plateau from forming in the total energy convergence figure as a function of energy cutoff. The Socorro DFTprogram with PAW functionality is used for these tests.The normalized logarithmic derivative and total energyconvergence scores are then summed to produce a finalfitness score.The Design Analysis Kit for Optimization and Teras-cale Applications (DAKOTA) offers convenient han-dling of input parameters, a parallel computation frame-work, and an interface to several optimization techniquesand packages. Comprehensive testing of DAKOTA algo-rithms for searching PAW parameter space, or perhapsnew algorithms, could reduce the time requirements fordiscovering optimal PAW functions, but we find that theevolutionary algorithm of the coliny package is suffi-cient. In the coliny ea algorithm, after a random groupof initial PAW parameter sets are evaluated, future pa-rameter sets are chosen as combinations of well-scoringsets from the previous iteration. In addition to the cross-ing of well-scoring parameter sets, random mutations arealso applied to original sets and to crossed sets, allowingfurther exploration into diverse areas in parameter space.In our practice, we typically use the ’two point’ methodof crossing with a ’cross over’ rate of 0.9, a mutation rateof 0.25 with the ’offset cauchy’ mutation method, and apopulation size of 100 and maximum number of itera-tions anywhere from 3,000 to 25,000, depending uponthe number of parameters involved. In the ’two point’method, optimal parent sets from the current generationare crossed by the division of parent sets into three re-gions, and the middle section is taken from one parentand the end sections from another parent. In the ’off-set cauchy’ method for mutation, a random number isgenerated according to a cauchy distribution with a meanof 0, and this is added to a selected parameter. Since theoptimization with DAKOTA is done with nominal test-ing, further tests of several of the top scoring PAW basissets are used to better gauge the total energy conver-gence and to verify the transferability of the PAW. Asa function of total energy cutoff, we look at the latticeconstant, bulk modulus, and total energy convergenceproperties to evaluate PAW basis set optimization. IV. RESULTSIV.1. Logarithmic Derivative Matching
The use of an evolutionary algorithm to optimize thenumerical matching of logarithmic derivatives results inmost cases with a nearly identical matching of loga-rithmic derivatives in the associated angular momentumchannels, a feat which is possible but painstaking whenoptimizing by hand. The case of the Ga PAW is shown infigure 3. Additionally, in many cases the matching of log-arithmic derivatives can be obtained without the use ofmatching radii R L,i as optimization parameters in eachangular momentum channel; these can often be set to thelocal potential matching radius, R P AW . The energy ofexpansion for a second projector function in each channel3
AEPAW
L=0 -4 -2 0 2 4-4-2024
L=1 -4 -2 0 2 4
Energy [Ry] -4-2024 d l n φ ( E ) / d r L=2 -4 -2 0 2 4-4-2024
L=3
FIG. 3. (Color online) The AE logarithmic derivatives ofthe Ga atom and its PAW approximation represent thewave function matching and the scattering properties of theatom. Excellent agreement in all angular momentum channelsL=0,1,2,3 is shown. then becomes the primary parameter. This reduces theparameter space significantly and eliminates the problemof kinks in projector functions, saving time both in thepreliminary inspection of logarithmic derivatives for therestriction of the parameter ranges and in the DAKOTAoptimization.
IV.2. DAKOTA Total Score Optimization
The final score used for a fitness evaluation is a com-bination of scores from the logarithmic derivative match-ing and from the total energy convergence, as discussedabove in section III. Figure 4 shows the scoring of atypical DAKOTA run of around 5,000 parameter sets forthe Ga PAW, using only E L,i with L=0,1,2 and E Loc as DAKOTA parameters while keeping all R L,i equal to R P AW . If an appropriate normalization for both the log-arithmic derivatives and the total energy convergence isused, the fitness score should converge to around 1.0. Inthis particular optimization, an ideal total energy conver-gence of 1 mRy was used for normalizing the total energyscore, the difference between the 25 Rydberg and 100 Ry-dberg cutoff calculations. An ideal logarithmic derivativematching of 0.01 was used to normalize the logarithmicderivative score. In this case, the total energy conver-gence normalization was too ambitious, and the optimalfitness score appears to converge to around 6.15, but theorder of magnitude is acceptable. The PAW with thebest score, known internally as Ga Rp2.1et001 5095 wasused in figures 1, 2, and 3, with all matching radii set to2.1 Bohr, E s found to be 4.287 Rydbergs, E p found tobe 4.647 Rydbergs, E d found to be 0.797 Rydbergs, and E loc found to be 2.248 Rydbergs. The standard devia-tion for the top ten scoring sets are 0.36, 1.79, 0.03, and0.09 for E s , E p , E d , and E l oc , respectively. PAW Iteration Number D a ko t a f it n e ss ( b e s t o f b i n s i ze ) Dakota optimization of Gallium PAW
Using zinc-blende GaN for E_total convergence
FIG. 4. DAKOTA total score optimization. The best score ofeach group of 100 parameter sets is plotted on the y-axis. Alower score represents better logarithmic derivative matchingand better total energy convergence.
IV.3. Variability of convergence and completenessof the Ga PAW
Total energy convergence rates and physical propertiesof PAW sets with the top ten DAKOTA scores in thistrial run appear to be nearly identical, as can be seen infigures 5, 6, and 7. This is due partly to a large measure ofsimilarity in parameter values of the optimized PAW sets,but also due to the effectiveness of the PAW method andan insensitivity of calculated properties such as latticeconstants to the exact parameters of the construction.
20 40 60 80 100
Energy Cut [Rydberg] T o t a l E n e r gy c onv e r g e n ce [ R ydb e r g ] ID 1666ID 3019ID 3643ID 3798ID 4249ID 4511ID 4965ID 4972ID 5005ID 5095
Gallium Total Energy Convergence
Top 10 results from Rp2.1et001, GaN with soft N
FIG. 5. Total energy convergence rates for the top ten resultsfrom a particular run are indistinguishable. Matching radii R L,i and R PAW are all equal to 2.1 Bohr. A soft NitrogenPAW is used here so that the total energy convergence isdetermined by the Ga PAW.
In the evaluation of the lattice parameter for GalliumNitride (GaN) in the zinc-blende structure, the R P AW parameter for Ga had a noticeable effect on the lat-tice constant, as can be seen in figure 8, with increasingmatching radii leading to smaller lattice constants, butwith a total range of only 0.5% in the lattice constant.We have similarly generated a series of N PAW sets withincreasing R P AW parameters, but this had negligible ef-fect on the GaN lattice constant.4
Energy Cutoff [Rydberg] L a tti ce C on s t a n t [ B oh r] ID 1666ID 3019ID 3643ID 3798ID 4249ID 4511ID 4965ID 4972ID 5005ID 5095
Gallium Lattice Constant Convergence
Top 10 results from Rp2.1et001, (zb) GaN with hard N
FIG. 6. GaN zinc-blende lattice constants as functions ofenergy cutoff are indistinguishable. Matching radii R L,i and R PAW are all equal to 2.1 Bohr. A hard Nitrogen PAW isused here so any variation in physical properties should beattributed to Ga PAW differences.
20 40 60 80 100
Energy Cutoff [Rydberg] B u l k M odu l u s [ M B a r] ID = 1666ID = 3019ID = 3643ID = 3798ID = 4249ID = 4511ID = 4965ID = 4972ID = 5005ID = 5095
Gallium Nitride Bulk Modulus
Top 10 results from Rp2.1et001, (zb) GaN with hard N
FIG. 7. GaN zinc-blende bulk moduli as functions of energycutoff are virtually indistinguishable. Matching radii R L,i and R PAW are all equal to 2.1 Bohr. A hard Nitrogen PAW is usedhere so any variation in physical properties can be attributedto Ga PAW differences.
Gallium Rpaw [Bohr] L a tti ce c on s t a n t [ B oh r] PAW (this work)LAPW Gallium Nitride Lattice Constant Variation
FIG. 8. The GaN zinc-blende lattice constant varies with the R PAW matching radius. The N PAW used in these calcu-lations has a matching radius of R PAW = 1.1 Bohr so thatthere is no overlap of spheres for any of the included points.A separate DAKOTA optimization was performed for eachvalue of R PAW . The maximum difference in lattice constantfrom R PAW =1.7 to R PAW =2.5 is less than 0.5%.
IV.4. Transferability
Transferability in the PAW method is effective dueto its AE treatment of wave functions and potentialswithin the core region, limited only by the complete-ness of basis functions, and by the inherent limitations ofthe frozen-core approximation and the DFT. To demon-strate the transferability of the Ga PAW, we have cal-culated lattice constants and bulk moduli for zinc-blendeGaN, zinc-blende GaAs, and zinc-blende GaP. The valuescalculated in the local density approximation (LDA) ,shown in Table I, are within 0.2% error in the latticeconstants in comparison with linearized augmented planewave (LAPW) values. Bulk moduli are within 1% error,except for zinc-blende GaP, which differs from LAPW values by nearly 2%. IV.5. R PAW as key to total energy convergence
While the E L,i parameters do affect total energy con-vergence, the dominant bottleneck parameter is the lo-cal potential matching radius, R P AW . For a series ofPAW basis sets with increasing R P AW , figure 9 showsa clear correspondence between the magnitude of thematching radius and the total energy convergence. Witha large matching radius, pseudo-basis functions and thelocal potential can be made smooth and therefore maybe expanded in a small number of plane waves. A smallamount of overlap of matching radii of atoms in a solidcalculation can have negligible effects, but in general thespheres should not overlap, limiting the size of the match-ing radius and the total energy convergence.
20 40 60 80 100
Energy Cutoff [Rydberg] T o t a l E n e r gy C onv e r g e n ce [ R ydb e r g ] Rpaw = 1.7Rpaw = 1.8Rpaw = 1.9Rpaw = 2.0Rpaw = 2.1Rpaw = 2.2Rpaw = 2.3Rpaw = 2.4Rpaw = 2.5
GaN Total Energy Convergence as a function of Gallium Rpaw
FIG. 9. (Color online) Correspondence between total energyconvergence and the matching radius R PAW . V. CONCLUSION
The procedure described in this paper automates theprocess of generating optimized PAW basis sets using an5
ABLE I. PAW calculations of zinc-blende lattice constants and bulk moduli of GaN, GaAs, and GaP were performed in theLDA with energy cutoffs of 100 Ry (although convergence is expected as in figures 6 and 7). Each of these used a single GaPAW, Ga Rp2.1et001 5095.Crystal PAW ID a [Bohr] B [MBar] LAPW a [Bohr] LAPW B [Bohr] a % error B % errorzinc-blende GaN N.febsec5953 8.439 2.03 8.447 2.05 -0.09 -0.9zinc-blende GaAs Arsenic d 1058 10.619 0.76 10.620 0.77 -0.02 -0.7zinc-blende GaP P.cea 22169 10.226 0.91 10.242 0.89 -0.16 1.89 a LAPW data from ’atompaw’ website . evolutionary algorithm to efficiently search a large pa-rameter space. In the example of the Ga PAW, the result-ing PAW matches all-electron scattering properties andis as efficient as the construction method and the partic-ular element permits. Efficiency and efficacy of the PAWwas confirmed in a handful of crystalline environments.This method will be of assistance in ongoing efforts toproduce efficient PAW sets for various DFT codes. VI. ACKNOWLEDGEMENTS
We acknowledge many discussions with NatalieHolzwarth and Marc Torrent, as well as for helpfulguidance given on the atompaw website . We wouldlike to thank Alan Tackett, Greg Walker, and RachaelHansel for an introduction to the DAKOTA program.Brian Adams helped in a critical way with details of theDAKOTA program. Sandia is a multiprogram labora-tory operated by Sandia Corporation, a Lockheed Mar-tin Company, for the United States Department of En-ergy’s National Nuclear Security Administration underContract No. DE-AC04-94AL85000. P. E. Bl¨ochl, Phys. Rev. B , 17953 (1994). P. Hohenberg and W. Kohn, Phys. Rev. , B864 (1964). W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965). G. Kresse and D. Joubert, Phys. Rev. B , 1758 (1999). N. A. W. Holzwarth, A. R. Tackett, and G. E. Matthews,Computer Physics Communications , 329 (2001),ISSN 0010-4655. M. Y. Chou, Phys. Rev. B , 11465 (1992). A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D.Joannopoulos, Phys. Rev. B , 1227 (1990). M. Torrent and N. Holzwarth,
A user’s guide for atompawcode (2007), . V. Heine, in
Solid State Physics , edited by F. Seitz andD. Turnbull (Academic Press, Inc, 1970), vol. 24. A. Kiejna, G. Kresse, J. Rogal, A. De Sarkar, K. Reuter,and M. Scheffler, Phys. Rev. B , 035404 (2006). X. Gonze, R. Stumpf, and M. Scheffler, Phys. Rev. B ,8503 (1991). X. Gonze, P. K¨ackell, and M. Scheffler, Phys. Rev. B ,12264 (1990). Socorro is developed at Sandia National Laboratories andavailable from http://dft.sandia.gov/Socorro/ . M. Eldred, S. Brown, B. Adams, D. Dunlavy, D. Gay,L. Swiler, A. Giunta, W. Hart, J.-P. Watson, J. Eddy,et al.,
Dakota, a multilevel parallel object-oriented frame-work for design optimization, parameter estimation, uncer-tainty quantification, and sensitivity analysis: Version 4.0users manual (2006), sandia Technical Report SAND2006-6337, Version 4.2, Updated November 2008. W. Hart,
The coliny optimization library (2004), http://software.sandia.gov/Acro/Coliny . J. P. Perdew and Y. Wang, Phys. Rev. B , 13244 (1992). N. Holzwarth et al., atompaw , ..