Accuracy of Partial Core Corrections using Fourier Transforms in

51 mins ago - Partial core corrections can be important in obtaining an accurate description of nonlinear exchange-correlation functionals and improvi...
1 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Accuracy of Partial Core Corrections Using Fourier Transforms in Pseudopotential−Density Functional Theory Weiwei Gao*,† and James R. Chelikowsky*,†,‡,§ †

Center for Computational Materials, Institute for Computational Engineering and Sciences, ‡Department of Physics, and Department of Chemical Engineering, The University of Texas at Austin, Austin, Texas 78712, United States

J. Chem. Theory Comput. Downloaded from pubs.acs.org by KAOHSIUNG MEDICAL UNIV on 11/22/18. For personal use only.

§

ABSTRACT: Partial core corrections can be important in obtaining an accurate description of nonlinear exchange-correlation functionals and improving the transferability of pseudopotentials. We show that a widely used procedure, which calculates partial core charge density, ρpartial core , in Fourier space and then converts it to real space with fast Fourier transforms, can lead to sizable numerical errors of exchange-correlation potentials in the vacuum region. Such errors occur in modeling lowdimensional materials or surfaces with supercells. The loss of accuracy originates from the slow-decaying feature of core charge density in reciprocal space. Numerical errors on the order of 1 eV in the Kohn−Sham energies of unoccupied states can occur in pseudopotential−density functional calculations. The direct calculation of the partial core charge in real space can avoid the numerical errors caused by Fourier transforms.

1. INTRODUCTION Pseudopotentials constructed within density functional theory (DFT) are popular theoretical tools for understanding and predicting material properties from first principles. In allelectron density function theory, the exchange-correlation energy Exc and potential Vxc are calculated as functionals of the ground-state electron charge density ρ(r) = ρcore(r) + ρval(r), where the total electron charge density ρ(r) is a summation of core charge ρcore and valence charge ρval. To avoid the computational difficulties caused by the large variations of ρcore and wave functions in the vicinity of nuclei, the pseudopotential theory was developed. Pseudopotentials include effects of both ions and core electrons in soft (i.e., smooth and slowly varying) potentials, which are optimized to reproduce allelectron calculations. In pseudopotential-DFT calculations, the exchange-correlation energy Exc and potential Vxc are treated as functionals of a smooth pseudovalence charge density ρps val ps given by ρval (r) = ∑ioccupied |ψips (r)|2 , where ψips (r) are pseudowave functions. This approximation is based on two assumptions: (1) The core charge density ρcore shows negligible changes under different chemical environments.1 (2) The spatial overlaps between core states and valence states are negligibly small such that Vxc[ρcore + ρval] ≅ Vxc[ρcore] + Vxc[ρval]. Assumption (1) is true for most materials, including gas, solids, nanocrystals, and molecules in a wide range of temperature and pressure, where core electrons are not strongly affected. However, assumption (2) need not hold in many systems, especially for spin-polarized systems and materials with transition metal elements or alkaline metal elements.2,3 In these systems, the contribution from the core states to exchange-correlation potential no longer remains constant as the chemical environment changes, since both Vxc and Exc are nonlinear functionals of ρ and there are sizable overlaps between ρcore and ρval. © XXXX American Chemical Society

The partial nonlinear core correction (NLCC) is an widely used method to account for the nonlinearity of exchangecorrelation potential and energy.2−4 With this method, one adds a partial core charge density ρpartial core to the pseudovalence charge density ρps in calculations of exchange-correlation val potential and energy partial ps Vxc(r ) = Vxc[ρcore + ρval ]

(1)

partial ps Exc = Exc[ρcore + ρval ]

(2)

ρpartial core (r)

is set to be the same as the all-electron core charge when the distance from the nucleus r is larger than a cutoff parameter rcore and a smooth function when r ≤ rcore. rcore is a parameter which needs to be carefully checked to reproduce the all-electron results.5 The advantages of using partial core corrections are twofold. First, partial core corrections are softer and require less computational cost than full core charge. Second, core corrections provide a faithful description of the nonlinearity of exchange-correlation effects and greatly improve the transferability of pseudopotentials. For periodic systems or confined systems modeled with periodic supercells, ionic potentials or core charges are commonly evaluated in reciprocal space and then Fourier transformed back to real-space grids.6 In general, given the distribution of atoms and the contribution to a quantity F from each atom, one wants to calculate quantity F(r) in a periodic system Natom

F(r) =

∑ ∑ F j(r − R − τj) j

R

(3)

Received: August 6, 2018 Published: November 6, 2018 A

DOI: 10.1021/acs.jctc.8b00820 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation where j ranges over all atoms in a unit cell, Fj is the contribution to the quantity F from atom j, and R runs over all periodic lattice vectors. F could be the local parts of pseudopotentials, partial core charge densities, or any other quantities of interest. The Fourier transform of F(r) is written as F(G) =

1 Ωcell

∫Ω

F(r)eiG·r d3r =

∑ F j(G)eiG·τ

j

cell

(4)

i

where Fj(G) = ∫ all space Fj(r)eiG·r d3r/Ωcell is the form factor of atom j and /Ωcell is the volume of one unit cell. With F(G) calculated for all G vectors satisfying |G| ≤ 2Ecut where Ecut is the kinetic energy cutoff, we perform a fast Fourier transform to obtain F(r) on a real-space grid. This procedure is usually much more efficient than the direct evaluation of eq 3 for the local part of pseudopotentials, because they decay as 1/r in real space and the summation over R converges slowly. However, we find that using this method (i.e., calculation in reciprocal space and then performing Fourier transform) to calculate ρpartial core (r) requires a more cautious treatment because of high-G Fourier components (i.e., ρpartial core (G) for large |G|) of the partial charge does not decay quickly in reciprocal space. The slow decay can result in numerical errors of the partial core charge and a finite value of the exchange-correlation potential in vacuum regions where both the charge and potentials are expected to be zero. In the following, we discuss the origin and repercussions of this numerical error and discuss possible solutions to the problem.

Figure 1. Comparison between the modulus of Fourier transform of partial core charge density |ρpartial core (G)| and the pseudocharge density |ρps val(G)| of a molybdenum atom. |ρ(G)| is inversely proportional to Ωcell. Here we use a 15.75 × 15.75 × 15.75 au3 supercell to calculate |ρ(G)|.

In comparison, the pseudovalence charge ρps val(G) is less than 10−12 au−3 when |G|2/2 > 300 Ry, which explains why the pseudocharge density converges faster with respect to planewave cutoff than partial core charge density. For partial core charge density, a typical plane-wave energy cutoff Ecut of 600 to −5 800 Ry leads to numerical error of δρpartial to 10−7 au−3 core = 10 partial if we use the Fourier transform to calculate ρcore in real-space grids. Such high precision is sufficient for calculations of bulk materials, where the pseudovalence charge density is significantly larger than δρpartial everywhere. However, for core low-dimension systems like 2D materials and nanocrystals, which are usually modeled in supercells with a large vacuum region, this precision can cause problems in calculating Vxc(r) where the valence charge is absent. To demonstrate how a small numerical error of charge density ρ can lead to a notable error in exchange-correlation potential Vxc, we show the LDA exchange-correlation potential Vxc as a function of charge density ρ in Figure 2a. As ρ varies from 10−8 to 10−5 au−3, which is around the typical precision of calculating core charge density with Fourier transform, the exchange-correlation potential Vxc changes from 0.1 to 1.0 eV. For example, in Figure 2c, we plot the in-plane average of partial core charge density ⟨ρpartial core ⟩ of MoSe2 along the zdirection, i.e., the direction perpendicular to the plane of the MoSe2 monolayer. If we use a plane-wave energy cutoff of 480 Ry for charge density, which is large enough to converge the DFT calculation without NLCC, we obtain a finite core charge density of approximately 10−5 au−3 in the vacuum region between MoSe2 slabs. Even with a high energy cutoff of 1000 Ry for charge density, the precision of calculated finite core charge density is still in the order of 10−6 au−3. As a result, the numerical error of calculating ρcore leads to a finite exchangecorrelation potential up to 0.8 eV in the interstitial region between layers, as shown in Figure 2d. In principle, Vxc should approach to zero as the position moves away from nuclei and charge density goes to zero. An error in Vxc directly affects the calculations of band structures. For example, we plot in Figure 3a the band structures of MoSe2 calculated with and without partial NLCC. The partial core charge densities ρcore(r) are evaluated by using two different methods. One method is the direct calculation of eq 3, which as expected gives a result agreeing well with that calculated without core correction, as shown by the green dots and black solid curves in Figure 3a. The other method is using

2. COMPUTATION DETAILS We use the plane-wave based code Quantum Espresso7,8 to calculate the results presented here and confirm our conclusions with another plane-wave code PARATEC9 and the real-space based code PARSEC.10−12 All our calculations were carried out in the local density approximation (LDA) of the exchange-correlation functional using the Perdew−Zunger parametrization.13 We use norm-conserving pseudopotentials generated with the Troullier−Martin’s method.14 Our conclusions also hold for generalized gradient approximation (GGA)15 and other types of pseudopotentials (as long as partial core corrections are included). In literature, rcore is usually chosen such that ρcore(rcore) = λρval(rcore), where λ ∈ [1,2]. Unless otherwise noted, we choose the cutoff radius rcore where ρcore(rcore) = 2ρval(rcore) (i.e., λ = 2) in our calculations. The parameters used for generating pseudopotentials are presented in Table 1. 3. RESULTS AND DISCUSSIONS First we show that ρpartial core (G) decays slowly in the reciprocal space. As shown in Figure 1, the modulus of ρpartial core (G) of Mo atom is on the order of 10−8 au−3 even for |G|2/2 ≈ 1200 Ry. Table 1. Matching Radii rcut for Orbitals and Parameters rcore for Partial Core Corrections Used for Generating Pseudopotentials element Mo Se Cr Cu

rcut for orbitals (au) 4d: 2.45; 5s: 4s: 1.90; 4p: 3d: 1.20; 4s: 3d: 1.20; 4s:

2.45; 1.90; 2.25; 2.20;

5p: 4d: 4p: 4p:

2.65 1.90 2.55 2.30

rcore (au) 1.44 1.15 1.20 0.37 B

DOI: 10.1021/acs.jctc.8b00820 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 2. (a) Exchange-correlation potential Vxc as a function of local charge density ρ. (b) Crystal structure of monolayer MoSe2 studied with a supercell with vacuum region separating the periodic images of monolayers. (c) In-plane average of partial core charge density of MoSe2, which is 1 partial partial ⟩ = A ∫ ρcore (x , y , z)| dx dy . A is the area of the unit cell. (d) In-plane average of exchange-correlation potential of MoSe2, defined as ⟨ρcore A

which is defined as ⟨Vxc⟩ =

1 A

∫A Vxc(x , y , z) dx dy . approximately as linear combinations of atomic orbitals localized around nuclei. ψ LUMO,Γ does not have any components in the interstitial region between MoSe2 layers. As a result, the energy of ψLUMO,Γ is not significantly affected by the errors of ρcore and Vxc. As another example, we study the work function of the copper [100] surface. We show the local part of total selfconsistent potential near the copper [100] surface in Figure 4.

Figure 3. (a) Comparison between the band structures of MoSe2 calculated with and without partial core correction. The partial core charge density is calculated comparatively in two different ways, namely, direct evaluation of ρcore in real space by using eq 3 and the Fourier transform method. The three right panels respectively show the in-plane average of modular square of (b) LUMO+8 at Γ point, (c) LUMO+4 at Γ point, and (d) LUMO at Γ point.

the Fourier transform to calculate partial core charges. This method produces a band structure which agrees well with the no-core-correction case in terms of occupied states and the first four conduction bands but shows significant differences for higher conduction bands. These affected bands are planewave-like states with large components in the vacuum region. We find using the Fourier transform to calculate ρpartial core results in errors as large as 0.65 eV in band structures of MoSe2. In the right panels of Figure 3, we plot the in-plane averaged modular square of wave functions for three representative unoccupied states. For the convenience of discussion, we label the ith state above the lowest unoccupied orbital at Γ point as ψLUMO+i,Γ. As clearly shown in Figure 3b,c, both ψLUMO+8,Γ and ψLUMO+4,Γ have substantial components located in the interstitial vacuum region between two slabs. They are more heavily affected by the error of Vxc in the vacuum region. This is reflected in the band structure in Figure 3, as we discussed previously. In comparison, ψLUMO,Γ can be considered

Figure 4. Average local part of the self-consistent potential ⟨Vxc + VHartree + Vion⟩ near the copper [100] surface. We compare the potential curves calculated with core correction using three different rcore parameters 0.80, 0.45, and 0.37 au, which correspond to conditions ρcore(rcore) = λρval(rcore) with λ = 1.0, 1.5, and 2.0, respectively.

The difference between the potential at the vacuum region and Fermi energy EFermi gives us the work function. Without including core corrections, our calculation converges with a plane-wave cutoff 150 Ry for wave functions and 600 Ry for charge density, and the calculated work function of the copper [100] surface is 4.8 eV, agreeing well with previous work.16 For the potential curves calculated with partial core corrections, as shown by the color lines in Figure 4, we use an enormous 2000 Ry cutoff for charge density and 500 Ry for wave functions. Unfortunately, such a large cutoff still cannot converge ρpartial core C

DOI: 10.1021/acs.jctc.8b00820 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation Table 2. Total Energy Etot and Forces between Two Cr Atoms Separated by Distance da Etot (Ry)

force (Ry/au)

d (au)

real-space

Fourier trans.

diff.

real-space

Fourier trans.

2.8 3.2 3.6 4.0

−58.04461 −57.96008 −57.82294 −57.68989

−58.04236 −57.95815 −57.82079 −57.68802

0.002 0.002 0.002 0.002

0.00518 −0.32287 −0.34759 −0.31334

0.00515 −0.32287 −0.34771 −0.31335

diff. −3 −6 1 −4

× × × ×

10−5 10−6 10−4 10−6

a We compare the results calculated with partial core corrections calculated with the real-space method (i.e., direct evaluation of eq 3) and Fourier transform.

and exchange correlation potential in vacuum. By comparing potential curves calculated with different rcore for partial core corrections and that without core correction, one can clearly see the numerical errors in Vxc become larger as rcore decreases and partial core charge becomes more localized. With rcore = 0.37 au, the calculated Fermi energy is affected by a few millielectronvolts, but the total local potential at vacuum is underestimated by 1.7 eV, which results in an underestimated work function of 3.1 eV. In comparison, if one uses rcore = 0.80 au and a large 2000 Ry cutoff for charge density, the calculated work function is 4.6 eV, which still has an error of −0.2 eV. Another common method is using the electrostatic potential (i.e., Hartree potential) instead of the total local self-consistent potential in vacuum to calculate work function.17 This can avoid the numerical errors of exchange-correlation potential caused by core corrections. However, the error of Vxc can still affect the plane-wave-like states with large components in vacuum. Our benchmarks demonstrate that this problem affects the total energy by less than 10 meV/atom and the forces less than 10−4 Ry/au, and both are negligible. For example, we calculate the total energy and forces of two chromium atoms separated by a distance d. The calculations are performed with a 31.5 × 31.5 × 31.5 au3 large supercell, and a plane-wave cutoff of 560 Ry for charge density. We compare the results obtained with partial core corrections calculated directly in real space and with Fourier transform, as shown in Table 2. The differences between the total energies calculated with two methods are about 0.002 Ry, and the differences of force range from 10−6 to 10−4 Ry/au. We stress that the problem caused by using Fourier transform and an insufficient plane-wave cutoff only appears in the calculations of systems that contain large vacuum regions, such as low-dimension materials and surfaces of bulk materials, and affects the energy of unoccupied states with large components in the vacuum region. The propagation of errors could also affect the calculations which involve optical transitions since the relative positions between occupied and unoccupied states are affected. However, the magnitudes of errors vary over a large range depending on systems. In Figure 5, we plot the partial core charge in reciprocal space ρpartial core (G) for a series of elements. As clearly shown in Figure 5, this issue is severe for first-row elements such as oxygen and nitrogen and third-row transition metals such as nickel and cobalt, −6 which have large ρpartial au−3 at highcore (G) in the order of 10 energy range. One may argue that it is unnecessary to use NLCC for first-row elements because they only have the highly localized 1s core orbital. This is true for many applications,2,4 but including NLCC for first-row elements is important for describing chemical reactions, radicals, and surfaces which have spin-polarized states.5 For alkaline metal elements, the problem is minor because rcore for these elements are relatively large and

Figure 5. Modulus of partial core charge density |ρpartial core (G)| (in reciprocal space) for several elements. −9 ρpartial au−3 at around 800 Ry, as shown in core (G) decreases to 10 the bottom panel of Figure 5. The numerical errors caused by Fourier transforms can be partially alleviated by using a smoother function for partial core charge inside rcore to suppress the high-G Fourier components. The original functional form3 of partial core correction is

l A sin(Br ) o o o , if r ≤ rcore o o partial r ρcore (r) = m o o o o o ρcore (r ) , if r > rcore (5) n where the core charge inside the sphere of radius rcore is replaced by a spherical Bessel function. Better functional n 2i 5 forms,2,5,18 such as ρpartial core (r < rcore) = exp (∑i cir ), have been proposed in order to guarantee the continuity of higher derivatives of ρpartial core at r = rcore. With these smoother functional forms, we can effectively reduce the plane-wave energy cutoff to converge calculations of core charge density in real space. Using Mo as an example, we plot its partial core charge |ρpartial core (G)| with different functional forms in Figure 6. It is 2 clear that using the function ρpartial core (r < rcore) ≡ exp(c0 + c1r + partial 4 c2r ) effectively reduces the magnitude of ρcore (G) for large |G| as compared to using the spherical Bessel function J0. However, even with smooth functions such as exp(c0 + c1r2 + c2r4), a large energy cutoff is still required to converge the D

DOI: 10.1021/acs.jctc.8b00820 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation Funding

We acknowledge support from a subaward from the Center for Computational Study of Excited-State Phenomena in Energy Materials at the Lawrence Berkeley National Laboratory, which is funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division under Contract No. DEAC02-05CH11231, as part of the Computational Materials Sciences Program. Computational resources are provided in part by the National Energy Research Scientific Computing Center (NERSC) and the Texas Advanced Computing Center (TACC). Figure 6. Comparison between different partial core charge density ρpartial core (G) of Molybdenum obtained with two different functions.

Notes

The authors declare no competing financial interest.

■

ρpartial core (r)

(1) Goedecker, S.; Maschke, K. Transferability of pseudopotentials. Phys. Rev. A: At., Mol., Opt. Phys. 1992, 45, 88−93. (2) Willand, A.; Kvashnin, Y. O.; Genovese, L.; Vázquez-Mayagoitia; Deb, A. K.; Sadeghi, A.; Deutsch, T.; Goedecker, S. Norm-conserving pseudopotentials with chemical accuracy compared to all-electron calculations. J. Chem. Phys. 2013, 138, 104109. (3) Louie, S. G.; Froyen, S.; Cohen, M. L. Nonlinear ionic pseudopotentials in spin-density-functional calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1982, 26, 1738−1742. (4) van Setten, M.; Giantomassi, M.; Bousquet, E.; Verstraete, M.; Hamann, D.; Gonze, X.; Rignanese, G.-M. The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table. Comput. Phys. Commun. 2018, 226, 39−54. (5) Porezag, D.; Pederson, M. R.; Liu, A. Y. Importance of nonlinear core corrections for density-functional based pseudopotential calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 60, 14132−14139. (6) Martin, R. M. Electronic structure basic theory and practical methods; Cambridge University Press: 2008; Chapter 12. (7) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Dal Corso, A.; de Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; MartinSamos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009, 21, 395502. (8) Giannozzi, P.; Andreussi, O.; Brumme, T.; Bunau, O.; Buongiorno Nardelli, M.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Cococcioni, M.; Colonna, N.; Carnimeo, I.; Dal Corso, A.; de Gironcoli, S.; Delugas, P.; DiStasio, R. A.; Ferretti, A.; Floris, A.; Fratesi, G.; Fugallo, G.; Gebauer, R.; Gerstmann, U.; Giustino, F.; Gorni, T.; Jia, J.; Kawamura, M.; Ko, H.-Y.; Kokalj, A.; Kücu̧ ̈kbenli, E.; Lazzeri, M.; Marsili, M.; Marzari, N.; Mauri, F.; Nguyen, N. L.; Nguyen, H.-V.; Otero-de-la Roza, A.; Paulatto, L.; Poncé, S.; Rocca, D.; Sabatini, R.; Santra, B.; Schlipf, M.; Seitsonen, A. P.; Smogunov, A.; Timrov, I.; Thonhauser, T.; Umari, P.; Vast, N.; Wu, X.; Baroni, S. Advanced capabilities for materials modelling with QUANTUM ESPRESSO. J. Phys.: Condens. Matter 2017, 29, 465901. (9) http://www.nersc.gov/users/software/applications/materialsscience/paratec/ (accessed July 20, 2018). (10) Kronik, L.; Makmal, A.; Tiago, M. L.; Alemany, M. M. G.; Jain, M.; Huang, X.; Saad, Y.; Chelikowsky, J. R. PARSEC the pseudopotential algorithm for real-space electronic structure calculations: recent advances and novel applications to nano-structures. Phys. Status Solidi B 2006, 243, 1063−1079. (11) Chelikowsky, J. R.; Kronik, L.; Vasiliev, I.; Jain, M.; Saad, Y. Using real space pseudopotentials for the electronic structure problem. Handbook of Numerical Analysis 2003, 10, 613−637.

calculation of using Fourier transforms, especially for some of the first-row elements and 3d transition metals, which require small rcore. In fact, an accurate and straightforward solution is directly calculating ρcore(r) in real space by using eq 3. We have demonstrated this solution in the monolayer MoSe2 example. There are a few reasons for doing this. First and most important, we can eliminate the errors of ρcore and Vxc due to the double Fourier transform and plane-wave cutoff. This method works for all functional forms of ρpartial and does not core require testing convergence and generating new pseudopotentials. Second, core charge corrections only affect the local part of the self-consistent potential in DFT calculations; thus, it is justified to only evaluate them in real-space grids and save the memory of storing them in plane-wave basis. Third, since the core charge decays fast in real space, the summation over R in eq 3 converges quickly. Last but not least, based on our experiences, we believe implementing this approach only needs minor changes in plane-wave codes.

4. CONCLUSION As various new algorithms are developed and implemented in DFT computer codes, first-principles methods based on DFT continue evolving to gain unprecedented accuracy and efficiency, and the precision and reproducibility19 of DFT calculations become increasingly important. We find a commonly used approach for calculating partial core charge density in real-space grids using the Fourier transform can lead to notable numerical errors in exchange-correlation potential of up to 1 eV, which is sufficiently large to cause reproducibility problems in pseudopotential−DFT calculations. The problem stems from the slow-decaying feature of partial core charge in reciprocal space. We demonstrate this problem with monolayer MoSe2 and also show it is hard to resolve by increasing planewave cutoff with the example of the copper surface. We suggest using a smoother function for partial core charge within the sphere of cutoff radius rcore or to simply perform direct calculation of partial core charge ρpartial core in a real-space grid to avoid this problem.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Authors

*(W.G.) E-mail: [email protected]. *(J.R.C.) E-mail: [email protected]. ORCID

Weiwei Gao: 0000-0003-4477-393X E

DOI: 10.1021/acs.jctc.8b00820 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation (12) Zhou, Y.; Saad, Y.; Tiago, M. L.; Chelikowsky, J. R. Phys. Rev. E 2006, 74, 066704. (13) Perdew, J. P.; Zunger, A. Self-interaction correction to densityfunctional approximations for many-electron systems. Phys. Rev. B: Condens. Matter Mater. Phys. 1981, 23, 5048−5079. (14) Troullier, N.; Martins, J. L. Efficient pseudopotentials for planewave calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 1993−2006. (15) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (16) Patra, A.; Bates, J. E.; Sun, J.; Perdew, J. P. Properties of real metallic surfaces: Effects of density functional semilocality and van der Waals nonlocality. Proc. Natl. Acad. Sci. U. S. A. 2017, 114, E9188− E9196. (17) Singh-Miller, N. E.; Marzari, N. Surface energies, work functions, and surface relaxations of low-index metallic surfaces from first principles. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 235407. (18) Fuchs, M.; Scheffler, M. Ab initio pseudopotentials for electronic structure calculations of poly-atomic systems using density-functional theory. Comput. Phys. Commun. 1999, 119, 67−98. (19) Lejaeghere, K.; Bihlmayer, G.; Björkman, T.; Blaha, P.; Blügel, S.; Blum, V.; Caliste, D.; Castelli, I. E.; Clark, S. J.; Dal Corso, A.; de Gironcoli, S.; Deutsch, T.; Dewhurst, J. K.; di Marco, I.; Draxl, C.; Dułak, M.; Eriksson, O.; Flores-Livas, J. A.; Garrity, K. F.; Genovese, L.; Giannozzi, P.; Giantomassi, M.; Goedecker, S.; Gonze, X.; Grånäs, O.; Gross, E. K. U.; Gulans, A.; Gygi, F.; Hamann, D. R.; Hasnip, P. J.; Holzwarth, N. A. W.; Iuşan, D.; Jochym, D. B.; Jollet, F.; Jones, D.; Kresse, G.; Koepernik, K.; Kücu̧ ̈kbenli, E.; Kvashnin, Y. O.; Locht, I. L. M.; Lubeck, S.; Marsman, M.; Marzari, N.; Nitzsche, U.; Nordström, L.; Ozaki, T.; Paulatto, L.; Pickard, C. J.; Poelmans, W.; Probert, M. I. J.; Refson, K.; Richter, M.; Rignanese, G.-M.; Saha, S.; Scheffler, M.; Schlipf, M.; Schwarz, K.; Sharma, S.; Tavazza, F.; Thunström, P.; Tkatchenko, A.; Torrent, M.; Vanderbilt, D.; van Setten, M. J.; Van Speybroeck, V.; Wills, J. M.; Yates, J. R.; Zhang, G.X.; Cottenier, S. Reproducibility in density functional theory calculations of solids. Science 2016, 351, aad3000.

F

DOI: 10.1021/acs.jctc.8b00820 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX