Author Archives: James M Goff

A brief introduction to the theory, parameterization and application of Projector Augmented Wave (PAW) datasets to plane wave DFT calculations

Leave a reply

Introduction

Many of the observable properties that makes materials interesting are linked to the electronic structure, so it is desireable with to perform DFT calculations with the most accurate computational methods that are affordable. Calculations within DFT scale as \( \mathcal{O}(n^3) \) for n electrons [1], so many technologically important materials, including transition metals such as Fe, are computationally demanding to calculate when considering all the electrons.This is exacerbated for core electrons, since the wave functions near the nuclei tend to oscillate, requring a large number of plane waves required to describe the wave functions near the nuclear cores in extended systems. Since many key properties are influenced by the outermost valance electrons that control bonding, it is sensible to fix the core electrons and treat them approximately to make calculations of the same level of theory more affordable or to make more complex calulations affordable. Treating core electrons as fixed is known as the frozen core approximation.

In this post, different electronic configurations and cutoff radii will be used to generate PAW datasets for Ti. Differences in the generated PAW datasets will be demonstrated by comparing the cutoff energy required to converge the energies of hcp Ti and TiO2 and by comparing the formation energy of TiO2 with calculated and experimental formation energies in the literature.

Treating core electrons

Core and valence electrons

Typically there are two main parameters to consider when treating electrons as core or valence that effect how the wave functions will be constructed: the core radius and the number of core vs. valence electrons. The core radius is the spherical distance from the center of the atom that forms the outer bound for the approximation to the wave functions. The number of electrons in the core vs. valance is simply the partition of how many electrons of charge density an atom can freely supply to the system. For example, Ti has 22 electrons that must be partitioned between core and valance. Both of these parameters are typically chosen by the person developing the new description.

A basic description of some comment methods for treating core electrons in extended systems follows.

Pseudopotentials

Pseudopotentials replace the core wave functions with a pseudo wave function with the goal of defining a smooth function to represent the effective potential of the core. The pseudo wave function are constructed so that the value and the derivative of the psuedo wave function matches the value and derivative of the true wave functions at the core radius.
The advantage of pseudopotentials is their formal simplicity, however the pseudo wave function prevents the recovery of the full core electron wave functions. James Goff published a post on pseudopotentials [2] that describes their construction in a detailed, introductory way.

LAPW

Another method for treating core electrons is the Linear Augmented Plane Wave (LAPW) method. The main idea of this method is that two partial waves, which can be used to describe the all electron wave functions, are joined are the core radius and matched to the all electron solution. Thus, the main advantage is the core wave functions can be recovered, but this method comes at a higher computational cost than the pseudopotential method. A more complete description of the psueodpotential and LAPW methods can be found in [3].

PAW

The Projector Augmented Wave (PAW) method presented by Blöchl [4] as a formalism that generalizes and combines favorable aspects of the pseudopotential and LAPW methods through projector functions. The basis of the PAW method is the combination of partial wave functions from the isolated atom with pseudo partial waves through projector functions that are constructed to match the all electron solution. In the PAW method, PAW datasets are generated that can recover the full core electron wave functions, like the LAPW method, but at a lower cost than the LAPW method, closer to the computational cost of pseudopotentials. More information can be found for the PAW method in [4-6] and more detailed description of the relationship between pseudopotentials, APW, and the PAW methods can be found in [4].

Calculation details

Calculations are performed using plane waves and the PBE functional implemented by Quantum Espresso [8] in the pwscf code. PAW datasets [4] are used to describe the core states. For Ti, PAW datastets are generated as part of this work with the ld1 code distributed with Quantum Espresso, based on the PS Library 1.0 parameterizations [9]. The PS Library 1.0 parameterization for O is used directly. The reference states for Ti and O are chosen as hcp Ti and the O2 molecule in a 10x10x12 box. The size of the O2 box was converged to within 0.1 mRy. Detailed discussion of the convergence of the box size of scope for this post.

The cutoff energy convergence of the different potentials are compared in the following sections below for Ti and TiO2. Oxygen cutoff energy was converged to 0.1 mRy at a wavefunction cutoff energy of 50 Ry and a kinetic energy cutoff of 200 Ry. For the calculation of formation energy, the maximum converged cutoff energy for each set of potentials is used across all structures.

The kpoint grids were converged to within 0.1 mRy are shown for each structure in the table below. Showing the detailed kpoints convergence tests are out of scope for this study. Even though the different PAW datasets can lead to different electronic structures and kpoints convergence, a single set of kpoints for each set of structures were used across all of the different PAW datasets for simplicity.

Generating PAW datasets for Ti

The ld1 program was used to generate the PAW datasets. The documetation can be found online [10]. The starting point for the tests here is the input file for generating a PAW dataset for the spin polarized, scalar relativistic Ti from the PS Library version 1.0.0 [9]. The input file for the PS Library starting point follows below. The settings here create a PAW dataset for Ti with Troullier-Martins pseudization [11] via Bessel functions. The the key inputs that we will change are the electronic configuration, `config` (and the corresponding psedized wave functions), and `rcore`, which is the matching radius in atomic units for the core charge smoothing. The bottom section describes the number of wave functions (6) that will be pseudized and which are described by the lines below. Taking the first line (3S) as the example, the first number is the principle quantum number (n=1, starting with 1 for the lowest s, 2 for the lowest p), then the angular momentum quantum number (l=0), followed by the number of occupying electrons, the energy used to pseudize the state (=0.00), the matching radius for the norm-conserving pseudopotential (=0.85), the matching radius for the ultrasoft pseudopotential (=1.30) and finally the total angular momentum (=0.0). Changing the norm-conserving or ultrasoft pseudopotential matching radii do not change the PAW dataset, since the radii are controlled by the input parameters above. Note that the energies for the unbound states (4p=4 Ry and unoccupied 3d=0.05 Ry) need to be non-zero.

&input
title='Ti',
zed=22.,
rel=1,
config='[Ar] 4s2 4p0 3d2',
iswitch=3,
dft='PBE'
/
&inputp
lpaw=.true.,
pseudotype=3,
file_pseudopw='Ti-PSL.UPF',
author='ADC',
lloc=-1,
rcloc=1.6,
which_augfun='PSQ',
rmatch_augfun_nc=.true.,
nlcc=.true.,
new_core_ps=.true.,
rcore=0.8,
tm=.true.
/
6
3S 1 0 2.00 0.00 0.85 1.30 0.0
4S 2 0 2.00 0.00 0.85 1.30 0.0
3P 2 1 6.00 0.00 0.80 1.60 0.0
4P 3 1 0.00 4.00 0.80 1.60 0.0
3D 3 2 2.00 0.00 0.80 1.50 0.0
3D 3 2 0.00 0.05 0.80 1.50 0.0

To demonstrate how the cutoff radius affects the calculation accuracy and expense (through the cutoff energy), two new PAW datasets will be generated with core radii of 0.4 a.u. (much smaller) and 2.0 a.u. (much larger). In addition, a different electronic configuration will be generated that does not pseudize the 3s or 3p core states. Note that these are note expected to be accurate, but instead to demonstrate how the parameters affect might affect calculations and impress the care and detail that should be given to generating proper potentials.

Energy cutoff convergence of Ti and TiO2

Formation energy of Ti and TiO2

References

[1] D.S. Sholl, J.A. Steckel, Density Functional Theory, John Wiley & Sons, Inc., Hoboken, NJ, USA, 2009. doi:10.1002/9780470447710.
[2] J. Goff, Transferability of Cu Pseudopotentials in Cu/CuO Systems, Density Funct. Theory Pract. Course. (2019). https://sites.psu.edu/dftap/2019/03/31/transferability-of-cu-pseudopotentials-in-cucuo-systems/.
[3] D.J. Singh, L. Nordstrom, Planewaves, Pseudopotentials and the LAPW Method, Second Edi, Springer US, 2006. doi:10.1007/978-0-387-29684-5.
[4] P.E. Blöchl, Projector augmented-wave method, Phys. Rev. B. 50 (1994) 17953–17979. doi:10.1103/PhysRevB.50.17953.
[5] X. Gonze, F. Finocchi, Pseudopotentials Plane Waves-Projector Augmented Waves: A Primer, Phys. Scr. T109 (2004) 40. doi:10.1238/physica.topical.109a00040.
[6] R.M. Martin, Electronic Structure: Basic Theory and Practical Methods, Cambridge University Press, 2004. doi:10.1017/CBO9780511805769.
[7] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77 (1996) 3865–3868. doi:10.1103/PhysRevLett.77.3865.
[8] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A.P. Seitsonen, A. Smogunov, P. Umari, R.M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials., J. Phys. Condens. Matter. 21 (2009) 395502. doi:10.1088/0953-8984/21/39/395502.
[9] A. Dal Corso, Pseudopotentials periodic table: From H to Pu, Comput. Mater. Sci. 95 (2014) 337–350. doi:10.1016/j.commatsci.2014.07.043.
[10] http://web.mit.edu/espresso_v6.1/amd64_ubuntu1404/qe-6.1/Doc/INPUT_LD1.html
[11] N. Troullier, J.L. Martins, Efficient pseudopotentials for plane-wave calculations, Phys. Rev. B. 43 (1991) 1993–2006. doi:10.1103/PhysRevB.43.1993.

 

Transferability of Cu Pseudopotentials in Cu/CuO Systems

1 Reply

Point of Study: 

The purpose of this work is to investigate the transferability of pseudopotentials in copper systems. When periodic atomic systems are modeled in Density Functional Theory (DFT), plane waves are often used as a basis set for wave functions. Though convenient for the periodicity of the system, large plane wave basis sets are needed to model the core of atoms where the wave function oscillates considerably.

Figure 1 ) To demonstrate the oscillations in the all electron wave function of an atom, a model radial wave function is plotted as the solid blue line. In the region closer to the nucleus (denoted as the core region in red), the wave function oscillates significantly over a short spatial extent in this region.

Approximations are often made to smooth out this portion of the wave function to reduce the need for massive PW basis sets. The pseudopotential approximation does just this; being constructed such that they reproduce the all electron wave function beyond a designated cutoff radius rc. (See figure 2) Pseudopotentials are defined by the potential of the core plus the screening from the core electrons. Though there are a number of formalisms to construct pseudopotentials,1–4 not all of these pseudopotentials are able to accurately reproduce the all electron wave-function in all environments. Most chemistry happens with the valence electrons, comprised mostly from the outer portion of the radial wave function, but this varies from system to system. The amount of electrons participating in bonding is system dependent,4 and pseudopotentials that work well for one system may not for one with different chemistry. This work investigates this ‘transferability’ in a Cu and Cu/O system.

Background:

Consider the following figure illustrating the comparison between pseudo-wave functions and the true all electron wave function reproduced from Methodology of Quantum Mechanics/Atomic Simulations.4

Figure 2 (left) Is a plot of the pseudo and all-electron wave functions. The pseudo wave function as the dashed line, and the all-electron is solid. Notice that the pseudo-wave function matches the all electron wave function at values at or above the core radius rc. (right) The screened pseudopotential is plotted (dashed) in comparison with the true unscreened, all-electron potential. Here the core radius is marked at the outermost node. The rc used for pseudization is not necessarily the true demarcation for core electrons, and its position can vary per pseudopotential formalism.

Pseudo wave functions (figure 2r) are constructed to obey the above principle constraint and others that depend on the method of generation (e.g. norm conservation). Simple pseudopotentials are constructed by inversion of the radial Kohn Sham equation after the pseudo wave function has been defined, with the screened all-electron potential. (See figure 3 below after introduction of the key equations) Usually this is done to match the single-atom all electron wave function in a certain electronic configuration from DFT. The radial Kohn-Sham equation is given by,

(1)

where, l is the angular quantum number, V is the Kohn-Sham potential, epsilon is the eigen energy, and R (r) is the radial wave function. Note that this is a function of principal quantum number as well as the angular quantum number but is left out for simplicity. Good pseudopotentials are expected to yield the same eigenvalues as the all-electron results, among other requirements the valence charge, and smoothness conditions outlined in reference 2. These vary per formalism, but it is typically required that:

(2)

The radial portion of the pseudo-wave function, R PP should match the all electron wave function R all e- outside the cutoff radius rc.

Figure 3 outlines the key steps in the generation of a pseudopotential. Note that this can be done with different electronic configurations. This is done for each angular momentum channel. See Kleinman and Bylander for further information on how this is done in practice.6 In the step of determining the cutoff radius, the core and valence electrons are partitioned; this should be done such that there is minimal core charge leakage. The Pseudo-Wave function yields a certain screened potential. The screening contributions from the Hartree and Exchange-Correlation potentials are removed to obtain the pseudopotential.

Note that the radial Kohn-Sham (eqn. 1) is a second order differential equation. When the energy and the potential are constant, the solution at a given position ro is uniquely determined by the wave function and its derivative. Or even more simply by the logarithmic derivative since it contains both:

(3)

The effect of evaluating this for different energies gives different (fixed) Kohn-Sham potentials. It offers a crude (but very fast) test of the transferability. If the Pseudo and all-electron wave function produce the same logarithmic derivative for more energies (and potentials implicitly) the pseudopotential used is more transferable. This is succinctly written:

(4)

This should be true when evaluated at ro greater than or equal to rc. High-quality pseudopotentials should reproduce the all-electron logarithmic derivative evaluated at energies near an eigenvalue of eqn 1. In this work, different Cu pseudopotentials commonly used in Quantum Espresso are tested with equation 4, with more rigorous and applied testing to follow in future posts.

Computational Methods

DFT calculations were performed using the Quantum Espresso 6.2 code. Though typically a plane wave code, the single atom calculations are done on a radial mesh using the ld1 module within Quantum Espresso. The single atom pseudopotential and all electron results were obtained on the default radial mesh designated in the pseudopotential files. The PSLib 1.0 PAW pseudopotential was used as the initial test pseudopotential due to its expected high transferability.7 This uses Cu in the configuration [Ar] 4s1.5   3d9.5 4p0.0 with minimum soft cutoff radii (l– dependent from local potential constructed from Bessel functions) of bohr, and a core cutoff radius of 1.0 bohr. This pseudopotential implements a scalar relativistic correction developed by Koelling and Harmon.8

Transferability is often a point of discussion during the generation of pseudopotentials; atoms do not stay in the same electronic configurations in all chemical environments. The pseudopotential is tested in different electronic configurations.

Results

Figure 3 (Top) Ground State Configuration: The pseudopotential reproduces the logarithmic derivative of the wave function well over a large energy range. Noting the energy scale is in Ry (~13.6 eV). The slight dip at the end of the plot corresponds to a node in the wave function. This becomes more evident when using another, slightly larger ro. Figure 3 (Bottom) Ground State at larger: The discontinuity in the derivative is due to a node in the pseudo and all electron wave functions.

In the results we see that the PAW pseudopotential provides an adequate description of the all electron wave function over a significant energy range. It is important to note that a deficiency of using logarithmic derivatives to test transferability is that the Kohn-Sham potential is considered at a fixed configuration. This changes as the chemical environment changes. In order to account for this to an extent, other electronic configurations are tested as well. (Including the configuration the pseudopotential was generated at) This was done for a range of ro values larger than the soft cutoff radius, and the results are appended. The overall results show that the logarithmic derivatives agree fairly well when varying the electronic configuration as well. This is important for Cu considering that it tends to lose electrons in many chemical environments.

Conclusions

The PAW pseudopotential tested is transferrable over a considerable range of energies and selected evaluation radii. The results would likely be that the pseudopotential would still reproduce reasonable energies in a chemically varied environment. (E.g. in CuO). Rigorous testing can be done by comparing the system energies produced by all-electron DFT simulations in different crystals to those of the pseudized cores. Due to the expense of all electron calculations for extended solids, a comparison will not be made with an all electron configuration. Instead, a comparison will be made for another copper pseudopotential with worse transferability seen in logarithmic derivatives. The two pseudopotentials will be analyzed in bulk Cu, CuO and compared to other studies in the literature.

  1. Hamann, D. R., Schlüter, M. & Chiang, C. Norm-Conserving Pseudopotentials. Phys. Rev. Lett. 43, 1494–1497 (1979).
  2. Troullier, N. & Martins, J. L. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B 43, 1993–2006 (1991).
  3. Kresse, G. & Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 59, 1758–1775 (1999).
  4. Umeno, Y., Shimada, T., Kinoshita, Y. & Kitamura, T. Methodology of Quantum Mechanics/Atomic Simulations. in Multiphysics in Nanostructures (eds. Umeno, Y., Shimada, T., Kinoshita, Y. & Kitamura, T.) 5–34 (Springer Japan, 2017). doi:10.1007/978-4-431-56573-4_2
  5. Dal Corso, A. Pseudopotentials periodic table: From H to Pu. Comput. Mater. Sci. 95, 337–350 (2014).
  6. Koelling, D. D. & Harmon, B. N. A technique for relativistic spin-polarised calculations. J. Phys. C Solid State Phys. 10, 3107–3114 (1977).

Appended Results:

Lattice Constants and Stability of FCC, SC, and HCP Platinum

Leave a reply

James Goff

Computational Methods

DFT simulations were performed using the CASTEP package implemented through MaterialsStudio70. The GGA was used for the density functional as implemented in Perdew-Burke-Ernzerhof.1 Ultra Soft Pseudopotentials (Pt_00.usp) were used to represent the atomic cores. These pseudopotentials have 10 valence electrons (5d96s1) in the configuration 1s22s22p63s23p63d104s24p64d105s25p64f145d96s1 with a cutoff radius of 2.4 au. A Monkhorst-pack grid of 10x10x10 was used to sample the Brillouin zone for primitive cubic cells.2 To maintain a consistent sampling of the Brillouin zone, the nearly equivalent kpoint density was achieved in HCP cells with a grid of 10x10x5. The electronic states were smoothed with 0.1 eV smearing. The convergence threshold for all SCF calculations was 10^-8 eV. It was found that a kinetic energy cutoff of 500 eV was required to converge energies to within 10 meV/atom. These results can be found in the appended convergence data. This cutoff was also selected such that the relative error in the pressure was < 0.1%. Cells were initially constructed such that the initial distance between Pt atoms was approximately 2.8 angstroms. To achieve this in HCP cells, the initial c/a ratio was approximately 1.75.

To obtain the lattice parameters for Platinum in FCC, SC, and HCP systems, the third-order Birch-Murnighan (BM) equation of state,

was used where Eo, Vo, and Bo are the system energy, system volume, and system bulk modulus at zero pressure, respectively. The coefficient, B’o , is the pressure derivative of the bulk modulus at constant temperature. These were the fitting parameters for the DFT data. This is outlined in Sholl et. al.3 The volume form of the equation is used so that Hexagonal systems can be fit as well. To fit the Equation of State (EOS), 5-8 SCF simulations were performed to obtain the system energy with varied lattice parameters. The range was constructed such that there were at least three data points to fit on either side energy vs volume minimum. The fits were obtained using Microsoft Excel’s least squares solver, and the BM parameters were optimized until the RMS error on the predicted energy was on the order of 10^-3 eV.

Results and Discussion

The minimum energies predicted with the BM equation for FCC, SC, and HCP Pt are 718.162, 717.612, and  718.104 eV, respectively. It is evident that FCC platinum is more stable than SC and HCP platinum. From figure 2b, the minimum energy for HCP Pt is slightly higher than that for FCC Pt by 0.06 eV. Experiment suggests that FCC Pt is the most stable under standard conditions with a lattice parameter of 2.772 Angstroms.4 The DFT results and BM fits for the FCC and SC lattices are shown in figure 1a and 1b Figure 2 outlines the DFT results and BM fits to HCP Pt.  The result obtained for the lattice constant of FCC Pt, at 2.828 Angstroms  is larger by 2 % than the experimental value. This error is likely due to ‘under-binding’ in solids common in PBE-GGA functionals. FCC is the most favorable lattice for Pt, the energy difference of 0.06 eV is well within the error present from the DFT simulations (based on 10 meV convergence criteria) and the RMS errors from the fits.

 

The initial c/a ratio of HCP Pt was tested with values of 1.65, 1.75, and 1.85, with the ratio 1.75 corresponding to distances between Pt of 2.82 angstroms. The BM equation was fit to the c/a ratio yielding the lowest energy; as discussed in 2b, this is the c/a ratio that yields Pt-Pt spacing most similar to that in the other Close-Packed structure (FCC).

 

 

Conclusions

The predicted structure for Pt is FCC with a primitive lattice constant of 2.828 Angstroms. While both HCP and FCC are significantly lower in energy than the SC lattice, the energy difference between the HCP and FCC structures is small. Only three c/a ratios were tested for HCP structures, and future work should test more c/a ratios near the ratio currently yielding the minimum. Though these lattices may be close in energy due to the previously mentioned similarities between FCC and HCP lattices.

References

  1. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865–3868 (1996).
  2. Monkhorst, H. J. & Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 13, 5188–5192 (1976).
  3. Sholl, D. & Steckel, J. A. Density Functional Theory: A Practical Introduction. (John Wiley & Sons, 2011).
  4. Callister, William D., and David G. Rethwisch. Materials science and engineering: an introduction. Vol. 7. New York: John wiley & sons, 2007.

Appended Convergence Data: