Projector Augmented Wave Method Incorporated into Gauss-Type

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Projector Augmented Wave Method Incorporated into Gausstype Atomic Orbital Based Density Functional Theory Xiao-Gen Xiong, and Takeshi Yanai J. Chem. Theory Comput., Just Accepted Manuscript • Publication Date (Web): 22 May 2017 Downloaded from http://pubs.acs.org on May 22, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Projector Augmented Wave Method Incorporated into Gauss-type Atomic Orbital Based Density Functional Theory Xiao-Gen Xiong∗,† and Takeshi Yanai∗,‡ †Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China ‡Department of Theoretical and Computational Molecular Science, Institute for Molecular Science, Okazaki, 444-8585 Aichi, Japan E-mail: [email protected]; [email protected] Abstract The Projector Augmented Wave (PAW) method developed by Blöchl is well recognized as an efficient, accurate pseudopotential approach in solid-state density functional theory (DFT) calculations with plane-wave basis. Here we present an approach to incorporate the PAW method into the Gauss-type function (GTF) based DFT implementation, which is widely used for molecular quantum chemistry calculations. The nodal and high-exponent GTF components of valence molecular orbitals (MOs) are removed or pseudized by the ultrasoft PAW treatment, while there is elaborate transparency to construct accurate and well-controlled pseudopotential from all-electron atomic description and to reconstruct all-electron form of valence MOs from the pseudo MOs. This smoothess should benefit the efficiency of GTF-based DFT calculations in terms of elimination of high-exponent primitive GTFs and reduction of grid points in the numerical quadrature. The processes of the PAW method are divided into basis-independent

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and -dependent parts. The former is carried out using the previously-developed PAW libraries libpaw and atompaw. The present scheme is implemented by incorporating libpaw into the conventional GTF-based DFT solver. The details of the formulations and implementations of GTF-related PAW procedures are presented. The test calculations are shown for illustrating the performance. With near-complete GTF basis at the cc-pVQZ level, the total energies obtained using our PAW method with suited frozen core treatments converge to those with the conventional all-electron GTF-based method with a rather small absolute error.

1

Introduction

The Kohn-Sham (KS) approach to the density function theory (DFT) is the most widely used as a workhorse method of quantum chemistry (QC) calculations. 1 Its numerical implementation for practical molecular calculations was studied in late 80’s or early 90’s in the pioneering works of Becke, 2–4 Pople, 5–7 etc. The KS equation in a finite orbital basis expansion with atom-centered Gauss-type functions (GTFs) can be formulated as an analog of the Roothaan-Hall equation, whose implementation was established long before. 5 This schematic analogy enables the KS module to be easily incorporated into the existing GTFbased quantum chemistry implementation. Therefore, there are a number of QC program packages using the GTF basis for molecular DFT calculations, such as gaussian, 5 orca, 8 molpro, 9 nwchem, 10 etc. Beside the historical background, the use of GTFs for orbital basis facilitates the all-electron treatment of molecular orbitals (MOs) entailing core states, and the analytic evaluation of electron repulsion integrals 11 including the Hartree-Fock exchange, which is required for popular, accurate hybrid functionals, such as B3LYP, 4 etc. The cusp structures arising in MOs have steeply-varying shapes related to exponential form in the vicinity of nuclei, but can be compactly represented by contractions of several primitive GTFs with high exponents. The plane-wave functions are widely used as one-electron basis of first-principles solid2 ACS Paragon Plus Environment

Page 2 of 47

Page 3 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


state calculations based on KS-DFT. 12 They are considered to be numerically efficient because of the periodicity for treating crystalline solids and the connectivity to the fast Fourier transformation (FFT) technique for evaluating electronic Coulomb interactions. It is, however, well recognized that the atom-centered cusp shapes in MOs cannot be readily represented with the plane-wave expansion because of their high-frequency complexity. Therefore, most of the plane-wave-based implementations use the pseudopotentials, 12 except some mixed-basis approaches that additionally use local basis. 13–15 The pseudopotential approach is built upon the frozen core approximation that eliminates atomic core states, resulting in effective removal of the nodal and high-frequency structures of valence orbitals. The effective core potential (ECP) method is a widely-used pseudopotential approach in quantum chemistry mainly for simplifying the treatment of core electrons and associated relativist effects of transition metals and heavier elements. 16–18 Although there have been various pseudopotential schemes developed, the Projector Augmented Wave (PAW) method introduced by Blöchl 19 is considered to be remarkably successful in plane-wave-based material calculations. It shares some features of Vanderbilt’s soft pseudopotential formalism 20 and Blöchl’s earlier work on generalized separable potentials. 21 One of the conceptual advantages of the PAW method over the other pseudopotential schemes is a well controlled approximation connecting all-electron formalism to pseudopotentials and auxiliary smooth wave functions, referred to as pseudo wave functions. It has been built into several plane-wave implementations, such as vasp, 22 abinit, 23 quantum espresso, 24 pwpaw, 25 etc, and also in real-space methods. 26–28 Recently, Rangel et al. developed a portable PAW library, called libpaw, 29 which is originally the plane-wave PAW module hard-wired to abinit, developed by Torrent et al. 30 It was modularized for the sake of its porting to the wavelet-based DFT code, bigdft, 29,31 in such a way that the resulting isolated library provides implementations of core PAW procedures with basis-independent data interfaces and thus potentially can be combined with any other basis used in different programs. The PAW atomic pseudopotential data



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

required by the libpaw library can be generated using the atomic code atompaw, which was developed by Holzwarth et al. 32 These programs are freely accessible and indeed serve as key computational infrastructures in the present study. In this work, we present an approach to incorporate the PAW method into the GTFbased molecular DFT implementation. It is hereafter referred to as GTF-PAW method. Its basic ansatz is to use GTFs in place of plane-waves for the basis expansion of pseudo wave functions. Our computational method is developed by partly using the existing PAW modules provided by the atompaw and libpaw libraries in order to carry out basis-independent PAW procedures; thus, our GTF-related implementation is complementary. However, various additional types of intermediate integrals arise in the GTF-based PAW procedures, and should be treated in a numerically optimal and safe way. The earlier pseudopotential schemes with GTF basis used in ECPs of quantum chemistry codes 16,17 and in the Gaussian and plane waves method (GPW) and its augmented extension (GAPW) method of the cp2k code 14,15 are categorized to the norm-conserving pseudopotential, 12,33 while the PAW method inherits Vanderbilt’s ultrasoft scheme. 21,22 Along a direction similar to our ansatz, the PAW method using local atom-centered numerical basis with the uniform real-space grids was developed in the DFT code gpaw based on the multigrid finitedifference approach. 27 Recently, Booth et al. introduced a method to construct periodic pseudized Gaussian basis within the plane-wave PAW framework as efficient basis for manybody calculations for periodic systems. 34 To the best of our knowledge, there is no earlier report that has shown a combination of the PAW method and traditional GTF-based DFT implementations. The central process of DFT calculation in the QC implementation is the numerical quadrature to evaluate exchange-correlation (XC) energy and potential matrix. It is done in most of the GTF-based QC programs by using the Becke’s multi-center fuzzy cell quadrature method. 2 The optimization of the radial quadrature on top it was intensively studied in order to efficiently describe atom-centered steep cusps and slowly-decaying tails of MOs


Page 4 of 47

Page 5 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


on equal footing. 2,6,35–37 The Lebdev grids 38 with octahedral symmetry (Oh ) is generally used for spherical integration. The present GTF-PAW implementation retains the use of the Becke quadrature scheme for the XC evaluation. When we use it for evaluating smoothed pseudo wave functions, there arises the numerical advantage over all-electron GTF calculations, which is potentially twofold: (1) elimination of high-exponent primitive GTFs, and (2) reduction of grid points in Becke quadrature. The PAW method is considered to be more efficient, accurate, and transparent than the norm-conserving pseudopotentials which encompass ECP, so that our approach can be a better alternative to the ECP method most widely-used in QC to date. The paper is organized as follows: In Sec. 2, the general theory of the GTF-PAW method is presented. In Sec. 3, we detail the numerically optimal and safe way to implement the PAW method within the conventional GTF-based framework. Numerical validations are shown in Sec. 4, followed by the conclusions in Sec. 5

2

Basic theory

Prior to describing the formulations of theory, let us here show a summary of the indices, functions, and symbols frequently used throughout this work in Table 1.

2.1

PAW-based Kohn-Sham method

We begin by giving a brief overview of the PAW-based DFT. The reader may consult Refs. 19,22,29,30,32,39,40 for further details. At the heart of the PAW method is to obtain the ˜ pi true valence one-electron wave function |Ψp i by transforming the pseudo wave function |Ψ with a predetermined operator T , as follows: ( ˜ pi = |Ψp i = T |Ψ

1+

X ai

) ˜ pi . |φai i − |φãi i h˜ pai | |Ψ


(1)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 1: Indices, functions, and symbols p a i, j ni , li , mi , L, M µ, ν, λ, ξ Ψp (r) ˜ p (r) Ψ φai (r) φãi (r) pãi (r) χµ (r) n ˜ (r) n ˆ (r) n ˜ c (r) n ˜ Zc (r) v H [..](r) v xc [..](r) E H [..] E xc [..] (a) gL (r) SLM (θ, ϕ) Ra a ) Z a (Zeff a rpaw a rshape a rcore h

Valence orbital index Atom index PAW atomic orbital indices Atomic quantum numbers Gauss-type atomic orbital indices True one-electron wave function Pseudo one-electron wave function True atomic orbital Pseudo atomic orbital Projector Gauss-type atomic orbital basis Pseudo density function Compensation charge Pseudo core density function n ˜ c (r) augmented with nuclear charge Hartree potential (eqn. (32)) Exchange-correlation potential Hartree energy Exchange-correlation energy Shape function Real spherical harmonics Atomic center (Effective) nuclear charge Radius for augmentation region Radius for compensation charge Radius for core density Mesh size of atomic fine mesh


Page 6 of 47

Page 7 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


˜ p i is obtained to be a numerically smooth function. Note that |Ψp i is We assume that |Ψ often referred to as “all-electron” even if the core-electron wave functions are discarded due to the frozen core treatment built into the present scheme. The operator T is expressed using the true and pseudo atomic partial waves, denoted φai and φãi , respectively, and the projector functions, pãi , which all are atom-centered functions on the atom a, as given by: (a)

φai (r) = φi (r − Ra ) , (a) φãi (r) = φ˜i (r − Ra ) ,

(2)

(a)

pãi (r) = p˜i (r − Ra ) . (a) (a) (a) The atomic functions, φi , φ˜i , and p˜i , are given as separated forms into the radial and

angular functions: (a)

(a) φi (r)

φ (r) = ni li Sli mi (θ, ϕ) r

(3)

(a) φñi li (r) Sli mi (θ, ϕ) r

(4)

(a) φ˜i (r) =

(a)

(a) p˜i (r)

p˜ (r) = ni li Sli mi (θ, ϕ) r

(5)

where r = r ( sin θ cos ϕ, sin θ sin ϕ, cos θ ) with 0 ≤ θ ≤ π and 0 ≤ ϕ ≤ 2π, and the angular part Sli mi (θ, ϕ) is the real spherical harmonics. We only access these radial functions within the so-called augmentation region Ωa , which is a finite atom-centered sphere volume (a) (a) (a) a (r < rpaw ). This is related to the fact that pñi li (r) and φni li (r) − φñi li (r) vanish outside Ωa .

In the actual implementation, they are provided as tabulated data, which are determined by the preceding atomic calculations. ˜ p i in the PAW-based The working equation to determine the pseudo wave function |Ψ



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 47

Kohn-Sham (KS) method is written as:

˜ Ψ ˜ p i = εp S| ˜Ψ ˜ pi , H|

(6)

˜ and orthonormality operator S˜ are given by: where the effective one-electron Hamiltonian H ˜ = − 1 ∇2 + v˜eff (r) + H 2

S˜ = 1 +

X aij

X aij

a h˜ paj | , |˜ pai iDij

|˜ pai isaij h˜ paj | .

(7)

(8)

The term v˜eff (r) is the smooth (or pseudo) local potential, the form of which will be shown later. The atomic overlap matrix saij is defined by:

saij = hφai |φaj i − hφãi |φãj i ,

(9)

and is calculated outside the self-consistent field (SCF) loop. The definition of the nonlocal a coefficient matrix Dij is written as:

a ˜ a |φã i , Dij = hφai |H a |φaj i − hφãi |H j

(10)

˜ a are the one-center effective Hamiltonians based on all-electron and pseudo where H a and H a electron densities, respectively; thus, Dij are optimized in response to the densities updated a every SCF cycle. The further details of the nonlocal coefficients Dij can be seen in Refs.

19,22,30,40. Because of the presence of the nonlocal terms, the core MOs do not arise in eigenstates resulting from the solution of eqn. (6).


Page 9 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.2

Gauss-type atomic orbital expansion

Here let us introduce a scheme to incorporate the PAW scheme into the molecular DFT implementation based on Gauss-type basis functions (GTFs). The central ansatz proposed in this study is to use a linear combination of Gauss-type atomic orbitals (AOs) for representing ˜ p (r). It is formulated as: the pseudo one-electron wave function Ψ ˜ pi = |Ψ

X µ

|χµ iCµp

(11)

where |χµ i are atom-centered GTFs. 11 This basis expansion is the same as done in the conventional all-electron GTF-based implementation. The main task of our PAW-based calculation is to determine the molecular orbital (MO) coefficients {Cµp } by solving eqn. (6) in an SCF manner; they represent pseudo MOs as eigenfunctions. The matrix form of eqn. (6) in the AO basis is thus written as:

˜ p = εp SC ˜ p HC

(12)

˜ µν = hχµ |H|χ ˜ νi H

(13)

˜ µν = hχµ |S|χ ˜ νi S

(14)

with the p-th eigenvector, Cp = {Cµp }. Given a set of the GTF-based AOs {χµ }, we introduce a matrix form of the projectors, denoted {˜ paiµ }, as the projection of pãi (r) onto the AOs, defined by: pãiµ

=

h˜ pai |χµ i

Z = Ωa

(a)

dr p˜i (r − Ra ) χµ (r) .

(15)

The matrix elements pãiµ are evaluated prior to the SCF cycles and stored into memory. Using ˜ µν (eqn. (13)) and S ˜ µν (eqn. (14)) are written in algebraic form them, the nonlocal terms in H



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 47

as follows: ˜ µν = Fµν + H

X

a a a p˜iµ p˜jν , Dij

(16)

saij pãiµ pãjν ,

(17)

aij

˜ µν = Sµν + S

X aij

with the Fock-like and overlap matrices in AO basis given by:

Fµν = χµ − 12 ∇2 + v˜eff (r) χν ,

(18)

Sµν = hχµ |χν i ,

(19)

respectively. The overlap of the projector and pseudo MO is calculated by a matrix multiplication: ˜ pi = h˜ pai |Ψ

X

pãiµ Cµp

(20)

µ

This expression lets the PAW transformation (eqn. (1)) be proportionally parameterized by the MO coefficients as:

|Ψp i =

with |δµ i =

P ai

X µ

{|χµ i + |δµ i} Cµp ,

|φai i − |φãi i pãiµ , which can be viewed as a correction to AO basis.


(21)

Page 11 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.3

Density matrices and functions

As similarly done in the conventional implementation, the AO-basis density matrix fµν is constructed as:

fµν =

X

fp Cµp Cνp

(22)

p

with the occupancies fp , in which there is no contribution from core states. The pseudo density function is built upon it as given by:

n ˜ (r) =

X

fµν χµ (r)χν (r)

(23)

µν

The one-center density matrix introduced in the PAW method can be calculated using fµν as: ρaij =

X p

˜ p ihΨ ˜ p |˜ fp h˜ pai |Ψ paj i =

X

fµν pãiµ pãjν

(24)

µν

The one-center density functions associated with atomic all-electron and pseudo partial waves are parameterized with ρaij as: n1a (r) =

X

ρaij φai (r)φaj (r) ,

ij 1a

n ˜ (r) =

(25)

X

ρaij φãi (r)φãj (r) ,

ij

respectively. The so-called compensation charge is given by n ˆ a (r) =

X ij

ˆ (a) (r − Ra ) ρaij Q ij


(26)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 47

with ˆ (a) (r) = Q ij

X

(a)

qijaLM gL (r)SLM (θ, ϕ) ,

(27)

LM (a) where qijaLM correspond to the moments of (φai φaj − φãi φãj ) and gL (r) are the so-called shape (a)

functions, whose L-th moments are unity by construction. Note that gL (r) is zero for (a)

a . Some functional forms of gL (r) have been proposed by earlier works 19,22 and r > rshape

implemented in atompaw and libpaw. 29,30,32 Users are assumed to choose the type of the ˆ (a) are predefined objects in form when constructing atomic data; therefore, the functions Q ij molecular calculation. In the frozen core approximation of the PAW method, electrostatic interaction of valence states with core states are accounted for by incorporating predetermined atomic core densities into the KS equation, while the core descriptions are left unrelaxed. Atomic data provided by atompaw entail tabulated radial functions of the true and smooth (or pseudo) core densities (a)

(a)

for each atom, denoted nc (r) and n ˜ c (r), respectively, assuming that they are spheres in real space, as given by a nac (r) = n(a) c (|r − R |) ,

n ˜ ac (r) (a)

=

n ˜ (a) c (|r

(28)

a

− R |) .

(a)

The radial data of nc and n ˜ c fully cover a range up to the radius where they are numerically a a vanishing (rcore ); thus, it is much larger than the augmentation radius rpaw . The sum of the

atomic core density and point nuclear charge is denoted naZc (r) (= −Z a δ(r − Ra ) + nac (r)), and its smooth description n ˜ aZc (r). The superposition of n â, n ˜ ac , and n ˜ aZc for polyatomic system is written as: n ˆ=

X a

n â,

n ˜c =

X a

n ˜ ac ,

n ˜ Zc =

X a


n ˜ aZc ,

(29)

Page 13 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the atom-wise sparsity of which can be easily exploited in the molecular implementation.

2.4

Effective potentials

Given the aforementioned density matrices and density functions, the following three kinds of the effective potentials constitute the main parts of the PAW-based KS equation:

v˜eff (r) = v H [˜ n+n ˆ+n ˜ Zc ](r) + v xc [˜ n(r) + n ˜ c (r)]

1a veff (r) = v H [n1a + naZc ](r) + v xc [n1a (r) + nac (r)] , 1a v˜eff (r)

H

1a

a

= v [˜ n +n ˆ +

n ˜ aZc ](r)

xc

1a

+ v [˜ n (r) +

(30)

(31)

n ˜ ac (r)] ,

where v H means the Hartree potential defined by Z

H

v [n](r) =

dr0

n(r0 ) , |r − r0 |

(32)

and v xc refers to the exchange correlation potential. This formalism forgoes explicit inclusion of the compensation charges n ˆ and n ˆ a as a part of the density function for v xc , as done in the original paper and also carefully analyzed in the work of Torrent et al. 39 a arising in the These potentials are all used for evaluating the nonlocal coefficients Dij a effective Hamiltonian. The formula of Dij consists of three integral elements based on v˜eff , 1a 1a veff , and v˜eff , respectively, as follows:

a a 1a 1a ˆ ij ˜ ij + Dij −D , Dij =D

a ˆ ij D

Z

(33)

ˆ (r − Ra ) , dr v˜eff (r) Q ij

(34)

1a 1a (r) φaj , Dij = φai − 12 ∇2 + veff E D ˜ 1a = φã − 1 ∇2 + v˜1a (r) φã . D ij i eff j 2

(35)

= Ωa

(a)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 47

1a 1a In contrast to the one-center potentials veff (r) and v˜eff (r) (eqn. (31)), the potential v˜eff (r)

(eqn. (30)) is a function globally describing the total system. Note that v˜eff (r) appears in a ˆ ij (eqn. (34)). the formation of Fµν (eqn. (18)) as well as D

2.5

Energy expression

The SCF solution of eqn. (12) leads to the optimized total energy: E = E˜ +

X a

a (E 1a − E˜ 1a ) + U ({Ra }, {Zeff }) ,

(36)

a where the last term is the effective nuclear energy based on effective nuclear charges Zeff .

˜ E 1a and E˜ 1a are formulated as: The three energy contributions E, E˜ =

X

fµν Tµν + E H [˜ n+n ˆ+n ˜ Zc ] + E xc [˜ n+n ˜c] ,

(37)

µν

E 1a =

X

E˜ 1a =

X

xc 1a 1a a ρaij Tija + E H [n1a + n1a Zc ] + E [n + nc ] + Tcore ,

ij

(38) ρaij T˜ija

H

1a

a

+ E [˜ n +n ˆ +

n ˜ 1a Zc ]

xc

1a

+ E [˜ n +

n ˜ 1a c ],

ij

respectively, where Tµν , Tija and T˜ija are the kinetic integral matrices in basis of GTFs {χµ }, atomic all-electron partial waves {φai }, and pseudo partial waves {φãi }, respectively; and a Tcore is the core-electron kinetic energy, which is provided by atomic-data calculation. The

Hartree energy functional E H [n] is evaluated using v H [n] (eqn. (32)) as E H [n] = 21 hn|v H [n]i, and E xc [n] is the exchange-correlation (XC) energy. In this expression, the kinetic and the Hartree energies of the core density, which are constant for any geometries, are included, so that the resultant energies can be in some sense directly compared with those obtained by the conventional all-electron calculations. Note that the previous solid-state developments discarded them as the zero of energy.


Page 15 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3

Implementation

The details of the computer implementation of the GTF-PAW method is described below, particularly focusing on the integrals involving GTFs. In the present study, our DFT implementation is limited to the use of the local-density approximation with the restricted closed-shell treatment, while the main objective of the research is to initiate and validate the basic numerical scheme to incorporate the PAW method into the GTF-based DFT implementation. As shown below, various types of intermediate integrals arise in the GTF-PAW method, and are treated with different numerical techniques in an optimal and safe manner.

3.1

Use of atompaw and libpaw libraries

The numerical processes of the PAW method are divided into basis-independent and dependent parts. For the former part of the processes, the present scheme in part makes direct use of the previously-developed PAW implementations, which are available in the atompaw 32 and libpaw 29 libraries. Prior to molecular DFT calculations, the atompaw module is executed to generate atomic data of every element for PAW calculations, which are stored to disk files. The libpaw library provides a basic infrastructure to restore and manipulate the atomic data and intermediate arrays and to access utility functions. We exploit the existing subroutines of libpaw in order to perform the basis-independent a PAW processes, which are mainly for obtaining Dij (eqn. (33)) and E 1a − E˜ 1a (eqn. (38));

these two are evaluated using the subroutines named pawdij and pawdenpot, respectively. Given that the function values v˜eff (r) (eqn. (30)) at real-space grid points and the matrix elements ρaij (eqn. (24)) are determined using our GTF-based implementation, the evaluation a of Dij and E 1a − E˜ 1a itself is generally independent of the GTF basis scheme.

Tables 2 and 3 summarize the most relevant data accessed from our implementation, which are provided by atompaw and libpaw libraries.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 47

Table 2: Atomic data provided by atompaw and most relevant to our implementation. (a)

pñi li (r)

eqn. (5)

(a) gL (r) (a) n ˜ c (r) a Tcore a Zeff Zc va (r)

eqn. eqn. eqn. eqn. eqn.

(27) (28) (38) (44) (45)

Table 3: Intermediate data obtained using the subroutines of the libpaw library. saij ρaij n ˆ a (r) aLM qij a Dij P 1a − E ˜ 1a ) a (E

3.2 3.2.1

eqn. eqn. eqn. eqn. eqn. eqn.

(9) (24) (26) (27) (33) (38)

Analytic molecular integrals Kinetic integral

The kinetic integrals in Gauss-type AO basis 1 Tµν = hχµ | − ∇2 |χν i 2

(39)

can be analytically evaluated using the well-established recursive relation formula, which may be found in Ref. 11.

3.2.2

Coulomb integral matrix

The AO-basis Coulomb integral matrix associated with the pseudo density n ˜ (r) (eqn. (23)) is computed as: Jµν = hχµ |v H [˜ n](r)|χν i X = fλξ (µν|λξ) , λξ


(40)

Page 17 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


using the so-called electron repulsion integrals (ERIs): Z (µν|λξ) =

dr1 dr2

χµ (r1 )χν (r1 )χλ (r2 )χξ (r2 ) . |r1 − r2 |

(41)

This four-index integral can be efficiently computed in an analytical manner within the conventional Gauss-type implementation. Various approximate approaches to further reduce its expense have been developed, such as resolution-of-identity, 41 pseudospectral, 42 fast-multiple 43,44 methods, etc. In the plane-wave implementation, the Coulomb integral is evaluated with the grid-based approach solving the Poisson equation via the fast Fourier transformation (FFT).

3.2.3

Long-range screened nuclear integral

We describe a part of the one-electron Hamiltonian matrix using the long-range screened nuclear attraction integral, which can be analytically evaluated:

Xµν = hχµ |v LR (r)|χν i

(42)

with the smoothed nuclear potential:

v LR (r) =

X a

vaLR (r)

vaLR (|r − Ra |) ,

Za = − eff erf r

r √ 2rl

(43)

.

(44)

It takes essentially the same form as the finite nuclear potential based on Gauss-type distribution, 45 and analytically accounts for a slowly-decaying long-range component of the a nuclear potential (−Zeff /r for large r).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.3 3.3.1

Page 18 of 47

Single-center integrals Short-range screened nuclear integrals

The screened nuclear potential v H [˜ nZc ](r) is constructed by the superposition of the predetermined atomic radial function vaZc (r) as follows: v H [˜ nZc ](r) =

X

v H [˜ naZc ](r) =

X a

a

vaZc (|r − Ra |) .

(45)

The atomic data file provided by atompaw contains the tabulated radial values of vaZc (r) a as a function of radial grids for the core sphere r < rcore , whereas it smoothly merges with a −Zeff /r outside the sphere.

Given the long-range potential vaLR (r) (eqn. (44)), its complement to the full description of the screened potential vaZc (r) (eqn. (45)) is defined by:

vaSR (r) = vaZc (r) − vaLR (r) ,

(46)

which is adjustably local and smooth around each atom with a short-range tail. 29 This locality is associated with the fact that vaZc (r) and vaLR (r) have the same asymptotic Coulomb tail. The matrix form of the short-range potential in AO basis is given by:

Yµν =

X

a Yµν ,

(47)

a

a Yµν

Z =

dr vaSR (|r − Ra |)χµ (r)χν (r)

(48)

a Because of the rapid decay of vaSR (r), Yµν is evaluated as a finite-volume integral written


Page 19 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in spherical coordinates: a Yµν

r0

Z

dr r2 vaSR (r)

=

π

Z

0

2π

Z dθ sin θ

0

dϕ Sli mi (θ, ϕ) (49)

0

a

a

× χµ (r + R )χν (r + R ) . This single-center integral is evaluated using the Gauss-Chebyshev and Lebedev quadrature methods for the radial and angular components, respectively. The radius r0 in the integral is chosen such that |vaSR (r0 )| < τ with a small threshold τ . The Gauss-Chebyshev grids are 1+x transformed in a nonlinear fashion from [−1 : 1] to [0 : r0 ] with t(x) = r0 r0 −(r , assuming 0 −2)x

r0 > 2. This transformation allows for generating more grids in vicinity of the atom center. The values of vaZc (r) on the grids are derived using the spline interpolation implemented in libpaw from its tabulated values in the atomic data. The matrices Xµν and Yµν are obtained using different integral schemes, but the total of them results in the AO matrix of the screened nuclear potential as follows:

Xµν + Yµν = hχµ |v H [˜ nZc ](r)|χν i .

(50)

This range separation treatment for v H [˜ nZc ] was introduced in the wavelet-based PAW implementation. 29 As suggested in Ref. 29, we use 0.2 to 1.0 for rl (the range separation parameter in eqn. (44)); the resulting values of Xµν + Yµν are (or should be) insensitive to rl . Note that the evaluations of Xµν and Yµν can be carried out outside the SCF loop. 3.3.2

Projection of pãi on AO basis

The projection integral pãiµ (eqn. (15)) is written in spherical coordinates as: pãiµ

Z = 0

rpaw

(a) dr r pñi li (r)

π

Z

Z dθ sin θ

0

2π

dϕ Sli mi (θ, ϕ) 0

× χµ (r + Ra ) .


(51)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 47

This integral is evaluated with the single-center numerical quadrature using the GaussLegendre and Lebedev schemes for the radial and angular integrations, respectively. This is also done outside the SCF loop.

3.3.3

Hartree potential v H [ˆ na ]

In our implementation, the Hartree potential with the single-center compensation charge n â (eqn. (26)) is evaluated on given real-space grids rg :

H

Z

a

dr0

v [ˆ n ](rg ) = Ωa

n ˆ a (r0 ) . |rg − r0 |

(52)

It is transformed to the single-center integral form: v H [ˆ na ](rg ) =

X ijLM

4π (a) ρaij qijaLM GL (rga )SLM (θga , ϕag ) , 2L + 1

(53)

where (rga , θga , ϕag ) is the spherical coordinate representation of rag (= rg − Ra ). The radial (a)

integral GL (r) is expressed as: (a) GL (r)

Z

rshape

=

(a)

dr0 gL (r0 )

0

r0L+2 rL+1

(54)

a for r ≥ rshape , and

(a) GL (r)

=

r

r0L+2 (a) dr0 gL (r0 ) L+1 r 0 Z rshape rL (a) + dr0 gL (r0 ) 0L−1 r r

Z

(55)

a for 0 ≤ r < rshape . These radial integrals are evaluated with the Gauss-Ledgendre quadrature

method. Note that eqn. (53) is derived using the Laplace expansion of 1/|rag − r0 |. In the plane-wave implementation, this Hartree potential is evaluated using the FFTbased scheme instead of the above approach. 20 ACS Paragon Plus Environment

Page 21 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.3.4

Core-electron Hartree energy

In eqn. (37), the core energy corresponds to E H [˜ nZc ], which is decomposed into the atomic a a , respectively: and E˜2c one- and two-electron contributions, denoted E˜1c

E˜c ≡ E H [˜ nZc ] =

X a a ) + E˜2c (E˜1c

(56)

a

They are given as the single-center Coulomb integrals: Z

a a E˜1c = −Zeff

a E˜2c

1 = 2

dr

n ˜ ac (r) , |r − Ra |

˜ ac (r2 ) n ˜ ac (r1 ) n dr1 dr2 . |r1 − r2 |

Z

(57)

(58)

and are simplified to the following integrals in the one- and two-dimensional radial coordinates over the finite range:

a a E˜1c = −4πZeff

a E˜2c = 8π 2

Z

rcore

dr r n ˜ (a) c (r) ,

(59)

0

rcore

Z

dr1 n ˜ (a) ˜c(a) (r1 ) , c (r1 ) v

(60)

0

respectively, where v˜c(a) (r1 )

Z

r1

dr2 r22 n ˜ (a) c (r2 )

=r1 0

+

r12

Z

(61)

rcore

dr2 r2 n ˜ (a) c (r2 ) .

r1

In these forms, the core energy is computed using the Gauss-Ledgendre quadrature method. As noted earlier, it is independent of molecular geometric parameters, and can be precalculated for each element.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.4

Page 22 of 47

Atomic fine-mesh numerical integrals (a)

Related to the given form of the shape functions gL (r), the functions n ˆ a (r) (eqn. (26)) and (a)

a ). They are Qij (r − Ra ) (eqn. (27)) vanish outside the short-radius sphere (|r − Ra | > rshape

numerically regulated to a certain degree, but are not as smooth as pseudo structures built ˜ p and density n into pseudo orbitals Ψ ˜. Kresse et al. introduced the double grid technique to perform numerical quadratures for integrands that have different shape-complexity. In this technique, the integrands involving (a)

n ˆ a or Qij are evaluated using the fine-resolution FFT grids, which are spherically confined a around each atom. within the radius rshape

Accordingly, for the numerical quadrature of some types of single-center integrals, we use fine uniform mesh grids within the sphere volume, as defined by

a , Θa = rijk + Ra rijk = (i, j, k) × h where |rijk | < rshape

(62)

where i, j, and k are integers. The mesh size h is a user-input parameter to specify the resolution. The center of Θa for every atom is set to the associated atom center in our approach, while it is not the case in the FFT- or wavelet-based implementations in which the mesh merges into the FFT grids. Using the numerical quadrature with the fine-mesh grids Θa , we evaluate three types of integrals as described below. 3.4.1

Nonlocal contribution

a ˆ ij The quadrature for the nonlocal elements D (eqn. (34)) is done using libpaw according to

the following formula: X (a) â = 1 D v˜eff (rG ) Qij (rG − Ra ) , ij h3 G∈Θa


(63)

Page 23 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for which our GTF-based module evaluates the effective potential (eqn. (30)) on mesh grid points rG as: v˜eff (rG ) =v H [˜ n](rG ) +

X (v H [ˆ na ](rG ) + v H [˜ naZc ](rG )) (64)

a

+ v xc [˜ n(rG ) + n ˜ c (rG )] . With eqn. (23), the Hartree potential v H [˜ n](rG ) is computed as: v H [˜ n](rG ) =

X

fµν AG µν ,

(65)

χµ (r0 )χν (r0 ) |rG − r0 |

(66)

µν

AG µν

Z =

dr0

where the intermediate array AG µν is obtained using the extant Gaussian integral module based on the analytic formula, which is conventionally used for the analytic evaluation of the nuclear attraction integral. 11 The computation of v H [ˆ na ](rG ) is written earlier in the ˜ c (rG ) are section 3.3.3. The screened nuclear potential v H [˜ naZc ](rG ) and pseudo core density n obtained by the spline interpolation of the atomic data, as also indicated in the section 3.3.1.

3.4.2

Fock matrix contribution

The AO-basis Coulomb matrix associated with n ˆ is evaluated as: Rµν = hχµ |v H [ˆ n](r)|χν i =

1 X X a n ˆ (rG )AG µν , 3 h a G∈Θa

where n ˆ a (r) and AG µν are given in eqn. (26) and (66), respectively.


(67)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.4.3

Page 24 of 47

Energy contribution

The energy contribution associated with the compensation charge n ˆ to E˜ (eqn. (37)) is calculated as: E˜ 0 =

X

fµν Rµν +

µν

×

X 1 a0

2

1 X X a n ˆ (rG ) h3 a G∈Θa

H

a0

v [ˆ n ](rG ) + v

H

(68)

0 [˜ naZc ](rG )

It corresponds to Hartree energy represented as (ˆ n)(˜ n) + (ˆ n)(ˆ n)/2 + (ˆ n)(˜ nZc ) in the notation introduced in Ref. 22. The energy contributions other than E˜ 0 are computed using analytic methods and Becke’s quadrature scheme, as written later.

3.5

Numerical integrals using Becke grids

As done in the conventional Gaussian-based DFT implementation, we use Becke’s fuzzy cell xc and XC energy quadrature method for evaluating the XC potential matrix in AO basis Fµν

in eqn. (37). xc Fµν = hχµ |v xc [˜ n+n ˜ c ]|χν i X = ωg χµ (rg )χν (rg )v xc [˜ n(rg ) + n ˜ c (rg )] ,

(69)

g

E xc [˜ n+n ˜c] =

X

ωg εxc [˜ n(rg ) + n ˜ c (rg )] (˜ n(rg ) + n ˜ c (rg )) .

(70)

g

The pseudo valence density n ˜ and core density n ˜ c are evaluated on the Becke’s multi-center grids rg . In constructing Becke grids, we uses Lebedev method as the de facto standard angular quadrature scheme. Two different types of radial quadrature are used for our DFT calculations with and without the PAW method, respectively. In the both types, the original 24 ACS Paragon Plus Environment

Page 25 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


quadrature grids Rg and weights Wg generated in the interval [−1, 1] are transformed to the interval [0, ∞] with the logarithmic mapping, 35–37 1 2 rg = ln , ln 2 1 − Rg 1 Wg wgrad = . ln 2 1 − Rg

(71)

For the generation of Rg and Wg , the standard Gauss-Chebyshev (GC) quadrature of the second kind is employed for PAW-based calculations as the efficient scheme confirmed to perform the best in our numerical tests. However, for the non-PAW calculations, the modified GC quadrature of the second kind 37 is used because it works better than the standard GC in our test cases and is also used as default in other programs, such as orca. 8

3.6

Fock matrix and Hartree energy

Consolidating all the intermediate matrices described above, the Fock matrix Fµν (eqn. (18)) is obtained to be:

xc Fµν = Tµν + Xµν + Yµν + Jµν + Rµν + Fµν .

(72)

The total Hartree energy is calculated using E˜c (eqn. (56)), E˜ 0 (eqn. (68)), and 1 fµν Xµν + Yµν + Jµν 2

(73)

E H [˜ n+n ˆ+n ˜ Zc ] = E˜ 0 + E˜ 00 + E˜c .

(74)

E˜ 00 =

X µν

as given by

The energy components E˜ 00 and E˜c correspond to (˜ n)(˜ nZc ) + (˜ n)(˜ n)/2 and (˜ nZc )(˜ nZc )/2, respectively. 25 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.7

Page 26 of 47

Sketch of implementation

The overall program structure of our PAW-based DFT method is written as follows: 1. Initialize the libpaw module in which atomic data are restored from disk files. 2. Calculate pãiµ (eqn. (51)), Tµν (eqn. (39)), Xµν (eqn. (42)), Yµν (eqn. (47)), E˜c (eqn. (56)), ˜ µν (eqn. (17)). and S 3. Initialize MO coefficients Cµp . 4. Compute density matrices fµν (eqn. (22)) and ρaij (eqn. (24)) from Cµp . 5. Obtain

P

a (E

1a

− E˜ 1a ) (eqn. (38)) from ρaij using the libpaw module (pawdenpot).

6. With atomic fine-mesh grids {rG }: (a) Evaluate the following terms, • AG µν (eqn. (66))) • n ˜ (rG ) (eqn. (23)) • n ˆ a (rG ) (eqn. (26)) • v H [˜ n](rG ) (eqn. (65)) • v H [ˆ na ](rG ) (eqn. (53)) using fµν and ρaij , and obtain • n ˜ ac (rG ) • v H [˜ naZc ](rG ) directly from atomic data. They lead to the formation of v˜eff (rG ) (eqn. (64)). (b) Evaluate Rµν (eqn. (67)) and E˜ 0 (eqn. (68)) using {ˆ na (rG )} and {v H [ˆ na ](rG )}. a (c) Call the subroutine of the libpaw module (pawdij), which provides Dij (eqn. (33))

from {˜ veff (rG )} and ρaij . 26 ACS Paragon Plus Environment

Page 27 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


7. Compute Jµν (eqn. (40)) from fµν . 8. With Becke grids {rg }: xc Evaluate Fµν (eqn. (69)) and E xc [˜ n+n ˜ c ] (eqn. (70)).

9. Assemble the Fock matrix Fµν (eqn. (72)) and E˜ 00 (eqn. (73)). 10. Compute the total energy E (eqn. (36)). ˜ µν (eqn. (16)) and solve eqn. (12) to update Cµp and εp . 11. Construct H 12. Go to 4 unless the convergence criteria is met. Note that the most relevant steps added to the conventional GTF-based DFT code to incorporate the PAW method are the above steps 1, 2, 5, and 6.

4

Test calculations

The GTF-PAW method was tested by performing its LDA-level calculations on the closedshell atomic and molecular systems, the Mg atom, LiH, CrF6 , Cl2 , and [Cu(CN)2 ] – . They were all done with uncontracted (i.e., primitive) GTF basis. The SVWN5 functional was used for the LDA treatment in all the DFT calculations. In Table 4, the grid settings for the Becke multi-center quadrature (section 3.5) are summarized. Table 4: Grid levels of Becke multi-center quadrature used in the test calculations IntAcca Lebedev Mg LiH CrF6 Cl2 [Cu(CN)2 ] – 3.67 29 45 204 12 014 31 786 4.34 23 39 838 10 320 28 006 4.67 29 70 436 17 954 49 400 5.01 35 115 622 29 278 80 948 5.34 41 171 258 42 780 119 452 5.67 47 28 916 54 947 243 908 60 692 170 554 a The number of the radial grid points for the atom a is given as nrad (a) = −5(3 log − na + 8) with = 10−IntAcc where na corresponds to the row of the periodic table of a. Level A 3 4 5 6 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4.1

Mg atom

We tested two kinds of PAW-based frozen core (FC) treatments for the Mg atom, treating 1s and 1s2s2p as FC orbitals, respectively. The primitive s- and p-type GTFs of the cc-pVTZ basis set 46–48 was employed. In addition to the uncontracted GTFs (15s10p), we examined the truncated primitive sets using the exponents smaller than 20 and 2 a.u., resulting in (7s7p) and (5s5p), respectively, for the 1s and 1s2s2p FC treatments, respectively. Table 5 shows the resultant total energies along with the results obtained without PAW using the uncontracted and contracted GTFs. The grid level 7 was used as very accurate Becke grid (Table 4), and the mesh size h (eqn. (62)) was set to 0.15 bohr. With the 1s2s2p FC treatment, our PAW approach yielded the total energies with an error of 1 mEh relative to the atompaw result that is regarded as numerically exact. The basis truncation from (15s10p) to (5s5p) had only a minor impact thanks to the smooth nature of the pseudo orbitals. With the (15s10p) basis set, the total energy with PAW (199.139 05 Eh ) was lower than that without PAW (-199.138 08 Eh ), meaning that the former is better than the latter in the variational sense; this is because the core description provided by PAW atomic data was determined at basis-limit accuracy in atompaw. Table 5: Total energy (Eh ) of the Mg atom obtained using SVWN5 functional with and without PAW treatment. The GTF parameters of the cc-pVTZ basis functions were used for contracted and uncontracted s- and p-type GTFs. GTFs PAW Frozen Total energy (15s10p) with 1s -199.131 86 (7s7p)a with 1s -199.130 55 (15s10p) with 1s2s2p -199.139 05 (5s5p)b with 1s2s2p -199.138 85 (15s10p) without -199.138 08 (15s10p)→[5s4p] without -199.128 98 atompaw (exact) without -199.139 41 a cc-pVTZ’s primitive GTFs with the exponents α < 20.0. b cc-pVTZ’s primitive GTFs with the exponents α < 2.0. With the 1s FC approximation, the pseudo 2s, 2p, and 3s orbitals were determined as solution of the PAW-based KS equation. These pseudo AOs computed with (15s10p) 28 ACS Paragon Plus Environment

Page 28 of 47

Page 29 of 47

are displayed in Figure 1. They are transformed into the all-electron AOs via eqn. (1). For comparison, the all-electron and pseudo atomic orbitals determined by atompaw are included in Figure 1. The pseudo AOs represented with the GTFs show excellent agreement with the reference AOs obtained by atompaw, and indeed have reduced-nodal and ultrasoft structures. This holds more or less true for the pseudo AOs calculated with the truncated basis (7s7p), as shown in Figure S1. The errors of the total energies with the 1s FC treatment relative to the exact energy were much larger than those with the 1s2s2p FC treatment. This may be related to the fact that the small core treatment with 1s causes much larger numerical complexity in the compensation charge than the large core treatment with 1s2s2p. 2.0

2s

5

1.5 1.0

3

0.5

atomic orbital

4

2 1

1.4

2p

0.8

0.0 0.5

0.6 0.4 0.2

1.0

0.0

1

1.5

0.2

2

2.0

1

0 x / a.u.

1

2

exact AO (atompaw) pseudo AO (atompaw) transformed orbital (15s10p) pseudo orbital (15s10p)

1.0

0

2

3s

1.2

atomic orbital

6

atomic orbital

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2

1

0 x / a.u.

1

2

0.4

2

1

0 x / a.u.

1

2

Figure 1: One-dimensional plots of exact and pseudo AOs obtained from atomic calculations with atompaw, and the pseudo AO and its transformed AO via eqn. (1) determined using the GTF-PAW method with the (15s10p) GTF basis. These AOs were obtained with the a 1s frozen core approximation and SVWN5 functional. The augmentation region rpaw set to 1.707 a.u. is indicated by the gray region.

4.2

LiH without frozen core approximation

Even with no use of FC approximation, the PAW scheme imposes pseudization on orbitals. A single-point calculation was performed on the LiH molecule (Li-H = 3 bohr) without any orbitals treated as FCs. The atomic orbitals of the hydrogen atom anyway cannot be frozen. In the atomic data construction, the 1s, 2s are 2p orbitals of the Li atom were all treated as valence orbitals. The Becke grid level 7 was employed, and the mesh size h was set to 0.15 bohr. 29 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 6 shows the total energies of all-electron DFT calculations with and without the PAW scheme. The MOs were expanded into uncontracted STO-3G 49 and cc-pVQZ basis sets, written as (6s3p/3s) and (12s6p3d2f 1g/6s3p2d1f ), respectively, where (X/Y ) means that X and Y are basis sets for Li and H, respectively. The all-electron calculation of LiH yields two doubly-occupied MOs, denoted Ψi (i = 0, 1) for the i-th lowest orbital states. For comparison, the table includes the results obtained using bigdft 31 with the mesh size hgrid set to 0.5 a.u. Table 6: Total energy (in Eh ) of LiH molecule obtained using SVWN5 functional with and without all-electron PAW treatment. GTFs PAW Total energy bigdft (hgrid=0.5) with -7.916 19 (6s3p/3s)a with -7.913 46 a (6s3p/3s) without -7.821 14 STO-3G without -7.792 28 (12s6p3d2f 1g/6s3p2d1f )b with -7.919 30 b (12s6p3d2f 1g/6s3p2d1f ) without -7.919 46 cc-pVQZ without -7.918 93 a Uncontracted STO-3G. b Uncontracted cc-pVQZ. The uncontracted STO-3G basis set performed poorly when used in the conventional DFT approach. Its combination with the PAW showed a significant improvement in terms of the lowering of the total energy. Figure 2 illustrates that the cusp shapes of the MOs are greatly regularized in the pseudo MOs. This explains that the poor basis set with PAW gave rather good total energy to 6 mEh accuracy compared to the near-complete-basis cc-pVQZ energy. The uncontracted cc-pVQZ basis with and without the PAW treatment gave similar total energies with an error of 0.16 mEh . This showed desirable convergence of total energy with near complete basis in our PAW method.


Page 30 of 47

Page 31 of 47

3.0

0.4

0

2.5

all-electron MO (cc-pVQZ) all-electron MO (6s3p/3s) transformed MO (6s3p/3s) pseudo MO (6s3p/3s)

1

0.3 0.2 molecular orbital

2.0 molecular orbital

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.5 1.0 0.5

0.1 0.0

Li

H

0.1 0.2 0.3

0.0 0.5

Li 3

2

0.4

H 1

0 x / a.u.

1

2

3

0.5

3

2

1

0 x / a.u.

1

2

3

Figure 2: One-dimensional plots of the two doubly-occupied orbitals of LiH (Ψ0 and Ψ1 ) obtained with the conventional DFT method using the cc-pVQZ and uncontracted STO-3G (6s3p/3s) basis sets along with the corresponding transformed and pseudo orbitals resulting from the PAW-based DFT calculation with the uncontracted STO-3G (6s3p/3s) basis set.

4.3

CrF6

We calculated potential energy curves (PECs) of the octahedral CrF6 molecule (Figure 3 (a)) near equilibrium as a function of the Cr-F bond distance, commonly set to all the Cr-F bonds. One of our interests is to what extent the smooth character of pseudo wave functions serves as safely reducing the number of grid points of Becke’s multi-center numerical quadrature. For checking this, various grid levels shown in Table 4 were tested for DFT calculations with and without the PAW scheme. The grid settings of the levels 3 to 7 are taken from the preset grids defined in the orca program package. 8 The level-7 grid is the finest, assumed to provide numerically-convergent results that can be used as reference. The level A is derived from the level 4 by aggressively reducing its radial grid points by 64%. The cc-pVDZ and ccpVTZ basis sets were examined, but their uncontracted forms were used for the PAW-based calculations. Unless otherwise noted, the atomic fine-mesh size h used was 0.25 bohr. The PAW atomic data were constructed with 1s2s2p and 1s frozen for Cr and F, respectively, and with 3s3p3d4s4p and 2s2p treated as valence AOs of Cr and F, respectively. The PECs are plotted in Figures S2 and 4 for cc-pVDZ and cc-pVTZ, respectively. The errors in the total energies obtained with the level-X grid, denoted EX (R), are estimated relative to the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 47

level-7 total energies shifted by the difference of the minima of the PECs of the levels X and min 7, denoted EX and E7min , respectively, as follows:

min − E7min ] ∆E(R) = [EX (R) − E7 (R)] − [EX

(75)

With this error estimation, the parallelity of the PECs can be visually checked in the plots. The equilibrium bond lengths estimated by the fifth-order polynomial fitting are listed in Table 7. The curvatures of the PECs at the equilibrium corresponds to the harmonic frequencies of the PECs where we assumed the PECs to be the ones of a diatom with a reduced mass of 10 amu.

(a)

(b)

F F F

Cr F

F F

N C

Cu C N R

Cu-C = 1.8061 Å C≡N = 1.2278 Å

Figure 3: Molecular structures of (a) CrF6 and (b) [Cu(CN)2 ] – considered in this study. Overall, the PAW-based method provides numerically-convergent results with much lowerlevel grids than the conventional method. The striking difference arises between the PAWbased and conventional calculations when using the level-A grid. The pseudo MOs are smoothed and numerically lessened by the ultrasoft PAW treatment; therefore, the rather coarse radial grid of the level A performs well to describe them with high accuracy, while it causes large errors with the conventional method. The grid-level sensitivity was more or less similar for the cc-pVDZ and cc-pVTZ basis sets. In order to check the mesh-size dependence of PAW-based calculations, the PECs obtained using two kinds of atomic fine-mesh sizes, 0.15 and 0.25 bohr, were compared on the top of the level-7 grid (Figure 5). The relative differences in the total energies between the two meshes are less than 0.8 mEh and approximately constant across the PECs. Thus, the 32 ACS Paragon Plus Environment

Page 33 of 47

120.0

(a)

Level of Grids

A 3 4 5 6 7

-1637.905

∆E / 10−4 a.u.

E / a.u.

90.0 60.0

-1637.925 -1637.945 1.64

1.68

30.0

1.72 1.76 ˚ R(Cr-F) / A

1.80

0.0 -30.0 1.64 120.0

1.66

1.68

1.70

1.72

1.74

˚ R(Cr-F) / A

1.76

1.78

1.80

(b)

Level of Grids

A 3 4 5 6 7

-1637.710 E / a.u.

90.0

∆E / 10−4 a.u.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


60.0

-1637.730 -1637.750 1.64

1.68

30.0

1.72 1.76 ˚ R(Cr-F) / A

1.80

0.0 -30.0 1.64

1.66

1.68

1.70

1.72

1.74

˚ R(Cr-F) / A

1.76

1.78

1.80

Figure 4: Relative errors of the total energies defined by eqn. (75) for CrF6 calculated as a function of R(Cr-F) at the SVWN5 level of theory without PAW (upper panel) and with PAW (lower panel) using various levels of Becke multi-center quadrature grids, specified in Table 4. The PECs in absolute value are shown in each inset. The contracted and uncontracted basis of cc-pVTZ were used without and with the PAW method, respectively.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 7: Equilibrium bond lengths (Å) and curvatures ωe (cm−1 ) of CrF6 calculated from the PECs obtained at the SVWN5/cc-pVnZ level of theory (n = D and T) using various Becke grid levels with and without the PAW treatment. The details to obtain ωe are written in the main text. Level cc-pVDZ A 3 4 5 6 7 cc-pVTZ A 3 4 5 6 7

without PAWa R(Cr-F) ωe 1.7079 1.7237 1.7209 1.7211 1.7223 1.7216

2534 2272 2379 2375 2292 2343

with PAWb R(Cr-F) ωe 1.7178 1.7148 1.7176 1.7171 1.7172 1.7172

2426 2460 2435 2440 2439 2439

1.7079 2730 1.7228 1.7215 2314 1.7237 1.7171 2449 1.7222 1.7174 2420 1.7226 1.7178 2394 1.7225 1.7173 2418 1.7225 a Contracted basis. b Uncontracted basis.

2371 2371 2385 2380 2381 2381

chemical descriptions should not be so susceptible to the choice of h, which mainly affects local augmentation descriptions. The pseudization of MOs enables to remove high-exponent GTF components in the MO representations. We performed PAW calculations with three sets of the truncated primitive GTFs derived from cc-pVTZ. They are constructed by choosing the exponents α less than 100, 50, and 20 a.u., by which the number of primitive AOs were reduced from 383 to 334, 319, and 297, respectively. Figure 6 displays the PECs obtained with these calculations for the level-7 grid, confirming that the tested truncations had negligible impacts on the PEC shapes. The PECs for the other grid levels are shown in Figures S4.

4.4

Cl2

We next turn to the PEC calculations for the Cl2 molecule of near equilibrium geometries. The cc-pVDZ, cc-pVTZ, and cc-pVQZ basis sets were employed; their uncontracted GTFs


Page 34 of 47

Page 35 of 47

25.0 -1637.720

Level 7 (h = 0.25 bohr) Level 7 (h = 0.15 bohr)

∆ E / 10−4 a.u.

E / a.u.

20.0

-1637.730

-1637.740

15.0

1.64

1.68

1.72

1.76

1.80

˚ R(Cr-F) / A

10.0

5.0

0.0 1.64

1.66

1.68

1.70

1.72

˚ R(Cr-F) / A

1.74

1.76

1.78

1.80

Figure 5: Relative errors of the PEC energies of CrF6 calculated at the SVWN5 functional level using the uncontracted cc-pVTZ basis and the level-7 grid with the atomic fine mesh size h set to 0.25 bohr with respect to those with h set to 0.15 bohr. The PECs in absolute value are shown in the inset.

9.0

Basis Set VTZ VTZ (α < 100) VTZ (α < 50) VTZ (α < 20)

6.0

-1637.700

E / a.u.

12.0

∆E / 10−5 a.u.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-1637.720

-1637.740

3.0

1.64

1.68

1.72 1.76 ˚ R(Cr-F) / A

1.80

0.0 -3.0 -6.0 1.64

1.66

1.68

1.70

1.72

˚ R(Cr-F) / A

1.74

1.76

1.78

1.80

Figure 6: Relative errors of the total energies defined by eqn. (75) for CrF6 calculated at the SVWN5 functional level with the uncontracted cc-pVTZ basis and its truncated sets that use exponents less than α. The PECs in absolute value are shown in the inset.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 47

were used for the conventional and PAW-based calculations. The mesh size h was set to 0.25 bohr. Atomic data was determined with 1s2s2p and 3s3p treated as FC and valence orbitals of the Cl atom, respectively. The PECs obtained with cc-pVDZ, cc-pVTZ, and cc-pVQZ were shown in Figures S5, S6, and 7, respectively. As discussed earlier, the pseudization greatly mitigated the gridlevel dependence of the total energies as well as the PEC shapes. Table 8 presents the bond lengths and total energies of equilibrium, and harmonic frequencies calculated from the PECs with the fifth-order polynomial fitting. With the PAW scheme, the bond lengths and harmonic frequencies were accurately predicted at the grid level 3 with an error of 0.0001 Å and 1 cm−1 , respectively. Without the PAW treatment, they exhibited slow and somewhat oscillating convergence with increasing grid level. With the high-level Becke grids, the PAW method provided the lowest total energy for the uncontracted cc-pVDZ and ccpVTZ basis sets, and nearly converged to the results obtained by the conventional method for the uncontracted cc-pVQZ basis set. Table 8: Equilibrium bond lengths (Å), harmonics frequencies ωe (cm−1 ), and minima of total energies Ee (a.u.) of Cl2 calculated from the PECs obtained at the SVWN5/cc-pVnZ level of theory (n = D, T, and Q) using various Becke grid levels with and without the PAW treatment. Level cc-pVDZ A 3 4 5 6 7 cc-pVTZ A 3 4 5 6 7 cc-pVQZ A 3 4 5 6 7

without PAWa R(Cl-Cl) ωe Ee

without PAWb R(Cl-Cl) ωe Ee

R(Cl-Cl)

with PAWb ωe

Ee

2.0181 2.0194 2.0199 2.0200 2.0198 2.0198

551 549 541 543 542 545

-917.40661 -917.39807 -917.39827 -917.39829 -917.39829 -917.39830

2.0209 2.0222 2.0227 2.0228 2.0226 2.0225

549 547 540 542 540 543

-917.43376 -917.42528 -917.42547 -917.42549 -917.42549 -917.42549

2.0082 2.0081 2.0082 2.0082 2.0082 2.0082

539 540 540 540 540 540

-917.44350 -917.44349 -917.44349 -917.44349 -917.44349 -917.44349

1.9930 1.9944 1.9939 1.9942 1.9939 1.9942

568 563 556 557 558 559

-917.44902 -917.44438 -917.44474 -917.44478 -917.44478 -917.44478

1.9937 1.9950 1.9945 1.9948 1.9945 1.9948

569 565 557 558 560 560

-917.46570 -917.46107 -917.46142 -917.46146 -917.46147 -917.46147

2.0021 2.0019 2.0019 2.0019 2.0019 2.0019

560 561 561 561 561 561

-917.46382 -917.46386 -917.46384 -917.46385 -917.46384 -917.46384

1.9870 1.9882 1.9877 1.9879 1.9877 1.9880

570 564 557 558 559 560

-917.45348 -917.45567 -917.45525 -917.45520 -917.45520 -917.45520

1.9874 570 -917.46850 1.9886 564 -917.47069 1.9880 557 -917.47027 1.9883 558 -917.47022 1.9881 559 -917.47022 1.9883 560 -917.47022 a Contracted basis. b Uncontracted basis.

1.9883 1.9884 1.9883 1.9883 1.9883 1.9883

563 563 563 563 563 563

-917.46991 -917.46991 -917.46990 -917.46991 -917.46990 -917.46990


Page 37 of 47

54.0

(a)

Level of Grids

A 3 4 5 6 7

E / a.u.

-917.452

∆E / 10−6 a.u.

36.0

-917.454

-917.456

18.0

1.94

1.98 2.02 ˚ R(Cl-Cl) / A

2.06

0.0

-18.0 1.94 54.0

1.96

1.98

2.00

˚ R(Cl-Cl) / A

2.02

2.04

2.06

(b)

Level of Grids

A 3 4 5 6 7

E / a.u.

-917.468

36.0

∆E / 10−6 a.u.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-917.470

-917.472

18.0

1.94

1.98 2.02 ˚ R(Cl-Cl) / A

2.06

0.0

-18.0 1.94

1.96

1.98

2.00

˚ R(Cl-Cl) / A

2.02

2.04

2.06

Figure 7: Relative errors of the total energies defined by eqn. (75) for Cl2 calculated as a function of the bond length at the SVWN5 level of theory without PAW (upper panel) and with PAW (lower panel) using various levels of Becke multi-center quadrature grids, specified in Table 4. The PECs in absolute value are shown in each inset. The contracted and uncontracted basis of cc-pVQZ were used without and with the PAW method, respectively.



Figure 8 shows the bond lengths predicted with the grid level 7 as a function of the basis set. With the cc-pVDZ basis set, the PAW scheme predicted the bond length with the smallest error relative to the basis-set-limit prediction. However, as the basis level was increased to cc-pVTZ, the PAW prediction improved over the cc-pVDZ counterpart but turned to be inferior to the results obtained by the conventional method. This error may be possibly ameliorated by reoptimizing the primitive exponents under the presence of the froze core approximation. 2.025 conventional a conventional b PAW b

2.020

˚ bond length / A

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.015 2.010 2.005 2.000 1.995 1.990 1.985

cc-pVDZ

cc-pVTZ

cc-pVQZ

basis set

Figure 8: Plots of equilibrium bond lengths as a function of the basis set levels obtained at the level-7 grid with and without the PAW treatment, as shown in Table 8.

4.5

[Cu(CN)2 ] –

Finally we tested [Cu(CN)2 ] – molecule and its dissociation limit CN – + CuCN. Here the geometries of CuCN and CN – were fixed in the calculations. The bond lengths R(Cu-C) and R(C-N) in these fixed fragments were set to 1.8061 and 1.2278 Å, respectively (Figure 3 (b)). The bonding energy De is defined as the energy difference between the anion complex [Cu(CN)]2 – and the two fragments CuCN and CN – . Uncontracted cc-pVDZ and cc-pVTZ basis sets were employed in the PAW calculations, where their contracted GTFs were used in the conventional DFT calculations. The atomic mesh size was set to 0.25 bohr as in the previous cases. The atomic data were constructed with 1s2s2p, 1s, and 1s electrons as frozen core for Cu, C, and N, respectively, and with 3s3p3d4s4p, 2s2p, and 2s2p electrons as valence electrons for Cu, C, and N, respectively. As listed in Table 9, the equilibrium bond 38 ACS Paragon Plus Environment

Page 38 of 47

Page 39 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


lengths R(Cu-N) are estimated by the fifth-order polynomial fitting. The curvatures of the PECs at the equilibrium corresponds to the harmonic frequencies of the diatomic PECs with a reduced mass of 1 amu. The the PAW scheme predicted the equilibrium bond lengths and harmonic frequencies with an error less than 0.0002 Å and 1 cm−1 at the level-3 grid with cc-pVDZ basis sets, compared to the results obtained with the most accurate level-7 grid. Even the rather coarse level-A grid predicted very accurate results. Conventional DFT calculations without PAW treatment show the similar features of oscillating convergence with increasing grid level as shown in Cl2 and CrF6 . Compared to cc-pVDZ basis sets, the stability of the PAW scheme with the grid setting has been enhanced when cc-pVTZ basis set was used in this case, while the conventional DFT still showed oscillating convergence with increasing grid level. From Table 9, we can see that the dissociation energy (De ) is less sensitive to the level of grid in both PAW scheme and conventional DFT, and the difference of calculated De between the two schemes are less than 0.1 eV. Table 9: Equilibrium bond lengths (Å), curvatures ωe (cm−1 ) of the [Cu(CN)2 ] – complex calculated at the SVWN5/cc-pVnZ level with and without PAW using various levels of Becke grids, along with the dissociation energies De (eV) corresponding to the reaction [Cu(CN)2 ] – → CN – + CuCN. Level cc-pVDZ A 3 4 5 6 7 cc-pVTZ A 3 4 5 6 7

without PAWa Re ωe De

with PAWb Re ωe De

1.8227 1.8264 1.8254 1.8257 1.8258 1.8260

2128 2096 2075 2068 2078 2067

5.129 5.126 5.127 5.127 5.127 5.127

1.8306 1.8304 1.8307 1.8306 1.8306 1.8306

2079 2077 2077 2077 2077 2077

5.024 5.024 5.024 5.024 5.024 5.024

1.8247 1.8278 1.8265 1.8268 1.8270 1.8271

2093 2090 2062 2058 2073 2059

4.706 4.706 4.705 4.706 4.706 4.706

1.8310 1.8309 1.8310 1.8309 1.8310 1.8310

2064 2064 2063 2063 2063 2063

4.655 4.656 4.655 4.655 4.655 4.655



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

5

Conclusions

We have presented a basic framework to incorporate the PAW method into the conventional GTF-based DFT solver in a straightforward manner. The formulas and techniques necessary to implement the present scheme have been described in detail. A number of additional types of intermediate integrals need to be implemented in the existing DFT module. Thanks to the access to the previously-developed PAW libraries, atompaw and libpaw, which provide fundamental infrastructures, our development mainly focused on the complementary implementations related to GTF basis. The illustrative molecular calculations showed that the PAW scheme and inherent ultrasoft treatment greatly regularize or smooth the cusp shapes of molecular orbitals represented with GTFs. This mainly benefits the efficiency of the DFT calculation in the sense that the primitive GTFs with high exponents can be truncated, and a much smaller number of Becke multi-center grid points deliver high accuracy and stability to integrate exchangecorrelation functionals than used in the conventional implementation. These features indeed underlie the great success of the plane-wave PAW implementations. Although the ECP for QC calculations is a well-established pseudopotential method, the PAW method has several appealing merits; particularly, well controlled connectivity to all-electron description. The PAW-related computation is readily divided into atom-wise tasks associated with local augmentation regions, so that its linear scaling implementation can be easily done. The construction of PAW atomic data is rather straightforward and automatic, and our implementation can directly use extant atomic data in abinit file format via libpaw. As illustrated in the test cases, the PAW scheme has transparency to all-electron calculations; for example, total energies are in some sense comparable. Thus, the errors in the PAW approximation can be controllably estimated in a relatively simple manner. Our test cases showed that when using the uncontracted basis at near-complete-basis cc-pVQZ level, the total energies with and without the PAW method shows good agreement with an error less than 1 mEh . This shows that the PAW method has desirable basis-set convergence behavior 40 ACS Paragon Plus Environment

Page 40 of 47

Page 41 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


to the all-electron description. There are a number of directions to extend the current study. The present implementation is limited to the LDA-level functional, while the incorporation of more sophisticated generalized gradient approximation and hybrid functionals is required for practical calculations. The energy gradients with respect to nuclear coordinates should be implemented, as routinely done in the plane-wave PAW method. In this work, the primitive GTFs of the existing contracted basis sets were used for the basis expansion of MOs; however, this is obviously inefficient. The introduction of contracted basis and the reoptimization of exponents should be investigated.

Acknowledgement This work was supported partially by Natural Science Foundation of China (21501189) and Grant-in-Aid for Invitational Fellowship for Research in Japan (L16525) for X.-G. X., and Grant-in-Aid for Scientific Research (B) and Scientific Research on Innovative Areas “Photosynergetics,” with JSPS KAKENHI Grant Numbers JP16H04101 and JP26107010, respectively, for T. Y.

References (1) Becke, A. D. Perspective: Fifty years of density-functional theory in chemical physics. J. Chem. Phys. 2014, 140, 18A301. (2) Becke, A. D. A multicenter numerical integration scheme for polyatomic molecules. J. Chem. Phys. 1988, 88, 2547–2553. (3) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098–3100.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(4) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648–5652. (5) Pople, J. A.; Gill, P. M.; Johnson, B. G. Kohn-Sham density-functional theory within a finite basis set. Chem. Phys. Lett. 1992, 199, 557–560. (6) Gill, P. M.; Johnson, B. G.; Pople, J. A. A standard grid for density functional calculations. Chem. Phys. Lett. 1993, 209, 506–512. (7) Johnson, B. G.; Gill, P. M. W.; Pople, J. A. The performance of a family of density functional methods. J. Chem. Phys. 1993, 98, 5612–5626. (8) Neese, F. The ORCA program system. WIREs: Comp. Mol. Sci. 2012, 2, 73–78. (9) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: a generalpurpose quantum chemistry program package. WIREs: Comp. Mol. Sci. 2012, 2, 242– 253. (10) Kendall, R. A.; Aprà, E.; Bernholdt, D. E.; Bylaska, E. J.; Dupuis, M.; Fann, G. I.; Harrison, R. J.; Ju, J.; Nichols, J. A.; Nieplocha, J.; Straatsma, T.; Windus, T. L.; Wong, A. T. High performance computational chemistry: An overview of NWChem a distributed parallel application. Comp. Phys. Comm. 2000, 128, 260–283. (11) Helgaker, T.; Jorgensen, P.; Olsen, J. Molecular Electronic-Structure Theory; John Wiley & Sons, 2014. (12) Martin, R. M. Electronic Structure: Basic Theory and Practical Methods; Cambridge Univ. Press: Cambridge, 2004. (13) Schwarz, K.; Blaha, P. Solid state calculations using WIEN2k. Comp. Mat. Sci. 2003, 28, 259–273. (14) Lippert, G.; Hutter, J.; Parrinello, M. A hybrid Gaussian and plane wave density functional scheme. Mol. Phys. 1997, 92, 477–488. 42 ACS Paragon Plus Environment

Page 42 of 47

Page 43 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(15) Lippert, G.; Hutter, J.; Parrinello, M. The Gaussian and augmented-plane-wave density functional method for ab initio molecular dynamics simulations. Theo. Chem. Acc. 1999, 103, 124–140. (16) Krauss, M.; Stevens, W. J. Effective Potentials in Molecular Quantum Chemistry. Ann. Rev. Phys. Chem. 1984, 35, 357–385. (17) Nicklass, A.; Dolg, M.; Stoll, H.; Preuss, H. Ab initio energy-adjusted pseudopotentials for the noble gases Ne through Xe: Calculation of atomic dipole and quadrupole polarizabilities. J. Chem. Phys. 1995, 102, 8942–8952. (18) Dolg, M.; Cao, X. Relativistic Pseudopotentials: Their Development and Scope of Applications. Chem. Rev. 2012, 112, 403–480, PMID: 21913696. (19) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 1994, 50, 17953–17979. (20) Vanderbilt, D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys. Rev. B 1990, 41, 7892–7895. (21) Blöchl, P. E. Generalized separable potentials for electronic-structure calculations. Phys. Rev. B 1990, 41, 5414–5416. (22) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmentedwave method. Phys. Rev. B 1999, 59, 1758–1775. (23) Gonze, X.; Amadon, B.; Anglade, P.-M.; Beuken, J.-M.; Bottin, F.; Boulanger, P.; Bruneval, F.; Caliste, D.; Caracas, R.; Côté, M.; Deutsch, T.; Genovese, L.; Ghosez, P.; Giantomassi, M.; Goedecker, S.; Hamann, D.; Hermet, P.; Jollet, F.; Jomard, G.; Leroux, S.; Mancini, M.; Mazevet, S.; Oliveira, M.; Onida, G.; Pouillon, Y.; Rangel, T.; Rignanese, G.-M.; Sangalli, D.; Shaltaf, R.; Torrent, M.; Verstraete, M.; Zerah, G.; Zwanziger, J. ABINIT: First-principles approach to material and nanosystem properties. Comp. Phys. Comm. 2009, 180, 2582–2615. 43 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(24) 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.; Martin-Samos, 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.: Cond. Matt. 2009, 21, 395502. (25) Tackett, A.; Holzwarth, N.; Matthews, G. A Projector Augmented Wave (PAW) code for electronic structure calculations, Part II: pwpaw for periodic solids in a plane wave basis. Comp. Phys. Comm. 2001, 135, 348–376. (26) Briggs, E. L.; Sullivan, D. J.; Bernholc, J. Real-space multigrid-based approach to large-scale electronic structure calculations. Phys. Rev. B 1996, 54, 14362–14375. (27) Larsen, A. H.; Vanin, M.; Mortensen, J. J.; Thygesen, K. S.; Jacobsen, K. W. Localized atomic basis set in the projector augmented wave method. Phys. Rev. B 2009, 80, 195112. (28) Hine, N. D. M. Linear-scaling density functional theory using the projector augmented wave method. J. Phys.: Cond. Matt. 2017, 29, 024001. (29) Rangel, T.; Caliste, D.; Genovese, L.; Torrent, M. A wavelet-based Projector Augmented-Wave (PAW) method: Reaching frozen-core all-electron precision with a systematic, adaptive and localized wavelet basis set. Comp. Phys. Comm. 2016, 208, 1–8. (30) Torrent, M.; Jollet, F.; Bottin, F.; Zérah, G.; Gonze, X. Implementation of the projector augmented-wave method in the ABINIT code: Application to the study of iron under pressure. Comp. Mat. Sci. 2008, 42, 337–351. 44 ACS Paragon Plus Environment

Page 44 of 47

Page 45 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(31) Mohr, S.; Ratcliff, L. E.; Boulanger, P.; Genovese, L.; Caliste, D.; Deutsch, T.; Goedecker, S. Daubechies wavelets for linear scaling density functional theory. J. Chem. Phys. 2014, 140, 204110. (32) Holzwarth, N.; Tackett, A.; Matthews, G. A Projector Augmented Wave (PAW) code for electronic structure calculations, Part I: atompaw for generating atom-centered functions. Comp. Phys. Comm. 2001, 135, 329–347. (33) Hamann, D. R.; Schlüter, M.; Chiang, C. Norm-Conserving Pseudopotentials. Phys. Rev. Lett. 1979, 43, 1494–1497. (34) Booth, G. H.; Tsatsoulis, T.; Chan, G. K.-L.; Grüneis, A. From plane waves to local Gaussians for the simulation of correlated periodic systems. J. Chem. Phys. 2016, 145, 084111. (35) Treutler, O.; Ahlrichs, R. Efficient molecular numerical integration schemes. J. Chem. Phys. 1995, 102, 346–354. (36) Mura, M. E.; Knowles, P. J. Improved radial grids for quadrature in molecular densityâĂŘfunctional calculations. J. Chem. Phys. 1996, 104, 9848–9858. (37) Krack, M.; Köster, A. M. An adaptive numerical integrator for molecular integrals. J. Chem. Phys. 1998, 108, 3226–3234. (38) Lebedev, V.; Laikov, D. A quadrature formula for the sphere of the 131st algebraic order of accuracy. Doklady. Mathematics. 1999; pp 477–481. (39) Torrent, M.; Holzwarth, N.; Jollet, F.; Harris, D.; Lepley, N.; Xu, X. Electronic structure packages: Two implementations of the projector augmented wave (PAW) formalism. Comp. Phys. Comm. 2010, 181, 1862–1867. (40) Holzwarth, N. A. W.; Matthews, G. E.; Dunning, R. B.; Tackett, A. R.; Zeng, Y. Comparison of the projector augmented-wave, pseudopotential, and linearized augmented45 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

plane-wave formalisms for density-functional calculations of solids. Phys. Rev. B 1997, 55, 2005–2017. (41) Vahtras, O.; Almlöf, J.; Feyereisen, M. Integral approximations for LCAO-SCF calculations. Chem. Phys. Lett. 1993, 213, 514–518. (42) Murphy, R. B.; Cao, Y.; Beachy, M. D.; Ringnalda, M. N.; Friesner, R. A. Efficient pseudospectral methods for density functional calculations. J. Chem. Phys. 2000, 112, 10131–10141. (43) White, C. A.; Johnson, B. G.; Gill, P. M.; Head-Gordon, M. The continuous fast multipole method. Chem. Phys. Lett. 1994, 230, 8–16. (44) Strain, M. C.; Scuseria, G. E.; Frisch, M. J. Achieving Linear Scaling for the Electronic Quantum Coulomb Problem. Science 1996, 271, 51–53. (45) Visscher, L.; Dyall, K. Dirac-Fock atomic electronic structure calculations using different nuclear charge distributions. At. Data Nucl. Data Tables 1997, 67, 207–224. (46) Dunning Jr, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (47) Dunning Jr, T. H.; Woon, D. E. Gaussian Basis Sets for Use in Correlated Molecular Calculations. III. The Second Row Atoms, Al- Ar. J. Chem. Phys 1993, 98, 1358–1371. (48) Balabanov, N. B.; Peterson, K. A. Systematically convergent basis sets for transition metals. I. All-electron correlation consistent basis sets for the 3 d elements Sc–Zn. J. Chem. Phys 2005, 123, 064107. (49) Hehre, W. J.; Stewart, R. F.; Pople, J. A. self-consistent molecular-orbital methods. i. use of gaussian expansions of Slater-type atomic orbitals. J. Chem. Phys. 1969, 51, 2657–2664.


Page 46 of 47

Page 47 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry

Projector Augmented Wave (PAW) method enables smoothing molecular orbitals represented with Gauss-type basis function.


Projector Augmented Wave Method Incorporated into Gauss-Type

Recommend Documents